summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs181
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
131import Data.Coerce (coerce) 137import Data.Coerce (coerce)
132import System.IO.Unsafe (unsafePerformIO) 138import System.IO.Unsafe (unsafePerformIO)
139import GHC.Generics
133 140
134import Numeric.LinearAlgebra.Devel (createVector) 141import 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
280constrName :: (HasConstructor (Rep a), Generic a)=> a -> String
281constrName = genericConstrName . from
282
283class HasConstructor (f :: * -> *) where
284 genericConstrName :: f x -> String
285
286instance HasConstructor f => HasConstructor (D1 c f) where
287 genericConstrName (M1 x) = genericConstrName x
288
289instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where
290 genericConstrName (L1 l) = genericConstrName l
291 genericConstrName (R1 r) = genericConstrName r
292
293instance Constructor c => HasConstructor (C1 c f) where
294 genericConstrName x = conName x
295
296instance Show ODEMethod where
297 show x = constrName x
298
299-- FIXME: We can probably do better here with generics
250getMethod :: ODEMethod -> Int 300getMethod :: ODEMethod -> Int
251getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2 301getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2
252getMethod (SDIRK_2_1_2') = sDIRK_2_1_2 302getMethod (SDIRK_2_1_2') = sDIRK_2_1_2
253getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2 303getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2
254getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2 304getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2
255getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2 305getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2
256getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2 306getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2
257getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3 307getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3
258getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3 308getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3
259getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3 309getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3
260getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3 310getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3
261getMethod (CASH_5_2_4 _) = cASH_5_2_4 311getMethod (CASH_5_2_4 _) = cASH_5_2_4
262getMethod (CASH_5_2_4') = cASH_5_2_4 312getMethod (CASH_5_2_4') = cASH_5_2_4
263getMethod (CASH_5_3_4 _) = cASH_5_3_4 313getMethod (CASH_5_3_4 _) = cASH_5_3_4
264getMethod (CASH_5_3_4') = cASH_5_3_4 314getMethod (CASH_5_3_4') = cASH_5_3_4
265getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4 315getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4
266getMethod (SDIRK_5_3_4') = sDIRK_5_3_4 316getMethod (SDIRK_5_3_4') = sDIRK_5_3_4
267getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4 317getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4
268getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4 318getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4
269getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4 319getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4
270getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4 320getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4
271getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5 321getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5
272getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5 322getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5
273getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5 323getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5
274getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5 324getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5
275getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5 325getMethod (HEUN_EULER_2_1_2 _) = hEUN_EULER_2_1_2
276getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5 326getMethod (HEUN_EULER_2_1_2') = hEUN_EULER_2_1_2
327getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3
328getMethod (BOGACKI_SHAMPINE_4_2_3') = bOGACKI_SHAMPINE_4_2_3
329getMethod (ARK324L2SA_ERK_4_2_3 _) = aRK324L2SA_ERK_4_2_3
330getMethod (ARK324L2SA_ERK_4_2_3') = aRK324L2SA_ERK_4_2_3
331getMethod (ZONNEVELD_5_3_4 _) = zONNEVELD_5_3_4
332getMethod (ZONNEVELD_5_3_4') = zONNEVELD_5_3_4
333getMethod (ARK436L2SA_ERK_6_3_4 _) = aRK436L2SA_ERK_6_3_4
334getMethod (ARK436L2SA_ERK_6_3_4') = aRK436L2SA_ERK_6_3_4
335getMethod (SAYFY_ABURUB_6_3_4 _) = sAYFY_ABURUB_6_3_4
336getMethod (SAYFY_ABURUB_6_3_4') = sAYFY_ABURUB_6_3_4
337getMethod (CASH_KARP_6_4_5 _) = cASH_KARP_6_4_5
338getMethod (CASH_KARP_6_4_5') = cASH_KARP_6_4_5
339getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5
340getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5
341getMethod (DORMAND_PRINCE_7_4_5 _) = dORMAND_PRINCE_7_4_5
342getMethod (DORMAND_PRINCE_7_4_5') = dORMAND_PRINCE_7_4_5
343getMethod (ARK548L2SA_ERK_8_4_5 _) = aRK548L2SA_ERK_8_4_5
344getMethod (ARK548L2SA_ERK_8_4_5') = aRK548L2SA_ERK_8_4_5
345getMethod (VERNER_8_5_6 _) = vERNER_8_5_6
346getMethod (VERNER_8_5_6') = vERNER_8_5_6
347getMethod (FEHLBERG_13_7_8 _) = fEHLBERG_13_7_8
348getMethod (FEHLBERG_13_7_8') = fEHLBERG_13_7_8
277 349
278getJacobian :: ODEMethod -> Maybe Jacobian 350getJacobian :: ODEMethod -> Maybe Jacobian
279getJacobian (SDIRK_2_1_2 j) = Just j 351getJacobian (SDIRK_2_1_2 j) = Just j
280getJacobian (SDIRK_2_1_2') = Nothing 352getJacobian (BILLINGTON_3_3_2 j) = Just j
281getJacobian (BILLINGTON_3_3_2 j) = Just j 353getJacobian (TRBDF2_3_3_2 j) = Just j
282getJacobian (BILLINGTON_3_3_2') = Nothing 354getJacobian (KVAERNO_4_2_3 j) = Just j
283getJacobian (TRBDF2_3_3_2 j) = Just j 355getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j
284getJacobian (TRBDF2_3_3_2') = Nothing 356getJacobian (CASH_5_2_4 j) = Just j
285getJacobian (KVAERNO_4_2_3 j) = Just j 357getJacobian (CASH_5_3_4 j) = Just j
286getJacobian (KVAERNO_4_2_3') = Nothing 358getJacobian (SDIRK_5_3_4 j) = Just j
287getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j 359getJacobian (KVAERNO_5_3_4 j) = Just j
288getJacobian (ARK324L2SA_DIRK_4_2_3') = Nothing 360getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j
289getJacobian (CASH_5_2_4 j) = Just j 361getJacobian (KVAERNO_7_4_5 j) = Just j
290getJacobian (CASH_5_2_4') = Nothing 362getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j
291getJacobian (CASH_5_3_4 j) = Just j 363getJacobian (HEUN_EULER_2_1_2 j) = Just j
292getJacobian (CASH_5_3_4') = Nothing 364getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j
293getJacobian (SDIRK_5_3_4 j) = Just j 365getJacobian (ARK324L2SA_ERK_4_2_3 j) = Just j
294getJacobian (SDIRK_5_3_4') = Nothing 366getJacobian (ZONNEVELD_5_3_4 j) = Just j
295getJacobian (KVAERNO_5_3_4 j) = Just j 367getJacobian (ARK436L2SA_ERK_6_3_4 j) = Just j
296getJacobian (KVAERNO_5_3_4') = Nothing 368getJacobian (SAYFY_ABURUB_6_3_4 j) = Just j
297getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j 369getJacobian (CASH_KARP_6_4_5 j) = Just j
298getJacobian (ARK436L2SA_DIRK_6_3_4') = Nothing 370getJacobian (FEHLBERG_6_4_5 j) = Just j
299getJacobian (KVAERNO_7_4_5 j) = Just j 371getJacobian (DORMAND_PRINCE_7_4_5 j) = Just j
300getJacobian (KVAERNO_7_4_5') = Nothing 372getJacobian (ARK548L2SA_ERK_8_4_5 j) = Just j
301getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j 373getJacobian (VERNER_8_5_6 j) = Just j
302getJacobian (ARK548L2SA_DIRK_8_4_5') = Nothing 374getJacobian (FEHLBERG_13_7_8 j) = Just j
303getJacobian (FEHLBERG_6_4_5 j) = Just j 375getJacobian _ = Nothing
304getJacobian (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.
307odeSolveV 378odeSolveV
@@ -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;
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 12:08:02 +0000