Support RKF45

3 files changed, 44 insertions, 9 deletions

diff --git a/packages/sundials/src/Arkode.hsc b/packages/sundials/src/Arkode.hsc
index 2904c36..4ef95e2 100644
--- a/packages/sundials/src/Arkode.hsc
+++ b/packages/sundials/src/Arkode.hsc
@@ -46,6 +46,8 @@ arkSMax :: Int
 arkSMax = #const ARK_S_MAX
 
 -- FIXME: We could just use inline-c instead
+
+-- /* Butcher table accessors -- implicit */
 sDIRK_2_1_2 :: Int
 sDIRK_2_1_2 = #const SDIRK_2_1_2
 -- #define BILLINGTON_3_3_2        13
@@ -67,3 +69,26 @@ sDIRK_5_3_4 = #const SDIRK_5_3_4
 -- #define DEFAULT_DIRK_3          ARK324L2SA_DIRK_4_2_3
 -- #define DEFAULT_DIRK_4          SDIRK_5_3_4
 -- #define DEFAULT_DIRK_5          ARK548L2SA_DIRK_8_4_5
+
+-- /* Butcher table accessors -- explicit */
+-- #define HEUN_EULER_2_1_2         0
+-- #define BOGACKI_SHAMPINE_4_2_3   1
+-- #define ARK324L2SA_ERK_4_2_3     2
+-- #define ZONNEVELD_5_3_4          3
+-- #define ARK436L2SA_ERK_6_3_4     4
+-- #define SAYFY_ABURUB_6_3_4       5
+-- #define CASH_KARP_6_4_5          6
+fEHLBERG_6_4_5 :: Int
+fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5
+-- #define FEHLBERG_6_4_5           7
+-- #define DORMAND_PRINCE_7_4_5     8
+-- #define ARK548L2SA_ERK_8_4_5     9
+-- #define VERNER_8_5_6            10
+-- #define FEHLBERG_13_7_8         11
+
+-- #define DEFAULT_ERK_2           HEUN_EULER_2_1_2
+-- #define DEFAULT_ERK_3           BOGACKI_SHAMPINE_4_2_3
+-- #define DEFAULT_ERK_4           ZONNEVELD_5_3_4
+-- #define DEFAULT_ERK_5           CASH_KARP_6_4_5
+-- #define DEFAULT_ERK_6           VERNER_8_5_6
+-- #define DEFAULT_ERK_8           FEHLBERG_13_7_8
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index e96522f..d8f3f3d 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -128,6 +128,10 @@ main = do
   putStrLn $ show resB
   putStrLn $ butcherTableauTex resB
 
+  -- let resC = butcherTable (FEHLBERG_6_4_5 undefined)
+  -- putStrLn $ show resC
+  -- putStrLn $ butcherTableauTex resC
+
   let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
   renderRasterific "diagrams/brusselator.png"
                    (D.dims2D 500.0 500.0)
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 444138d..380b1d6 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -137,7 +137,7 @@ import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
                                                 size, subVector)
 
 import qualified Types as T
-import           Arkode (sDIRK_2_1_2, kVAERNO_4_2_3, sDIRK_5_3_4)
+import           Arkode (sDIRK_2_1_2, kVAERNO_4_2_3, sDIRK_5_3_4, fEHLBERG_6_4_5)
 import qualified Arkode as B
 
 import Debug.Trace
@@ -225,18 +225,24 @@ data ODEMethod = SDIRK_2_1_2   Jacobian
                | KVAERNO_4_2_3 Jacobian
                | SDIRK_5_3_4   Jacobian
                | SDIRK_5_3_4'
+               | FEHLBERG_6_4_5 Jacobian
+               | FEHLBERG_6_4_5'
 
 getMethod :: ODEMethod -> Int
-getMethod (SDIRK_2_1_2   _) = sDIRK_2_1_2
-getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3
-getMethod (SDIRK_5_3_4   _) = sDIRK_5_3_4
-getMethod (SDIRK_5_3_4'   ) = sDIRK_5_3_4
+getMethod (SDIRK_2_1_2    _) = sDIRK_2_1_2
+getMethod (KVAERNO_4_2_3  _) = kVAERNO_4_2_3
+getMethod (SDIRK_5_3_4    _) = sDIRK_5_3_4
+getMethod (SDIRK_5_3_4'    ) = sDIRK_5_3_4
+getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5
+getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5
 
 getJacobian :: ODEMethod -> Maybe Jacobian
-getJacobian (SDIRK_2_1_2   j) = Just j
-getJacobian (KVAERNO_4_2_3 j) = Just j
-getJacobian (SDIRK_5_3_4   j) = Just j
-getJacobian (SDIRK_5_3_4'   ) = Nothing
+getJacobian (SDIRK_2_1_2    j) = Just j
+getJacobian (KVAERNO_4_2_3  j) = Just j
+getJacobian (SDIRK_5_3_4    j) = Just j
+getJacobian (SDIRK_5_3_4'    ) = Nothing
+getJacobian (FEHLBERG_6_4_5 j) = Just j
+getJacobian (FEHLBERG_6_4_5' ) = Nothing
 
 -- | A version of 'odeSolveVWith' with reasonable default step control.
 odeSolveV