diff options
author | Dominic Steinitz <dominic@steinitz.org> | 2018-05-02 14:42:43 +0100 |
---|---|---|
committer | Dominic Steinitz <dominic@steinitz.org> | 2018-05-02 14:42:43 +0100 |
commit | 4ba859636396d211637b5507f19722b6953656a5 (patch) | |
tree | 9493c4851e6141a400e6345efe59a07197709f63 /packages/sundials/src/Numeric/Sundials/CVode/ODE.hs | |
parent | 149dedfc6ec8dea039a4df7ad1d31880820c52eb (diff) |
Add more options
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/CVode/ODE.hs')
-rw-r--r-- | packages/sundials/src/Numeric/Sundials/CVode/ODE.hs | 23 |
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 | ||
74 | import qualified Language.C.Inline as C | 73 | import 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. |
128 | odeSolveV | 127 | odeSolveV |
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 = | |||
161 | odeSolveVWith :: | 160 | odeSolveVWith :: |
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 | ||
186 | odeSolveVWith' :: | 186 | odeSolveVWith' :: |
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 |
200 | odeSolveVWith' opts method control initStepSize f y0 tt = | 200 | odeSolveVWith' 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 | ||
223 | solveOdeC :: | 223 | solveOdeC :: |
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 |
234 | solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts = | 235 | solveOdeC 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 */ |