diff options
-rw-r--r-- | packages/sundials/hmatrix-sundials.cabal | 8 | ||||
-rw-r--r-- | packages/sundials/src/Main.hs | 3 | ||||
-rw-r--r-- | packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs | 19 |
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 | ||
4 | name: hmatrix-sundials | 4 | name: hmatrix-sundials |
5 | version: 0.1.0.0 | 5 | version: 0.1.0.0 |
6 | -- synopsis: | 6 | -- synopsis: |
7 | -- description: | 7 | -- description: |
8 | homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials | 8 | homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials |
9 | license: BSD3 | 9 | license: BSD3 |
10 | license-file: LICENSE | 10 | license-file: LICENSE |
11 | author: Dominic Steinitz | 11 | author: Dominic Steinitz |
12 | maintainer: dominic@steinitz.org | 12 | maintainer: dominic@steinitz.org |
13 | -- copyright: | 13 | -- copyright: |
14 | category: Math | 14 | category: Math |
15 | build-type: Simple | 15 | build-type: Simple |
16 | extra-source-files: ChangeLog.md | 16 | extra-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>" | |||
58 | C.include "<math.h>" | 58 | C.include "<math.h>" |
59 | C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. | 59 | C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. |
60 | C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros | 60 | C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros |
61 | C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix | 61 | C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix |
62 | C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver | 62 | C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver |
63 | C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface | 63 | C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface |
64 | C.include "<sundials/sundials_types.h>" -- definition of type realtype | 64 | C.include "<sundials/sundials_types.h>" -- definition of type realtype |
65 | C.include "<sundials/sundials_math.h>" | 65 | C.include "<sundials/sundials_math.h>" |
66 | C.include "../../../helpers.h" | 66 | C.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 | ||
129 | solveOde :: | 129 | solveOde :: |
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 :: | |||
134 | solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of | 134 | solveOde 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 | ||
138 | solveOdeC :: | 138 | solveOdeC :: |
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 | |||