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.hs153
1 files changed, 25 insertions, 128 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
index abe1bfe..d7a2b53 100644
--- a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
@@ -61,46 +61,6 @@
61-- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) 61-- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
62-- @ 62-- @
63-- 63--
64-- KVAERNO_4_2_3
65--
66-- \[
67-- \begin{array}{c|cccc}
68-- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
69-- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\
70-- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
71-- 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
72-- \hline
73-- & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
74-- & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
75-- \end{array}
76-- \]
77--
78-- SDIRK_2_1_2
79--
80-- \[
81-- \begin{array}{c|cc}
82-- 1.0 & 1.0 & 0.0 \\
83-- 0.0 & -1.0 & 1.0 \\
84-- \hline
85-- & 0.5 & 0.5 \\
86-- & 1.0 & 0.0 \\
87-- \end{array}
88-- \]
89--
90-- SDIRK_5_3_4
91--
92-- \[
93-- \begin{array}{c|ccccc}
94-- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\
95-- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\
96-- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\
97-- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\
98-- 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
99-- \hline
100-- & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
101-- & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\
102-- \end{array}
103-- \]
104----------------------------------------------------------------------------- 64-----------------------------------------------------------------------------
105module Numeric.Sundials.CVode.ODE ( odeSolve 65module Numeric.Sundials.CVode.ODE ( odeSolve
106 , odeSolveV 66 , odeSolveV
@@ -109,7 +69,6 @@ module Numeric.Sundials.CVode.ODE ( odeSolve
109 , ODEMethod(..) 69 , ODEMethod(..)
110 , StepControl(..) 70 , StepControl(..)
111 , Jacobian 71 , Jacobian
112 , SundialsDiagnostics(..)
113 ) where 72 ) where
114 73
115import qualified Language.C.Inline as C 74import qualified Language.C.Inline as C
@@ -118,13 +77,10 @@ import qualified Language.C.Inline.Unsafe as CU
118import Data.Monoid ((<>)) 77import Data.Monoid ((<>))
119import Data.Maybe (isJust) 78import Data.Maybe (isJust)
120 79
121import Foreign.C.Types 80import Foreign.C.Types (CDouble, CInt, CLong)
122import Foreign.Ptr (Ptr) 81import Foreign.Ptr (Ptr)
123import Foreign.ForeignPtr (newForeignPtr_)
124import Foreign.Storable (Storable)
125 82
126import qualified Data.Vector.Storable as V 83import qualified Data.Vector.Storable as V
127import qualified Data.Vector.Storable.Mutable as VM
128 84
129import Data.Coerce (coerce) 85import Data.Coerce (coerce)
130import System.IO.Unsafe (unsafePerformIO) 86import System.IO.Unsafe (unsafePerformIO)
@@ -136,7 +92,7 @@ import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
136 92
137import qualified Types as T 93import qualified Types as T
138import Arkode (cV_ADAMS, cV_BDF) 94import Arkode (cV_ADAMS, cV_BDF)
139import qualified Arkode as B 95import qualified Numeric.Sundials.ODEOpts as SO
140 96
141 97
142C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx) 98C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
@@ -155,65 +111,6 @@ C.include "../../../helpers.h"
155C.include "Arkode_hsc.h" 111C.include "Arkode_hsc.h"
156 112
157 113
158getDataFromContents :: Int -> Ptr T.SunVector -> IO (V.Vector CDouble)
159getDataFromContents len ptr = do
160 qtr <- B.getContentPtr ptr
161 rtr <- B.getData qtr
162 vectorFromC len rtr
163
164-- FIXME: Potentially an instance of Storable
165_getMatrixDataFromContents :: Ptr T.SunMatrix -> IO T.SunMatrix
166_getMatrixDataFromContents ptr = do
167 qtr <- B.getContentMatrixPtr ptr
168 rs <- B.getNRows qtr
169 cs <- B.getNCols qtr
170 rtr <- B.getMatrixData qtr
171 vs <- vectorFromC (fromIntegral $ rs * cs) rtr
172 return $ T.SunMatrix { T.rows = rs, T.cols = cs, T.vals = vs }
173
174putMatrixDataFromContents :: T.SunMatrix -> Ptr T.SunMatrix -> IO ()
175putMatrixDataFromContents mat ptr = do
176 let rs = T.rows mat
177 cs = T.cols mat
178 vs = T.vals mat
179 qtr <- B.getContentMatrixPtr ptr
180 B.putNRows rs qtr
181 B.putNCols cs qtr
182 rtr <- B.getMatrixData qtr
183 vectorToC vs (fromIntegral $ rs * cs) rtr
184-- FIXME: END
185
186putDataInContents :: Storable a => V.Vector a -> Int -> Ptr b -> IO ()
187putDataInContents vec len ptr = do
188 qtr <- B.getContentPtr ptr
189 rtr <- B.getData qtr
190 vectorToC vec len rtr
191
192-- Utils
193
194vectorFromC :: Storable a => Int -> Ptr a -> IO (V.Vector a)
195vectorFromC len ptr = do
196 ptr' <- newForeignPtr_ ptr
197 V.freeze $ VM.unsafeFromForeignPtr0 ptr' len
198
199vectorToC :: Storable a => V.Vector a -> Int -> Ptr a -> IO ()
200vectorToC vec len ptr = do
201 ptr' <- newForeignPtr_ ptr
202 V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
203
204data SundialsDiagnostics = SundialsDiagnostics {
205 aRKodeGetNumSteps :: Int
206 , aRKodeGetNumStepAttempts :: Int
207 , aRKodeGetNumRhsEvals_fe :: Int
208 , aRKodeGetNumRhsEvals_fi :: Int
209 , aRKodeGetNumLinSolvSetups :: Int
210 , aRKodeGetNumErrTestFails :: Int
211 , aRKodeGetNumNonlinSolvIters :: Int
212 , aRKodeGetNumNonlinSolvConvFails :: Int
213 , aRKDlsGetNumJacEvals :: Int
214 , aRKDlsGetNumRhsEvals :: Int
215 } deriving Show
216
217type Jacobian = Double -> Vector Double -> Matrix Double 114type Jacobian = Double -> Vector Double -> Matrix Double
218 115
219-- | Stepping functions 116-- | Stepping functions
@@ -243,7 +140,7 @@ odeSolveV
243 -> Vector Double -- ^ desired solution times 140 -> Vector Double -- ^ desired solution times
244 -> Matrix Double -- ^ solution 141 -> Matrix Double -- ^ solution
245odeSolveV meth hi epsAbs epsRel f y0 ts = 142odeSolveV meth hi epsAbs epsRel f y0 ts =
246 case odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts of 143 case odeSolveVWith' meth (X epsAbs epsRel) hi g y0 ts of
247 Left c -> error $ show c -- FIXME 144 Left c -> error $ show c -- FIXME
248 -- FIXME: Can we do better than using lists? 145 -- FIXME: Can we do better than using lists?
249 Right (v, _d) -> (nR >< nC) (V.toList v) 146 Right (v, _d) -> (nR >< nC) (V.toList v)
@@ -263,7 +160,7 @@ odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y
263 -> Matrix Double -- ^ solution 160 -> Matrix Double -- ^ solution
264odeSolve f y0 ts = 161odeSolve f y0 ts =
265 -- FIXME: These tolerances are different from the ones in GSL 162 -- FIXME: These tolerances are different from the ones in GSL
266 case odeSolveVWith BDF (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts) of 163 case odeSolveVWith' BDF (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts) of
267 Left c -> error $ show c -- FIXME 164 Left c -> error $ show c -- FIXME
268 Right (v, _d) -> (nR >< nC) (V.toList v) 165 Right (v, _d) -> (nR >< nC) (V.toList v)
269 where 166 where
@@ -272,7 +169,7 @@ odeSolve f y0 ts =
272 nC = length y0 169 nC = length y0
273 g t x0 = V.fromList $ f t (V.toList x0) 170 g t x0 = V.fromList $ f t (V.toList x0)
274 171
275odeSolveVWith' :: 172odeSolveVWith ::
276 ODEMethod 173 ODEMethod
277 -> StepControl 174 -> StepControl
278 -> Maybe Double -- ^ initial step size - by default, ARKode 175 -> Maybe Double -- ^ initial step size - by default, ARKode
@@ -285,15 +182,15 @@ odeSolveVWith' ::
285 -> V.Vector Double -- ^ Initial conditions 182 -> V.Vector Double -- ^ Initial conditions
286 -> V.Vector Double -- ^ Desired solution times 183 -> V.Vector Double -- ^ Desired solution times
287 -> Matrix Double -- ^ Error code or solution 184 -> Matrix Double -- ^ Error code or solution
288odeSolveVWith' method control initStepSize f y0 tt = 185odeSolveVWith method control initStepSize f y0 tt =
289 case odeSolveVWith method control initStepSize f y0 tt of 186 case odeSolveVWith' method control initStepSize f y0 tt of
290 Left c -> error $ show c -- FIXME 187 Left c -> error $ show c -- FIXME
291 Right (v, _d) -> (nR >< nC) (V.toList v) 188 Right (v, _d) -> (nR >< nC) (V.toList v)
292 where 189 where
293 nR = V.length tt 190 nR = V.length tt
294 nC = V.length y0 191 nC = V.length y0
295 192
296odeSolveVWith :: 193odeSolveVWith' ::
297 ODEMethod 194 ODEMethod
298 -> StepControl 195 -> StepControl
299 -> Maybe Double -- ^ initial step size - by default, ARKode 196 -> Maybe Double -- ^ initial step size - by default, ARKode
@@ -305,8 +202,8 @@ odeSolveVWith ::
305 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 202 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
306 -> V.Vector Double -- ^ Initial conditions 203 -> V.Vector Double -- ^ Initial conditions
307 -> V.Vector Double -- ^ Desired solution times 204 -> V.Vector Double -- ^ Desired solution times
308 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution 205 -> Either Int ((V.Vector Double), SO.SundialsDiagnostics) -- ^ Error code or solution
309odeSolveVWith method control initStepSize f y0 tt = 206odeSolveVWith' method control initStepSize f y0 tt =
310 case solveOdeC (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) 207 case solveOdeC (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
311 (coerce f) (coerce y0) (coerce tt) of 208 (coerce f) (coerce y0) (coerce tt) of
312 Left c -> Left $ fromIntegral c 209 Left c -> Left $ fromIntegral c
@@ -335,7 +232,7 @@ solveOdeC ::
335 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 232 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
336 -> V.Vector CDouble -- ^ Initial conditions 233 -> V.Vector CDouble -- ^ Initial conditions
337 -> V.Vector CDouble -- ^ Desired solution times 234 -> V.Vector CDouble -- ^ Desired solution times
338 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution 235 -> Either CInt ((V.Vector CDouble), SO.SundialsDiagnostics) -- ^ Error code or solution
339solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO $ do 236solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO $ do
340 237
341 let isInitStepSize :: CInt 238 let isInitStepSize :: CInt
@@ -366,9 +263,9 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
366 funIO x y f _ptr = do 263 funIO x y f _ptr = do
367 -- Convert the pointer we get from C (y) to a vector, and then 264 -- Convert the pointer we get from C (y) to a vector, and then
368 -- apply the user-supplied function. 265 -- apply the user-supplied function.
369 fImm <- fun x <$> getDataFromContents dim y 266 fImm <- fun x <$> SO.getDataFromContents dim y
370 -- Fill in the provided pointer with the resulting vector. 267 -- Fill in the provided pointer with the resulting vector.
371 putDataInContents fImm dim f 268 SO.putDataInContents fImm dim f
372 -- FIXME: I don't understand what this comment means 269 -- FIXME: I don't understand what this comment means
373 -- Unsafe since the function will be called many times. 270 -- Unsafe since the function will be called many times.
374 [CU.exp| int{ 0 } |] 271 [CU.exp| int{ 0 } |]
@@ -380,8 +277,8 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
380 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do 277 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
381 case jacH of 278 case jacH of
382 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined" 279 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
383 Just jacI -> do j <- jacI t <$> getDataFromContents dim y 280 Just jacI -> do j <- jacI t <$> SO.getDataFromContents dim y
384 putMatrixDataFromContents j jacS 281 SO.putMatrixDataFromContents j jacS
385 -- FIXME: I don't understand what this comment means 282 -- FIXME: I don't understand what this comment means
386 -- Unsafe since the function will be called many times. 283 -- Unsafe since the function will be called many times.
387 [CU.exp| int{ 0 } |] 284 [CU.exp| int{ 0 } |]
@@ -541,16 +438,16 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
541 if res == 0 438 if res == 0
542 then do 439 then do
543 preD <- V.freeze diagMut 440 preD <- V.freeze diagMut
544 let d = SundialsDiagnostics (fromIntegral $ preD V.!0) 441 let d = SO.SundialsDiagnostics (fromIntegral $ preD V.!0)
545 (fromIntegral $ preD V.!1) 442 (fromIntegral $ preD V.!1)
546 (fromIntegral $ preD V.!2) 443 (fromIntegral $ preD V.!2)
547 (fromIntegral $ preD V.!3) 444 (fromIntegral $ preD V.!3)
548 (fromIntegral $ preD V.!4) 445 (fromIntegral $ preD V.!4)
549 (fromIntegral $ preD V.!5) 446 (fromIntegral $ preD V.!5)
550 (fromIntegral $ preD V.!6) 447 (fromIntegral $ preD V.!6)
551 (fromIntegral $ preD V.!7) 448 (fromIntegral $ preD V.!7)
552 (fromIntegral $ preD V.!8) 449 (fromIntegral $ preD V.!8)
553 (fromIntegral $ preD V.!9) 450 (fromIntegral $ preD V.!9)
554 m <- V.freeze qMatMut 451 m <- V.freeze qMatMut
555 return $ Right (m, d) 452 return $ Right (m, d)
556 else do 453 else do
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 05:28:49 +0000