summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/sundials/hmatrix-sundials.cabal8
-rw-r--r--packages/sundials/src/Main.hs3
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs19
3 files changed, 14 insertions, 16 deletions
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
index dbe50e0..3304cb8 100644
--- a/packages/sundials/hmatrix-sundials.cabal
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -1,16 +1,16 @@
1-- Initial sundials.cabal generated by cabal init. For further 1-- Initial sundials.cabal generated by cabal init. For further
2-- documentation, see http://haskell.org/cabal/users-guide/ 2-- documentation, see http://haskell.org/cabal/users-guide/
3 3
4name: hmatrix-sundials 4name: hmatrix-sundials
5version: 0.1.0.0 5version: 0.1.0.0
6-- synopsis: 6-- synopsis:
7-- description: 7-- description:
8homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials 8homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials
9license: BSD3 9license: BSD3
10license-file: LICENSE 10license-file: LICENSE
11author: Dominic Steinitz 11author: Dominic Steinitz
12maintainer: dominic@steinitz.org 12maintainer: dominic@steinitz.org
13-- copyright: 13-- copyright:
14category: Math 14category: Math
15build-type: Simple 15build-type: Simple
16extra-source-files: ChangeLog.md 16extra-source-files: ChangeLog.md
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 895e610..c81d1a3 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -76,7 +76,7 @@ main = do
76 let res = btGet 76 let res = btGet
77 putStrLn $ show res 77 putStrLn $ show res
78 putStrLn $ butcherTableauTex res 78 putStrLn $ butcherTableauTex res
79 79
80 let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) 80 let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
81 putStrLn $ show res 81 putStrLn $ show res
82 renderRasterific "diagrams/brusselator.png" 82 renderRasterific "diagrams/brusselator.png"
@@ -88,4 +88,3 @@ main = do
88 renderRasterific "diagrams/stiffish.png" 88 renderRasterific "diagrams/stiffish.png"
89 (D.dims2D 500.0 500.0) 89 (D.dims2D 500.0 500.0)
90 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res)) 90 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res))
91
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 76ed61b..ff4ede8 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -58,10 +58,10 @@ C.include "<stdio.h>"
58C.include "<math.h>" 58C.include "<math.h>"
59C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. 59C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts.
60C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros 60C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros
61C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix 61C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
62C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver 62C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
63C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface 63C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface
64C.include "<sundials/sundials_types.h>" -- definition of type realtype 64C.include "<sundials/sundials_types.h>" -- definition of type realtype
65C.include "<sundials/sundials_math.h>" 65C.include "<sundials/sundials_math.h>"
66C.include "../../../helpers.h" 66C.include "../../../helpers.h"
67 67
@@ -126,7 +126,7 @@ odeSolve f y0 ts = case solveOde g (V.fromList y0) (V.fromList $ toList ts) of
126 nC = length y0 126 nC = length y0
127 g t x0 = V.fromList $ f t (V.toList x0) 127 g t x0 = V.fromList $ f t (V.toList x0)
128 128
129solveOde :: 129solveOde ::
130 (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 130 (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
131 -> V.Vector Double -- ^ Initial conditions 131 -> V.Vector Double -- ^ Initial conditions
132 -> V.Vector Double -- ^ Desired solution times 132 -> V.Vector Double -- ^ Desired solution times
@@ -134,7 +134,7 @@ solveOde ::
134solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of 134solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of
135 Left c -> Left $ fromIntegral c 135 Left c -> Left $ fromIntegral c
136 Right (v, d) -> Right (coerce v, d) 136 Right (v, d) -> Right (coerce v, d)
137 137
138solveOdeC :: 138solveOdeC ::
139 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 139 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
140 -> V.Vector CDouble -- ^ Initial conditions 140 -> V.Vector CDouble -- ^ Initial conditions
@@ -303,13 +303,13 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
303 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS); 303 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
304 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1); 304 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
305 ($vec-ptr:(long int *diagMut))[9] = ncfn; 305 ($vec-ptr:(long int *diagMut))[9] = ncfn;
306 306
307 /* Clean up and return */ 307 /* Clean up and return */
308 N_VDestroy(y); /* Free y vector */ 308 N_VDestroy(y); /* Free y vector */
309 ARKodeFree(&arkode_mem); /* Free integrator memory */ 309 ARKodeFree(&arkode_mem); /* Free integrator memory */
310 SUNLinSolFree(LS); /* Free linear solver */ 310 SUNLinSolFree(LS); /* Free linear solver */
311 SUNMatDestroy(A); /* Free A matrix */ 311 SUNMatDestroy(A); /* Free A matrix */
312 312
313 return flag; 313 return flag;
314 } |] 314 } |]
315 if res == 0 315 if res == 0
@@ -423,7 +423,7 @@ getButcherTable = unsafePerformIO $ do
423 /* Clean up and return */ 423 /* Clean up and return */
424 N_VDestroy(y); /* Free y vector */ 424 N_VDestroy(y); /* Free y vector */
425 ARKodeFree(&arkode_mem); /* Free integrator memory */ 425 ARKodeFree(&arkode_mem); /* Free integrator memory */
426 426
427 return flag; 427 return flag;
428 } |] 428 } |]
429 if res == 0 429 if res == 0
@@ -433,4 +433,3 @@ getButcherTable = unsafePerformIO $ do
433 return $ Right (x, y) 433 return $ Right (x, y)
434 else do 434 else do
435 return $ Left res 435 return $ Left res
436