summaryrefslogtreecommitdiff
path: root/packages/sundials/src
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials/src')
-rw-r--r--packages/sundials/src/Main.hs34
-rw-r--r--packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs (renamed from packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs)31
2 files changed, 31 insertions, 34 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 5e51372..d1f6755 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -1,6 +1,5 @@
1{-# OPTIONS_GHC -Wall #-} 1module Main where
2 2
3import qualified Data.Vector.Storable as V
4import Numeric.Sundials.Arkode.ODE 3import Numeric.Sundials.Arkode.ODE
5import Numeric.LinearAlgebra 4import Numeric.LinearAlgebra
6 5
@@ -9,11 +8,10 @@ import qualified Diagrams.Prelude as D
9import Diagrams.Backend.Rasterific 8import Diagrams.Backend.Rasterific
10 9
11import Control.Lens 10import Control.Lens
12import Data.List (zip4)
13 11
14 12brusselator :: a -> [Double] -> [Double]
15brusselator _t x = [ a - (w + 1) * u + v * u^2 13brusselator _t x = [ a - (w + 1) * u + v * u^(2::Int)
16 , w * u - v * u^2 14 , w * u - v * u^(2::Int)
17 , (b - w) / eps - w * u 15 , (b - w) / eps - w * u
18 ] 16 ]
19 where 17 where
@@ -24,6 +22,7 @@ brusselator _t x = [ a - (w + 1) * u + v * u^2
24 v = x !! 1 22 v = x !! 1
25 w = x !! 2 23 w = x !! 2
26 24
25stiffish :: Double -> [Double] -> [Double]
27stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] 26stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
28 where 27 where
29 lamda = -100.0 28 lamda = -100.0
@@ -45,15 +44,14 @@ kSaxis xs = P.r2Axis &~ do
45 44
46main :: IO () 45main :: IO ()
47main = do 46main = do
48 let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) 47 do let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
49 putStrLn $ show res 48 putStrLn $ show res
50 renderRasterific "diagrams/brusselator.png" 49 renderRasterific "diagrams/brusselator.png"
51 (D.dims2D 500.0 500.0) 50 (D.dims2D 500.0 500.0)
52 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res)) 51 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res))
53 52
54 let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) 53 do let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
55 putStrLn $ show res 54 putStrLn $ show res
56 renderRasterific "diagrams/stiffish.png" 55 renderRasterific "diagrams/stiffish.png"
57 (D.dims2D 500.0 500.0) 56 (D.dims2D 500.0 500.0)
58 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res)) 57 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res))
59
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs
index f432951..6d9a1b2 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs
@@ -1,14 +1,14 @@
1{-# OPTIONS_GHC -Wall #-}
2
3{-# LANGUAGE QuasiQuotes #-}
4{-# LANGUAGE TemplateHaskell #-}
5{-# LANGUAGE MultiWayIf #-} 1{-# LANGUAGE MultiWayIf #-}
6{-# LANGUAGE OverloadedStrings #-} 2{-# LANGUAGE OverloadedStrings #-}
3{-# LANGUAGE QuasiQuotes #-}
7{-# LANGUAGE ScopedTypeVariables #-} 4{-# LANGUAGE ScopedTypeVariables #-}
5{-# LANGUAGE TemplateHaskell #-}
8 6
9module Numeric.Sundials.Arkode.ODE ( solveOde 7module Numeric.Sundials.Arkode.ODE
10 , odeSolve 8 ( SundialsDiagnostics(..)
11 ) where 9 , solveOde
10 , odeSolve
11 ) where
12 12
13import qualified Language.C.Inline as C 13import qualified Language.C.Inline as C
14import qualified Language.C.Inline.Unsafe as CU 14import qualified Language.C.Inline.Unsafe as CU
@@ -40,10 +40,10 @@ C.include "<stdio.h>"
40C.include "<math.h>" 40C.include "<math.h>"
41C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. 41C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts.
42C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros 42C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros
43C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix 43C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
44C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver 44C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
45C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface 45C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface
46C.include "<sundials/sundials_types.h>" -- definition of type realtype 46C.include "<sundials/sundials_types.h>" -- definition of type realtype
47C.include "<sundials/sundials_math.h>" 47C.include "<sundials/sundials_math.h>"
48C.include "../../../helpers.h" 48C.include "../../../helpers.h"
49 49
@@ -108,7 +108,7 @@ odeSolve f y0 ts = case solveOde g (V.fromList y0) (V.fromList $ toList ts) of
108 nC = length y0 108 nC = length y0
109 g t x0 = V.fromList $ f t (V.toList x0) 109 g t x0 = V.fromList $ f t (V.toList x0)
110 110
111solveOde :: 111solveOde ::
112 (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 112 (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
113 -> V.Vector Double -- ^ Initial conditions 113 -> V.Vector Double -- ^ Initial conditions
114 -> V.Vector Double -- ^ Desired solution times 114 -> V.Vector Double -- ^ Desired solution times
@@ -116,7 +116,7 @@ solveOde ::
116solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of 116solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of
117 Left c -> Left $ fromIntegral c 117 Left c -> Left $ fromIntegral c
118 Right (v, d) -> Right (coerce v, d) 118 Right (v, d) -> Right (coerce v, d)
119 119
120solveOdeC :: 120solveOdeC ::
121 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 121 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
122 -> V.Vector CDouble -- ^ Initial conditions 122 -> V.Vector CDouble -- ^ Initial conditions
@@ -285,13 +285,13 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
285 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS); 285 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
286 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1); 286 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
287 ($vec-ptr:(long int *diagMut))[9] = ncfn; 287 ($vec-ptr:(long int *diagMut))[9] = ncfn;
288 288
289 /* Clean up and return */ 289 /* Clean up and return */
290 N_VDestroy(y); /* Free y vector */ 290 N_VDestroy(y); /* Free y vector */
291 ARKodeFree(&arkode_mem); /* Free integrator memory */ 291 ARKodeFree(&arkode_mem); /* Free integrator memory */
292 SUNLinSolFree(LS); /* Free linear solver */ 292 SUNLinSolFree(LS); /* Free linear solver */
293 SUNMatDestroy(A); /* Free A matrix */ 293 SUNMatDestroy(A); /* Free A matrix */
294 294
295 return flag; 295 return flag;
296 } |] 296 } |]
297 if res == 0 297 if res == 0
@@ -311,4 +311,3 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
311 return $ Right (m, d) 311 return $ Right (m, d)
312 else do 312 else do
313 return $ Left res 313 return $ Left res
314