summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-04-05 18:28:56 +0100
committerDominic Steinitz <dominic@steinitz.org>2018-04-05 18:28:56 +0100
commit59a0413a83a9bcee93e3f0761cae6fdda2a98933 (patch)
treed706524e98e62bae13283af05385d126bfb6890c /packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
parent6fcc1b01cecc88f1a8eb1608667368c7e72048aa (diff)
Get closer to the hmatrix-gsl interface
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs')
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs42
1 files changed, 24 insertions, 18 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 0973c82..4270a13 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -185,14 +185,22 @@ data SundialsDiagnostics = SundialsDiagnostics {
185type Jacobian = Double -> Vector Double -> Matrix Double 185type Jacobian = Double -> Vector Double -> Matrix Double
186 186
187-- | Stepping functions 187-- | Stepping functions
188data ODEMethod = SDIRK_2_1_2 Jacobian 188data ODEMethod = SDIRK_2_1_2 Jacobian
189 | KVAERNO_4_2_3 Jacobian 189 | KVAERNO_4_2_3 Jacobian
190 | SDIRK_5_3_4 Jacobian 190 | SDIRK_5_3_4 Jacobian
191 | SDIRK_5_3_4'
191 192
192getMethod :: ODEMethod -> Int 193getMethod :: ODEMethod -> Int
193getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2 194getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2
194getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3 195getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3
195getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4 196getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4
197getMethod (SDIRK_5_3_4' ) = sDIRK_5_3_4
198
199getJacobian :: ODEMethod -> Maybe Jacobian
200getJacobian (SDIRK_2_1_2 j) = Just j
201getJacobian (KVAERNO_4_2_3 j) = Just j
202getJacobian (SDIRK_5_3_4 j) = Just j
203getJacobian (SDIRK_5_3_4' ) = Nothing
196 204
197-- | A version of 'odeSolveVWith' with reasonable default step control. 205-- | A version of 'odeSolveVWith' with reasonable default step control.
198odeSolveV 206odeSolveV
@@ -215,7 +223,7 @@ odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y
215 -> Vector Double -- ^ desired solution times 223 -> Vector Double -- ^ desired solution times
216 -> Matrix Double -- ^ solution 224 -> Matrix Double -- ^ solution
217odeSolve f y0 ts = 225odeSolve f y0 ts =
218 case odeSolveVWith (SDIRK_5_3_4 undefined) Nothing 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of 226 case odeSolveVWith SDIRK_5_3_4' 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of
219 Left c -> error $ show c -- FIXME 227 Left c -> error $ show c -- FIXME
220 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v) 228 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v)
221 where 229 where
@@ -225,13 +233,12 @@ odeSolve f y0 ts =
225 g t x0 = V.fromList $ f t (V.toList x0) 233 g t x0 = V.fromList $ f t (V.toList x0)
226 234
227odeSolve' :: ODEMethod 235odeSolve' :: ODEMethod
228 -> (Double -> Vector Double -> Matrix Double)
229 -> (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 236 -> (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
230 -> [Double] -- ^ initial conditions 237 -> [Double] -- ^ initial conditions
231 -> Vector Double -- ^ desired solution times 238 -> Vector Double -- ^ desired solution times
232 -> Matrix Double -- ^ solution 239 -> Matrix Double -- ^ solution
233odeSolve' method jac f y0 ts = 240odeSolve' method f y0 ts =
234 case odeSolveVWith method (pure jac') 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of 241 case odeSolveVWith method 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of
235 Left c -> error $ show c -- FIXME 242 Left c -> error $ show c -- FIXME
236 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v) 243 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v)
237 where 244 where
@@ -239,30 +246,29 @@ odeSolve' method jac f y0 ts =
239 nR = length us 246 nR = length us
240 nC = length y0 247 nC = length y0
241 g t x0 = V.fromList $ f t (V.toList x0) 248 g t x0 = V.fromList $ f t (V.toList x0)
242 jac' t v = foo $ jac t (V.fromList $ toList v)
243 foo m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
244 where
245 nr = fromIntegral $ rows m
246 nc = fromIntegral $ cols m
247 vs = V.fromList $ map coerce $ concat $ toLists m
248 249
249odeSolveVWith :: 250odeSolveVWith ::
250 ODEMethod 251 ODEMethod
251 -> (Maybe (Double -> V.Vector Double -> T.SunMatrix))
252 -> Double 252 -> Double
253 -> Double 253 -> Double
254 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 254 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
255 -> V.Vector Double -- ^ Initial conditions 255 -> V.Vector Double -- ^ Initial conditions
256 -> V.Vector Double -- ^ Desired solution times 256 -> V.Vector Double -- ^ Desired solution times
257 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution 257 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution
258odeSolveVWith method jac relTol absTol f y0 tt = 258odeSolveVWith method relTol absTol f y0 tt =
259 case solveOdeC (fromIntegral $ getMethod method) jacH (CDouble relTol) (CDouble absTol) 259 case solveOdeC (fromIntegral $ getMethod method) jacH (CDouble relTol) (CDouble absTol)
260 (coerce f) (coerce y0) (coerce tt) of 260 (coerce f) (coerce y0) (coerce tt) of
261 Left c -> Left $ fromIntegral c 261 Left c -> Left $ fromIntegral c
262 Right (v, d) -> Right (coerce v, d) 262 Right (v, d) -> Right (coerce v, d)
263 where 263 where
264 jacH :: Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix) 264 -- FIXME: Can we do better than going via a list?
265 jacH = fmap (\g -> (\t v -> g (coerce t) (coerce v))) jac 265 jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce $ V.fromList $ toList v)) $
266 getJacobian method
267 matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
268 where
269 nr = fromIntegral $ rows m
270 nc = fromIntegral $ cols m
271 vs = V.fromList $ map coerce $ concat $ toLists m
266 272
267solveOdeC :: 273solveOdeC ::
268 CInt -> 274 CInt ->