summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/CVode/ODE.hs')
-rw-r--r--packages/sundials/src/Numeric/Sundials/CVode/ODE.hs23
1 files changed, 14 insertions, 9 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
index 159fbe2..a6f185e 100644
--- a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
@@ -68,7 +68,6 @@ module Numeric.Sundials.CVode.ODE ( odeSolve
68 , odeSolveVWith' 68 , odeSolveVWith'
69 , ODEMethod(..) 69 , ODEMethod(..)
70 , StepControl(..) 70 , StepControl(..)
71 , Jacobian
72 ) where 71 ) where
73 72
74import qualified Language.C.Inline as C 73import qualified Language.C.Inline as C
@@ -127,7 +126,7 @@ getJacobian _ = Nothing
127-- | A version of 'odeSolveVWith' with reasonable default step control. 126-- | A version of 'odeSolveVWith' with reasonable default step control.
128odeSolveV 127odeSolveV
129 :: ODEMethod 128 :: ODEMethod
130 -> Maybe Double -- ^ initial step size - by default, ARKode 129 -> Maybe Double -- ^ initial step size - by default, CVode
131 -- estimates the initial step size to be the 130 -- estimates the initial step size to be the
132 -- solution \(h\) of the equation 131 -- solution \(h\) of the equation
133 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where 132 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
@@ -161,7 +160,7 @@ odeSolve f y0 ts =
161odeSolveVWith :: 160odeSolveVWith ::
162 ODEMethod 161 ODEMethod
163 -> StepControl 162 -> StepControl
164 -> Maybe Double -- ^ initial step size - by default, ARKode 163 -> Maybe Double -- ^ initial step size - by default, CVode
165 -- estimates the initial step size to be the 164 -- estimates the initial step size to be the
166 -- solution \(h\) of the equation 165 -- solution \(h\) of the equation
167 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where 166 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
@@ -181,13 +180,14 @@ odeSolveVWith method control initStepSize f y0 tt =
181 , relTol = error "relTol" 180 , relTol = error "relTol"
182 , absTols = error "absTol" 181 , absTols = error "absTol"
183 , initStep = error "initStep" 182 , initStep = error "initStep"
183 , maxFail = 10
184 } 184 }
185 185
186odeSolveVWith' :: 186odeSolveVWith' ::
187 ODEOpts 187 ODEOpts
188 -> ODEMethod 188 -> ODEMethod
189 -> StepControl 189 -> StepControl
190 -> Maybe Double -- ^ initial step size - by default, ARKode 190 -> Maybe Double -- ^ initial step size - by default, CVode
191 -- estimates the initial step size to be the 191 -- estimates the initial step size to be the
192 -- solution \(h\) of the equation 192 -- solution \(h\) of the equation
193 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where 193 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
@@ -198,13 +198,13 @@ odeSolveVWith' ::
198 -> V.Vector Double -- ^ Desired solution times 198 -> V.Vector Double -- ^ Desired solution times
199 -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution 199 -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
200odeSolveVWith' opts method control initStepSize f y0 tt = 200odeSolveVWith' opts method control initStepSize f y0 tt =
201 case solveOdeC (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) 201 case solveOdeC (fromIntegral $ maxFail opts)
202 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
202 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) 203 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
203 (coerce f) (coerce y0) (coerce tt) of 204 (coerce f) (coerce y0) (coerce tt) of
204 Left c -> Left $ fromIntegral c 205 Left c -> Left $ fromIntegral c
205 Right (v, d) -> Right (reshape nC (coerce v), d) 206 Right (v, d) -> Right (reshape l (coerce v), d)
206 where 207 where
207 nC = V.length y0
208 l = size y0 208 l = size y0
209 scise (X aTol rTol) = coerce (V.replicate l aTol, rTol) 209 scise (X aTol rTol) = coerce (V.replicate l aTol, rTol)
210 scise (X' aTol rTol) = coerce (V.replicate l aTol, rTol) 210 scise (X' aTol rTol) = coerce (V.replicate l aTol, rTol)
@@ -221,6 +221,7 @@ odeSolveVWith' opts method control initStepSize f y0 tt =
221 vs = V.fromList $ map coerce $ concat $ toLists m 221 vs = V.fromList $ map coerce $ concat $ toLists m
222 222
223solveOdeC :: 223solveOdeC ::
224 CInt ->
224 CLong -> 225 CLong ->
225 CDouble -> 226 CDouble ->
226 CInt -> 227 CInt ->
@@ -231,7 +232,8 @@ solveOdeC ::
231 -> V.Vector CDouble -- ^ Initial conditions 232 -> V.Vector CDouble -- ^ Initial conditions
232 -> V.Vector CDouble -- ^ Desired solution times 233 -> V.Vector CDouble -- ^ Desired solution times
233 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution 234 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
234solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts = 235solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
236 jacH (aTols, rTol) fun f0 ts =
235 unsafePerformIO $ do 237 unsafePerformIO $ do
236 238
237 let isInitStepSize :: CInt 239 let isInitStepSize :: CInt
@@ -243,6 +245,7 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
243 -- used :( 245 -- used :(
244 Nothing -> 0.0 246 Nothing -> 0.0
245 Just x -> x 247 Just x -> x
248
246 let dim = V.length f0 249 let dim = V.length f0
247 nEq :: CLong 250 nEq :: CLong
248 nEq = fromIntegral dim 251 nEq = fromIntegral dim
@@ -271,7 +274,7 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
271 IO CInt 274 IO CInt
272 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do 275 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
273 case jacH of 276 case jacH of
274 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined" 277 Nothing -> error "Numeric.Sundials.CVode.ODE: Jacobian not defined"
275 Just jacI -> do j <- jacI t <$> getDataFromContents dim y 278 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
276 poke jacS j 279 poke jacS j
277 -- FIXME: I don't understand what this comment means 280 -- FIXME: I don't understand what this comment means
@@ -326,6 +329,8 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
326 if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1; 329 if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1;
327 flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_)); 330 flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_));
328 if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1; 331 if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1;
332 flag = CVodeSetMaxErrTestFails(cvode_mem, $(int maxErrTestFails));
333 if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return 1;
329 334
330 /* Call CVodeSVtolerances to specify the scalar relative tolerance 335 /* Call CVodeSVtolerances to specify the scalar relative tolerance
331 * and vector absolute tolerances */ 336 * and vector absolute tolerances */