diff options
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials')
-rw-r--r-- | packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs | 181 |
1 files changed, 129 insertions, 52 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs index 1460680..e5a2e4d 100644 --- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs +++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs | |||
@@ -5,6 +5,12 @@ | |||
5 | {-# LANGUAGE MultiWayIf #-} | 5 | {-# LANGUAGE MultiWayIf #-} |
6 | {-# LANGUAGE OverloadedStrings #-} | 6 | {-# LANGUAGE OverloadedStrings #-} |
7 | {-# LANGUAGE ScopedTypeVariables #-} | 7 | {-# LANGUAGE ScopedTypeVariables #-} |
8 | {-# LANGUAGE DeriveGeneric #-} | ||
9 | {-# LANGUAGE TypeOperators #-} | ||
10 | {-# LANGUAGE KindSignatures #-} | ||
11 | {-# LANGUAGE TypeSynonymInstances #-} | ||
12 | {-# LANGUAGE FlexibleInstances #-} | ||
13 | {-# LANGUAGE FlexibleContexts #-} | ||
8 | 14 | ||
9 | ----------------------------------------------------------------------------- | 15 | ----------------------------------------------------------------------------- |
10 | -- | | 16 | -- | |
@@ -130,6 +136,7 @@ import qualified Data.Vector.Storable.Mutable as VM | |||
130 | 136 | ||
131 | import Data.Coerce (coerce) | 137 | import Data.Coerce (coerce) |
132 | import System.IO.Unsafe (unsafePerformIO) | 138 | import System.IO.Unsafe (unsafePerformIO) |
139 | import GHC.Generics | ||
133 | 140 | ||
134 | import Numeric.LinearAlgebra.Devel (createVector) | 141 | import Numeric.LinearAlgebra.Devel (createVector) |
135 | 142 | ||
@@ -244,64 +251,128 @@ data ODEMethod = SDIRK_2_1_2 Jacobian | |||
244 | | KVAERNO_7_4_5' | 251 | | KVAERNO_7_4_5' |
245 | | ARK548L2SA_DIRK_8_4_5 Jacobian | 252 | | ARK548L2SA_DIRK_8_4_5 Jacobian |
246 | | ARK548L2SA_DIRK_8_4_5' | 253 | | ARK548L2SA_DIRK_8_4_5' |
254 | | HEUN_EULER_2_1_2 Jacobian | ||
255 | | HEUN_EULER_2_1_2' | ||
256 | | BOGACKI_SHAMPINE_4_2_3 Jacobian | ||
257 | | BOGACKI_SHAMPINE_4_2_3' | ||
258 | | ARK324L2SA_ERK_4_2_3 Jacobian | ||
259 | | ARK324L2SA_ERK_4_2_3' | ||
260 | | ZONNEVELD_5_3_4 Jacobian | ||
261 | | ZONNEVELD_5_3_4' | ||
262 | | ARK436L2SA_ERK_6_3_4 Jacobian | ||
263 | | ARK436L2SA_ERK_6_3_4' | ||
264 | | SAYFY_ABURUB_6_3_4 Jacobian | ||
265 | | SAYFY_ABURUB_6_3_4' | ||
266 | | CASH_KARP_6_4_5 Jacobian | ||
267 | | CASH_KARP_6_4_5' | ||
247 | | FEHLBERG_6_4_5 Jacobian | 268 | | FEHLBERG_6_4_5 Jacobian |
248 | | FEHLBERG_6_4_5' | 269 | | FEHLBERG_6_4_5' |
270 | | DORMAND_PRINCE_7_4_5 Jacobian | ||
271 | | DORMAND_PRINCE_7_4_5' | ||
272 | | ARK548L2SA_ERK_8_4_5 Jacobian | ||
273 | | ARK548L2SA_ERK_8_4_5' | ||
274 | | VERNER_8_5_6 Jacobian | ||
275 | | VERNER_8_5_6' | ||
276 | | FEHLBERG_13_7_8 Jacobian | ||
277 | | FEHLBERG_13_7_8' | ||
278 | deriving Generic | ||
249 | 279 | ||
280 | constrName :: (HasConstructor (Rep a), Generic a)=> a -> String | ||
281 | constrName = genericConstrName . from | ||
282 | |||
283 | class HasConstructor (f :: * -> *) where | ||
284 | genericConstrName :: f x -> String | ||
285 | |||
286 | instance HasConstructor f => HasConstructor (D1 c f) where | ||
287 | genericConstrName (M1 x) = genericConstrName x | ||
288 | |||
289 | instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where | ||
290 | genericConstrName (L1 l) = genericConstrName l | ||
291 | genericConstrName (R1 r) = genericConstrName r | ||
292 | |||
293 | instance Constructor c => HasConstructor (C1 c f) where | ||
294 | genericConstrName x = conName x | ||
295 | |||
296 | instance Show ODEMethod where | ||
297 | show x = constrName x | ||
298 | |||
299 | -- FIXME: We can probably do better here with generics | ||
250 | getMethod :: ODEMethod -> Int | 300 | getMethod :: ODEMethod -> Int |
251 | getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2 | 301 | getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2 |
252 | getMethod (SDIRK_2_1_2') = sDIRK_2_1_2 | 302 | getMethod (SDIRK_2_1_2') = sDIRK_2_1_2 |
253 | getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2 | 303 | getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2 |
254 | getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2 | 304 | getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2 |
255 | getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2 | 305 | getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2 |
256 | getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2 | 306 | getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2 |
257 | getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3 | 307 | getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3 |
258 | getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3 | 308 | getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3 |
259 | getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3 | 309 | getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3 |
260 | getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3 | 310 | getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3 |
261 | getMethod (CASH_5_2_4 _) = cASH_5_2_4 | 311 | getMethod (CASH_5_2_4 _) = cASH_5_2_4 |
262 | getMethod (CASH_5_2_4') = cASH_5_2_4 | 312 | getMethod (CASH_5_2_4') = cASH_5_2_4 |
263 | getMethod (CASH_5_3_4 _) = cASH_5_3_4 | 313 | getMethod (CASH_5_3_4 _) = cASH_5_3_4 |
264 | getMethod (CASH_5_3_4') = cASH_5_3_4 | 314 | getMethod (CASH_5_3_4') = cASH_5_3_4 |
265 | getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4 | 315 | getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4 |
266 | getMethod (SDIRK_5_3_4') = sDIRK_5_3_4 | 316 | getMethod (SDIRK_5_3_4') = sDIRK_5_3_4 |
267 | getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4 | 317 | getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4 |
268 | getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4 | 318 | getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4 |
269 | getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4 | 319 | getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4 |
270 | getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4 | 320 | getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4 |
271 | getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5 | 321 | getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5 |
272 | getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5 | 322 | getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5 |
273 | getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5 | 323 | getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5 |
274 | getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5 | 324 | getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5 |
275 | getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5 | 325 | getMethod (HEUN_EULER_2_1_2 _) = hEUN_EULER_2_1_2 |
276 | getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5 | 326 | getMethod (HEUN_EULER_2_1_2') = hEUN_EULER_2_1_2 |
327 | getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3 | ||
328 | getMethod (BOGACKI_SHAMPINE_4_2_3') = bOGACKI_SHAMPINE_4_2_3 | ||
329 | getMethod (ARK324L2SA_ERK_4_2_3 _) = aRK324L2SA_ERK_4_2_3 | ||
330 | getMethod (ARK324L2SA_ERK_4_2_3') = aRK324L2SA_ERK_4_2_3 | ||
331 | getMethod (ZONNEVELD_5_3_4 _) = zONNEVELD_5_3_4 | ||
332 | getMethod (ZONNEVELD_5_3_4') = zONNEVELD_5_3_4 | ||
333 | getMethod (ARK436L2SA_ERK_6_3_4 _) = aRK436L2SA_ERK_6_3_4 | ||
334 | getMethod (ARK436L2SA_ERK_6_3_4') = aRK436L2SA_ERK_6_3_4 | ||
335 | getMethod (SAYFY_ABURUB_6_3_4 _) = sAYFY_ABURUB_6_3_4 | ||
336 | getMethod (SAYFY_ABURUB_6_3_4') = sAYFY_ABURUB_6_3_4 | ||
337 | getMethod (CASH_KARP_6_4_5 _) = cASH_KARP_6_4_5 | ||
338 | getMethod (CASH_KARP_6_4_5') = cASH_KARP_6_4_5 | ||
339 | getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5 | ||
340 | getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5 | ||
341 | getMethod (DORMAND_PRINCE_7_4_5 _) = dORMAND_PRINCE_7_4_5 | ||
342 | getMethod (DORMAND_PRINCE_7_4_5') = dORMAND_PRINCE_7_4_5 | ||
343 | getMethod (ARK548L2SA_ERK_8_4_5 _) = aRK548L2SA_ERK_8_4_5 | ||
344 | getMethod (ARK548L2SA_ERK_8_4_5') = aRK548L2SA_ERK_8_4_5 | ||
345 | getMethod (VERNER_8_5_6 _) = vERNER_8_5_6 | ||
346 | getMethod (VERNER_8_5_6') = vERNER_8_5_6 | ||
347 | getMethod (FEHLBERG_13_7_8 _) = fEHLBERG_13_7_8 | ||
348 | getMethod (FEHLBERG_13_7_8') = fEHLBERG_13_7_8 | ||
277 | 349 | ||
278 | getJacobian :: ODEMethod -> Maybe Jacobian | 350 | getJacobian :: ODEMethod -> Maybe Jacobian |
279 | getJacobian (SDIRK_2_1_2 j) = Just j | 351 | getJacobian (SDIRK_2_1_2 j) = Just j |
280 | getJacobian (SDIRK_2_1_2') = Nothing | 352 | getJacobian (BILLINGTON_3_3_2 j) = Just j |
281 | getJacobian (BILLINGTON_3_3_2 j) = Just j | 353 | getJacobian (TRBDF2_3_3_2 j) = Just j |
282 | getJacobian (BILLINGTON_3_3_2') = Nothing | 354 | getJacobian (KVAERNO_4_2_3 j) = Just j |
283 | getJacobian (TRBDF2_3_3_2 j) = Just j | 355 | getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j |
284 | getJacobian (TRBDF2_3_3_2') = Nothing | 356 | getJacobian (CASH_5_2_4 j) = Just j |
285 | getJacobian (KVAERNO_4_2_3 j) = Just j | 357 | getJacobian (CASH_5_3_4 j) = Just j |
286 | getJacobian (KVAERNO_4_2_3') = Nothing | 358 | getJacobian (SDIRK_5_3_4 j) = Just j |
287 | getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j | 359 | getJacobian (KVAERNO_5_3_4 j) = Just j |
288 | getJacobian (ARK324L2SA_DIRK_4_2_3') = Nothing | 360 | getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j |
289 | getJacobian (CASH_5_2_4 j) = Just j | 361 | getJacobian (KVAERNO_7_4_5 j) = Just j |
290 | getJacobian (CASH_5_2_4') = Nothing | 362 | getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j |
291 | getJacobian (CASH_5_3_4 j) = Just j | 363 | getJacobian (HEUN_EULER_2_1_2 j) = Just j |
292 | getJacobian (CASH_5_3_4') = Nothing | 364 | getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j |
293 | getJacobian (SDIRK_5_3_4 j) = Just j | 365 | getJacobian (ARK324L2SA_ERK_4_2_3 j) = Just j |
294 | getJacobian (SDIRK_5_3_4') = Nothing | 366 | getJacobian (ZONNEVELD_5_3_4 j) = Just j |
295 | getJacobian (KVAERNO_5_3_4 j) = Just j | 367 | getJacobian (ARK436L2SA_ERK_6_3_4 j) = Just j |
296 | getJacobian (KVAERNO_5_3_4') = Nothing | 368 | getJacobian (SAYFY_ABURUB_6_3_4 j) = Just j |
297 | getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j | 369 | getJacobian (CASH_KARP_6_4_5 j) = Just j |
298 | getJacobian (ARK436L2SA_DIRK_6_3_4') = Nothing | 370 | getJacobian (FEHLBERG_6_4_5 j) = Just j |
299 | getJacobian (KVAERNO_7_4_5 j) = Just j | 371 | getJacobian (DORMAND_PRINCE_7_4_5 j) = Just j |
300 | getJacobian (KVAERNO_7_4_5') = Nothing | 372 | getJacobian (ARK548L2SA_ERK_8_4_5 j) = Just j |
301 | getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j | 373 | getJacobian (VERNER_8_5_6 j) = Just j |
302 | getJacobian (ARK548L2SA_DIRK_8_4_5') = Nothing | 374 | getJacobian (FEHLBERG_13_7_8 j) = Just j |
303 | getJacobian (FEHLBERG_6_4_5 j) = Just j | 375 | getJacobian _ = Nothing |
304 | getJacobian (FEHLBERG_6_4_5' ) = Nothing | ||
305 | 376 | ||
306 | -- | A version of 'odeSolveVWith' with reasonable default step control. | 377 | -- | A version of 'odeSolveVWith' with reasonable default step control. |
307 | odeSolveV | 378 | odeSolveV |
@@ -515,6 +586,12 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO | |||
515 | if (check_flag(&flag, "ARKodeInit", 1)) return 1; | 586 | if (check_flag(&flag, "ARKodeInit", 1)) return 1; |
516 | } | 587 | } |
517 | 588 | ||
589 | /* FIXME: A hack for initial testing */ | ||
590 | flag = ARKodeSetMinStep(arkode_mem, 1.0e-12); | ||
591 | if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1; | ||
592 | flag = ARKodeSetMaxNumSteps(arkode_mem, 10000); | ||
593 | if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1; | ||
594 | |||
518 | /* Set routines */ | 595 | /* Set routines */ |
519 | flag = ARKodeSVtolerances(arkode_mem, $(double relTol), tv); | 596 | flag = ARKodeSVtolerances(arkode_mem, $(double relTol), tv); |
520 | if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1; | 597 | if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1; |