diff options
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/CVode/ODE.hs')
-rw-r--r-- | packages/sundials/src/Numeric/Sundials/CVode/ODE.hs | 45 |
1 files changed, 20 insertions, 25 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs index a6f185e..ad7cf51 100644 --- a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs +++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs | |||
@@ -172,7 +172,7 @@ odeSolveVWith :: | |||
172 | -> Matrix Double -- ^ Error code or solution | 172 | -> Matrix Double -- ^ Error code or solution |
173 | odeSolveVWith method control initStepSize f y0 tt = | 173 | odeSolveVWith method control initStepSize f y0 tt = |
174 | case odeSolveVWith' opts method control initStepSize f y0 tt of | 174 | case odeSolveVWith' opts method control initStepSize f y0 tt of |
175 | Left c -> error $ show c -- FIXME | 175 | Left (c, _v) -> error $ show c -- FIXME |
176 | Right (v, _d) -> v | 176 | Right (v, _d) -> v |
177 | where | 177 | where |
178 | opts = ODEOpts { maxNumSteps = 10000 | 178 | opts = ODEOpts { maxNumSteps = 10000 |
@@ -196,13 +196,13 @@ odeSolveVWith' :: | |||
196 | -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) | 196 | -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) |
197 | -> V.Vector Double -- ^ Initial conditions | 197 | -> V.Vector Double -- ^ Initial conditions |
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 (Matrix Double, 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 $ maxFail opts) | 201 | case solveOdeC (fromIntegral $ maxFail opts) |
202 | (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) | 202 | (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) |
203 | (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) | 203 | (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) |
204 | (coerce f) (coerce y0) (coerce tt) of | 204 | (coerce f) (coerce y0) (coerce tt) of |
205 | Left c -> Left $ fromIntegral c | 205 | Left (v, c) -> Left (reshape l (coerce v), fromIntegral c) |
206 | Right (v, d) -> Right (reshape l (coerce v), d) | 206 | Right (v, d) -> Right (reshape l (coerce v), d) |
207 | where | 207 | where |
208 | l = size y0 | 208 | l = size y0 |
@@ -231,7 +231,8 @@ solveOdeC :: | |||
231 | (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) | 231 | (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) |
232 | -> V.Vector CDouble -- ^ Initial conditions | 232 | -> V.Vector CDouble -- ^ Initial conditions |
233 | -> V.Vector CDouble -- ^ Desired solution times | 233 | -> V.Vector CDouble -- ^ Desired solution times |
234 | -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution | 234 | -> Either (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or |
235 | -- solution and diagnostics | ||
235 | solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize | 236 | solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize |
236 | jacH (aTols, rTol) fun f0 ts = | 237 | jacH (aTols, rTol) fun f0 ts = |
237 | unsafePerformIO $ do | 238 | unsafePerformIO $ do |
@@ -367,22 +368,16 @@ solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize | |||
367 | } | 368 | } |
368 | 369 | ||
369 | /* Main time-stepping loop: calls CVode to perform the integration */ | 370 | /* Main time-stepping loop: calls CVode to perform the integration */ |
370 | /* Stops when the final time has been reached */ | 371 | /* Stops when the final time has been reached */ |
371 | for (i = 1; i < $(int nTs); i++) { | 372 | for (i = 1; i < $(int nTs); i++) { |
372 | 373 | ||
373 | flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[i], y, &t, CV_NORMAL); /* call integrator */ | 374 | flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[i], y, &t, CV_NORMAL); /* call integrator */ |
374 | if (check_flag(&flag, "CVode", 1)) break; | 375 | if (check_flag(&flag, "CVode solver failure, stopping integration", 1)) return 1; |
375 | 376 | ||
376 | /* Store the results for Haskell */ | 377 | /* Store the results for Haskell */ |
377 | for (j = 0; j < NEQ; j++) { | 378 | for (j = 0; j < NEQ; j++) { |
378 | ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j); | 379 | ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j); |
379 | } | 380 | } |
380 | |||
381 | /* unsuccessful solve: break */ | ||
382 | if (flag < 0) { | ||
383 | fprintf(stderr,"Solver failure, stopping integration\n"); | ||
384 | break; | ||
385 | } | ||
386 | } | 381 | } |
387 | 382 | ||
388 | /* Get some final statistics on how the solve progressed */ | 383 | /* Get some final statistics on how the solve progressed */ |
@@ -434,23 +429,23 @@ solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize | |||
434 | 429 | ||
435 | return flag; | 430 | return flag; |
436 | } |] | 431 | } |] |
432 | preD <- V.freeze diagMut | ||
433 | let d = SundialsDiagnostics (fromIntegral $ preD V.!0) | ||
434 | (fromIntegral $ preD V.!1) | ||
435 | (fromIntegral $ preD V.!2) | ||
436 | (fromIntegral $ preD V.!3) | ||
437 | (fromIntegral $ preD V.!4) | ||
438 | (fromIntegral $ preD V.!5) | ||
439 | (fromIntegral $ preD V.!6) | ||
440 | (fromIntegral $ preD V.!7) | ||
441 | (fromIntegral $ preD V.!8) | ||
442 | (fromIntegral $ preD V.!9) | ||
443 | m <- V.freeze qMatMut | ||
437 | if res == 0 | 444 | if res == 0 |
438 | then do | 445 | then do |
439 | preD <- V.freeze diagMut | ||
440 | let d = SundialsDiagnostics (fromIntegral $ preD V.!0) | ||
441 | (fromIntegral $ preD V.!1) | ||
442 | (fromIntegral $ preD V.!2) | ||
443 | (fromIntegral $ preD V.!3) | ||
444 | (fromIntegral $ preD V.!4) | ||
445 | (fromIntegral $ preD V.!5) | ||
446 | (fromIntegral $ preD V.!6) | ||
447 | (fromIntegral $ preD V.!7) | ||
448 | (fromIntegral $ preD V.!8) | ||
449 | (fromIntegral $ preD V.!9) | ||
450 | m <- V.freeze qMatMut | ||
451 | return $ Right (m, d) | 446 | return $ Right (m, d) |
452 | else do | 447 | else do |
453 | return $ Left res | 448 | return $ Left (m, res) |
454 | 449 | ||
455 | -- | Adaptive step-size control | 450 | -- | Adaptive step-size control |
456 | -- functions. | 451 | -- functions. |