summaryrefslogtreecommitdiff
path: root/packages/sundials/src
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials/src')
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs47
-rw-r--r--packages/sundials/src/Numeric/Sundials/CVode/ODE.hs45
2 files changed, 40 insertions, 52 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index fafc237..48ac887 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -428,7 +428,7 @@ odeSolveVWith ::
428 -> Matrix Double -- ^ Error code or solution 428 -> Matrix Double -- ^ Error code or solution
429odeSolveVWith method control initStepSize f y0 tt = 429odeSolveVWith method control initStepSize f y0 tt =
430 case odeSolveVWith' opts method control initStepSize f y0 tt of 430 case odeSolveVWith' opts method control initStepSize f y0 tt of
431 Left c -> error $ show c -- FIXME 431 Left (c, _v) -> error $ show c -- FIXME
432 Right (v, _d) -> v 432 Right (v, _d) -> v
433 where 433 where
434 opts = ODEOpts { maxNumSteps = 10000 434 opts = ODEOpts { maxNumSteps = 10000
@@ -452,14 +452,14 @@ odeSolveVWith' ::
452 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 452 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
453 -> V.Vector Double -- ^ Initial conditions 453 -> V.Vector Double -- ^ Initial conditions
454 -> V.Vector Double -- ^ Desired solution times 454 -> V.Vector Double -- ^ Desired solution times
455 -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution 455 -> Either (Matrix Double, Int) (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
456odeSolveVWith' opts method control initStepSize f y0 tt = 456odeSolveVWith' opts method control initStepSize f y0 tt =
457 case solveOdeC (fromIntegral $ maxFail opts) 457 case solveOdeC (fromIntegral $ maxFail opts)
458 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) 458 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
459 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) 459 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
460 (coerce f) (coerce y0) (coerce tt) of 460 (coerce f) (coerce y0) (coerce tt) of
461 Left c -> Left $ fromIntegral c 461 Left (v, c) -> Left (reshape l (coerce v), fromIntegral c)
462 Right (v, d) -> Right (reshape l (coerce v), d) 462 Right (v, d) -> Right (reshape l (coerce v), d)
463 where 463 where
464 l = size y0 464 l = size y0
465 scise (X aTol rTol) = coerce (V.replicate l aTol, rTol) 465 scise (X aTol rTol) = coerce (V.replicate l aTol, rTol)
@@ -487,7 +487,8 @@ solveOdeC ::
487 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 487 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
488 -> V.Vector CDouble -- ^ Initial conditions 488 -> V.Vector CDouble -- ^ Initial conditions
489 -> V.Vector CDouble -- ^ Desired solution times 489 -> V.Vector CDouble -- ^ Desired solution times
490 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution 490 -> Either (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or
491 -- solution and diagnostics
491solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize 492solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
492 jacH (aTols, rTol) fun f0 ts = unsafePerformIO $ do 493 jacH (aTols, rTol) fun f0 ts = unsafePerformIO $ do
493 494
@@ -506,8 +507,6 @@ solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
506 nEq = fromIntegral dim 507 nEq = fromIntegral dim
507 nTs :: CInt 508 nTs :: CInt
508 nTs = fromIntegral $ V.length ts 509 nTs = fromIntegral $ V.length ts
509 -- FIXME: I believe this gets taken from the ghc heap and so should
510 -- be subject to garbage collection.
511 quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs)) 510 quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
512 qMatMut <- V.thaw quasiMatrixRes 511 qMatMut <- V.thaw quasiMatrixRes
513 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME 512 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
@@ -644,18 +643,12 @@ solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
644 for (i = 1; i < $(int nTs); i++) { 643 for (i = 1; i < $(int nTs); i++) {
645 644
646 flag = ARKode(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */ 645 flag = ARKode(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */
647 if (check_flag(&flag, "ARKode", 1)) break; 646 if (check_flag(&flag, "ARKode solver failure, stopping integration", 1)) return 1;
648 647
649 /* Store the results for Haskell */ 648 /* Store the results for Haskell */
650 for (j = 0; j < NEQ; j++) { 649 for (j = 0; j < NEQ; j++) {
651 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j); 650 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
652 } 651 }
653
654 /* unsuccessful solve: break */
655 if (flag < 0) {
656 fprintf(stderr,"Solver failure, stopping integration\n");
657 break;
658 }
659 } 652 }
660 653
661 /* Get some final statistics on how the solve progressed */ 654 /* Get some final statistics on how the solve progressed */
@@ -706,23 +699,23 @@ solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
706 699
707 return flag; 700 return flag;
708 } |] 701 } |]
702 preD <- V.freeze diagMut
703 let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
704 (fromIntegral $ preD V.!1)
705 (fromIntegral $ preD V.!2)
706 (fromIntegral $ preD V.!3)
707 (fromIntegral $ preD V.!4)
708 (fromIntegral $ preD V.!5)
709 (fromIntegral $ preD V.!6)
710 (fromIntegral $ preD V.!7)
711 (fromIntegral $ preD V.!8)
712 (fromIntegral $ preD V.!9)
713 m <- V.freeze qMatMut
709 if res == 0 714 if res == 0
710 then do 715 then do
711 preD <- V.freeze diagMut
712 let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
713 (fromIntegral $ preD V.!1)
714 (fromIntegral $ preD V.!2)
715 (fromIntegral $ preD V.!3)
716 (fromIntegral $ preD V.!4)
717 (fromIntegral $ preD V.!5)
718 (fromIntegral $ preD V.!6)
719 (fromIntegral $ preD V.!7)
720 (fromIntegral $ preD V.!8)
721 (fromIntegral $ preD V.!9)
722 m <- V.freeze qMatMut
723 return $ Right (m, d) 716 return $ Right (m, d)
724 else do 717 else do
725 return $ Left res 718 return $ Left (m, res)
726 719
727data ButcherTable = ButcherTable { am :: Matrix Double 720data ButcherTable = ButcherTable { am :: Matrix Double
728 , cv :: Vector Double 721 , cv :: Vector Double
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
173odeSolveVWith method control initStepSize f y0 tt = 173odeSolveVWith 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
200odeSolveVWith' opts method control initStepSize f y0 tt = 200odeSolveVWith' 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
235solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize 236solveOdeC 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.
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-14 09:51:48 +0000