summaryrefslogtreecommitdiff
path: root/packages/sundials
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials')
-rw-r--r--packages/sundials/src/Main.hs14
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs45
2 files changed, 36 insertions, 23 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index b950036..9978aa5 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -127,15 +127,15 @@ main = do
127 -- putStrLn $ show res 127 -- putStrLn $ show res
128 -- putStrLn $ butcherTableauTex res 128 -- putStrLn $ butcherTableauTex res
129 129
130 let res1 = odeSolve' (SDIRK_5_3_4 brussJac) brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 100.0]) 130 let res1 = odeSolve' (SDIRK_5_3_4 brussJac) brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
131 renderRasterific "diagrams/brusselator.png" 131 renderRasterific "diagrams/brusselator.png"
132 (D.dims2D 500.0 500.0) 132 (D.dims2D 500.0 500.0)
133 (renderAxis $ lSaxis $ [0.0, 0.1 .. 100.0]:(toLists $ tr res1)) 133 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
134 134
135 let res1a = odeSolve' (SDIRK_5_3_4') brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 100.0]) 135 let res1a = odeSolve' (SDIRK_5_3_4') brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
136 renderRasterific "diagrams/brusselatorA.png" 136 renderRasterific "diagrams/brusselatorA.png"
137 (D.dims2D 500.0 500.0) 137 (D.dims2D 500.0 500.0)
138 (renderAxis $ lSaxis $ [0.0, 0.1 .. 100.0]:(toLists $ tr res1a)) 138 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))
139 139
140 let res2 = odeSolve' (SDIRK_5_3_4 stiffJac) stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) 140 let res2 = odeSolve' (SDIRK_5_3_4 stiffJac) stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
141 putStrLn $ show res2 141 putStrLn $ show res2
@@ -143,17 +143,17 @@ main = do
143 (D.dims2D 500.0 500.0) 143 (D.dims2D 500.0 500.0)
144 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) 144 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2))
145 145
146 let res3 = odeSolve' (SDIRK_5_3_4 lorenzJac) lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 1000.0]) 146 let res3 = odeSolve' (SDIRK_5_3_4 lorenzJac) lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
147 putStrLn $ show $ last ((toLists $ tr res3)!!0) 147 putStrLn $ show $ last ((toLists $ tr res3)!!0)
148 148
149 let res3 = odeSolve' (SDIRK_5_3_4') lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 1000.0]) 149 let res3 = odeSolve' (SDIRK_5_3_4') lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
150 putStrLn $ show $ last ((toLists $ tr res3)!!0) 150 putStrLn $ show $ last ((toLists $ tr res3)!!0)
151 151
152 renderRasterific "diagrams/lorenz.png" 152 renderRasterific "diagrams/lorenz.png"
153 (D.dims2D 500.0 500.0) 153 (D.dims2D 500.0 500.0)
154 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1)) 154 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))
155 155
156 let res3a = odeSolve' (SDIRK_5_3_4') lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 30.0]) 156 let res3a = odeSolve' (SDIRK_5_3_4') lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
157 renderRasterific "diagrams/lorenzA.png" 157 renderRasterific "diagrams/lorenzA.png"
158 (D.dims2D 500.0 500.0) 158 (D.dims2D 500.0 500.0)
159 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3a)!!1)) 159 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3a)!!1))
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 4270a13..a7e302d 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -76,6 +76,7 @@ module Numeric.Sundials.ARKode.ODE ( odeSolve
76 , btGet 76 , btGet
77 , ODEMethod(..) 77 , ODEMethod(..)
78 , StepControl(..) 78 , StepControl(..)
79 , SundialsDiagnostics(..)
79 ) where 80 ) where
80 81
81import qualified Language.C.Inline as C 82import qualified Language.C.Inline as C
@@ -98,7 +99,8 @@ import System.IO.Unsafe (unsafePerformIO)
98import Numeric.LinearAlgebra.Devel (createVector) 99import Numeric.LinearAlgebra.Devel (createVector)
99 100
100import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><), 101import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
101 subMatrix, rows, cols, toLists) 102 subMatrix, rows, cols, toLists,
103 size)
102 104
103import qualified Types as T 105import qualified Types as T
104import Arkode (sDIRK_2_1_2, kVAERNO_4_2_3, sDIRK_5_3_4) 106import Arkode (sDIRK_2_1_2, kVAERNO_4_2_3, sDIRK_5_3_4)
@@ -170,16 +172,16 @@ vectorToC vec len ptr = do
170 V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec 172 V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
171 173
172data SundialsDiagnostics = SundialsDiagnostics { 174data SundialsDiagnostics = SundialsDiagnostics {
173 _aRKodeGetNumSteps :: Int 175 aRKodeGetNumSteps :: Int
174 , _aRKodeGetNumStepAttempts :: Int 176 , aRKodeGetNumStepAttempts :: Int
175 , _aRKodeGetNumRhsEvals_fe :: Int 177 , aRKodeGetNumRhsEvals_fe :: Int
176 , _aRKodeGetNumRhsEvals_fi :: Int 178 , aRKodeGetNumRhsEvals_fi :: Int
177 , _aRKodeGetNumLinSolvSetups :: Int 179 , aRKodeGetNumLinSolvSetups :: Int
178 , _aRKodeGetNumErrTestFails :: Int 180 , aRKodeGetNumErrTestFails :: Int
179 , _aRKodeGetNumNonlinSolvIters :: Int 181 , aRKodeGetNumNonlinSolvIters :: Int
180 , _aRKodeGetNumNonlinSolvConvFails :: Int 182 , aRKodeGetNumNonlinSolvConvFails :: Int
181 , _aRKDlsGetNumJacEvals :: Int 183 , aRKDlsGetNumJacEvals :: Int
182 , _aRKDlsGetNumRhsEvals :: Int 184 , aRKDlsGetNumRhsEvals :: Int
183 } deriving Show 185 } deriving Show
184 186
185type Jacobian = Double -> Vector Double -> Matrix Double 187type Jacobian = Double -> Vector Double -> Matrix Double
@@ -208,11 +210,20 @@ odeSolveV
208 -> Double -- ^ initial step size 210 -> Double -- ^ initial step size
209 -> Double -- ^ absolute tolerance for the state vector 211 -> Double -- ^ absolute tolerance for the state vector
210 -> Double -- ^ relative tolerance for the state vector 212 -> Double -- ^ relative tolerance for the state vector
211 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) 213 -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
212 -> Vector Double -- ^ initial conditions 214 -> Vector Double -- ^ initial conditions
213 -> Vector Double -- ^ desired solution times 215 -> Vector Double -- ^ desired solution times
214 -> Matrix Double -- ^ solution 216 -> Matrix Double -- ^ solution
215odeSolveV _meth _hi _epsAbs _epsRel = undefined 217odeSolveV meth _hi epsAbs epsRel f y0 ts =
218 case odeSolveVWith meth (XX' epsAbs epsRel 1 1) epsAbs epsAbs g y0 ts of
219 Left c -> error $ show c -- FIXME
220 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v)
221 where
222 us = toList ts
223 nR = length us
224 nC = size y0
225 -- FIXME: via a list - really?
226 g t x0 = V.fromList $ toList $ f t x0
216 227
217-- | A version of 'odeSolveV' with reasonable default parameters and 228-- | A version of 'odeSolveV' with reasonable default parameters and
218-- system of equations defined using lists. FIXME: we should say 229-- system of equations defined using lists. FIXME: we should say
@@ -223,7 +234,8 @@ odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y
223 -> Vector Double -- ^ desired solution times 234 -> Vector Double -- ^ desired solution times
224 -> Matrix Double -- ^ solution 235 -> Matrix Double -- ^ solution
225odeSolve f y0 ts = 236odeSolve f y0 ts =
226 case odeSolveVWith SDIRK_5_3_4' 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of 237 -- FIXME: These tolerances are different from the ones in GSL
238 case odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1) 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of
227 Left c -> error $ show c -- FIXME 239 Left c -> error $ show c -- FIXME
228 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v) 240 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v)
229 where 241 where
@@ -238,7 +250,7 @@ odeSolve' :: ODEMethod
238 -> Vector Double -- ^ desired solution times 250 -> Vector Double -- ^ desired solution times
239 -> Matrix Double -- ^ solution 251 -> Matrix Double -- ^ solution
240odeSolve' method f y0 ts = 252odeSolve' method f y0 ts =
241 case odeSolveVWith method 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of 253 case odeSolveVWith method (XX' 1.0e-6 1.0e-10 1 1) 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of
242 Left c -> error $ show c -- FIXME 254 Left c -> error $ show c -- FIXME
243 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v) 255 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v)
244 where 256 where
@@ -249,13 +261,14 @@ odeSolve' method f y0 ts =
249 261
250odeSolveVWith :: 262odeSolveVWith ::
251 ODEMethod 263 ODEMethod
264 -> StepControl
252 -> Double 265 -> Double
253 -> Double 266 -> Double
254 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 267 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
255 -> V.Vector Double -- ^ Initial conditions 268 -> V.Vector Double -- ^ Initial conditions
256 -> V.Vector Double -- ^ Desired solution times 269 -> V.Vector Double -- ^ Desired solution times
257 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution 270 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution
258odeSolveVWith method relTol absTol f y0 tt = 271odeSolveVWith method _control relTol absTol f y0 tt =
259 case solveOdeC (fromIntegral $ getMethod method) jacH (CDouble relTol) (CDouble absTol) 272 case solveOdeC (fromIntegral $ getMethod method) jacH (CDouble relTol) (CDouble absTol)
260 (coerce f) (coerce y0) (coerce tt) of 273 (coerce f) (coerce y0) (coerce tt) of
261 Left c -> Left $ fromIntegral c 274 Left c -> Left $ fromIntegral c