diff options
Diffstat (limited to 'packages/sundials/src')
-rw-r--r-- | packages/sundials/src/Main.hs | 45 | ||||
-rw-r--r-- | packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs | 146 |
2 files changed, 124 insertions, 67 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs index d22fafa..b65d668 100644 --- a/packages/sundials/src/Main.hs +++ b/packages/sundials/src/Main.hs | |||
@@ -94,44 +94,39 @@ kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double | |||
94 | kSaxis xs = P.r2Axis &~ do | 94 | kSaxis xs = P.r2Axis &~ do |
95 | P.linePlot' xs | 95 | P.linePlot' xs |
96 | 96 | ||
97 | butcherTableauTex :: (Show a, Element a) => Matrix a -> String | 97 | butcherTableauTex :: ButcherTable -> String |
98 | butcherTableauTex m = render $ | 98 | butcherTableauTex (ButcherTable m c b b2) = |
99 | vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}") | 99 | render $ |
100 | , us | 100 | vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}") |
101 | , text "\\end{array}" | 101 | , us |
102 | ] | 102 | , text "\\hline" |
103 | , text bs <+> text "\\\\" | ||
104 | , text b2s <+> text "\\\\" | ||
105 | , text "\\end{array}" | ||
106 | ] | ||
103 | where | 107 | where |
104 | n = rows m | 108 | n = rows m |
105 | rs = toLists m | 109 | rs = toLists m |
106 | ss = map (\r -> intercalate " & " $ map show r) rs | 110 | ss = map (\r -> intercalate " & " $ map show r) rs |
107 | ts = zipWith (\i r -> "c_" ++ show i ++ " & " ++ r) [1..n] ss | 111 | ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss |
108 | us = vcat $ map (\r -> text r <+> text "\\\\") ts | 112 | us = vcat $ map (\r -> text r <+> text "\\\\") ts |
113 | bs = " & " ++ (intercalate " & " $ map show $ toList b) | ||
114 | b2s = " & " ++ (intercalate " & " $ map show $ toList b2) | ||
109 | 115 | ||
110 | main :: IO () | 116 | main :: IO () |
111 | main = do | 117 | main = do |
112 | -- $$ | 118 | |
113 | -- \begin{array}{c|cccc} | 119 | let res = butcherTable (SDIRK_2_1_2 undefined) |
114 | -- c_1 & a_{11} & a_{12}& \dots & a_{1s}\\ | ||
115 | -- c_2 & a_{21} & a_{22}& \dots & a_{2s}\\ | ||
116 | -- \vdots & \vdots & \vdots& \ddots& \vdots\\ | ||
117 | -- c_s & a_{s1} & a_{s2}& \dots & a_{ss} \\ | ||
118 | -- \hline | ||
119 | -- & b_1 & b_2 & \dots & b_s\\ | ||
120 | -- & b^*_1 & b^*_2 & \dots & b^*_s\\ | ||
121 | -- \end{array} | ||
122 | -- $$ | ||
123 | |||
124 | let res = btGet (SDIRK_2_1_2 undefined) | ||
125 | putStrLn $ show res | 120 | putStrLn $ show res |
126 | putStrLn $ butcherTableauTex $ fst res | 121 | putStrLn $ butcherTableauTex res |
127 | 122 | ||
128 | let res = btGet (KVAERNO_4_2_3 undefined) | 123 | let res = butcherTable (KVAERNO_4_2_3 undefined) |
129 | putStrLn $ show res | 124 | putStrLn $ show res |
130 | putStrLn $ butcherTableauTex $ fst res | 125 | putStrLn $ butcherTableauTex res |
131 | 126 | ||
132 | let res = btGet (SDIRK_5_3_4 undefined) | 127 | let res = butcherTable (SDIRK_5_3_4 undefined) |
133 | putStrLn $ show res | 128 | putStrLn $ show res |
134 | putStrLn $ butcherTableauTex $ fst res | 129 | putStrLn $ butcherTableauTex res |
135 | 130 | ||
136 | let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | 131 | let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) |
137 | renderRasterific "diagrams/brusselator.png" | 132 | renderRasterific "diagrams/brusselator.png" |
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs index b6a59e2..67378cc 100644 --- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs +++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs | |||
@@ -21,30 +21,57 @@ | |||
21 | -- | 21 | -- |
22 | -- A simple example: | 22 | -- A simple example: |
23 | -- | 23 | -- |
24 | -- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>> | ||
25 | -- | ||
24 | -- @ | 26 | -- @ |
25 | -- import Numeric.Sundials.ARKode | 27 | -- import Numeric.Sundials.ARKode.ODE |
26 | -- import Numeric.LinearAlgebra | 28 | -- import Numeric.LinearAlgebra |
27 | -- import Graphics.Plot(mplot) | ||
28 | -- | 29 | -- |
29 | -- xdot t [x,v] = [v, -0.95*x - 0.1*v] | 30 | -- import Plots as P |
31 | -- import qualified Diagrams.Prelude as D | ||
32 | -- import Diagrams.Backend.Rasterific | ||
30 | -- | 33 | -- |
31 | -- ts = linspace 100 (0,20 :: Double) | 34 | -- brusselator :: Double -> [Double] -> [Double] |
35 | -- brusselator _t x = [ a - (w + 1) * u + v * u * u | ||
36 | -- , w * u - v * u * u | ||
37 | -- , (b - w) / eps - w * u | ||
38 | -- ] | ||
39 | -- where | ||
40 | -- a = 1.0 | ||
41 | -- b = 3.5 | ||
42 | -- eps = 5.0e-6 | ||
43 | -- u = x !! 0 | ||
44 | -- v = x !! 1 | ||
45 | -- w = x !! 2 | ||
32 | -- | 46 | -- |
33 | -- sol = odeSolve xdot [10,0] ts | 47 | -- lSaxis :: [[Double]] -> P.Axis B D.V2 Double |
48 | -- lSaxis xs = P.r2Axis &~ do | ||
49 | -- let ts = xs!!0 | ||
50 | -- us = xs!!1 | ||
51 | -- vs = xs!!2 | ||
52 | -- ws = xs!!3 | ||
53 | -- P.linePlot' $ zip ts us | ||
54 | -- P.linePlot' $ zip ts vs | ||
55 | -- P.linePlot' $ zip ts ws | ||
34 | -- | 56 | -- |
35 | -- main = mplot (ts : toColumns sol) | 57 | -- main = do |
58 | -- let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | ||
59 | -- renderRasterific "diagrams/brusselator.png" | ||
60 | -- (D.dims2D 500.0 500.0) | ||
61 | -- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) | ||
36 | -- @ | 62 | -- @ |
37 | -- | 63 | -- |
38 | -- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>> | ||
39 | -- | ||
40 | -- KVAERNO_4_2_3 | 64 | -- KVAERNO_4_2_3 |
41 | -- | 65 | -- |
42 | -- \[ | 66 | -- \[ |
43 | -- \begin{array}{c|cccc} | 67 | -- \begin{array}{c|cccc} |
44 | -- c_1 & 0.0 & 0.0 & 0.0 & 0.0 \\ | 68 | -- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ |
45 | -- c_2 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\ | 69 | -- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\ |
46 | -- c_3 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\ | 70 | -- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\ |
47 | -- c_4 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\ | 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 \\ | ||
48 | -- \end{array} | 75 | -- \end{array} |
49 | -- \] | 76 | -- \] |
50 | -- | 77 | -- |
@@ -52,8 +79,11 @@ | |||
52 | -- | 79 | -- |
53 | -- \[ | 80 | -- \[ |
54 | -- \begin{array}{c|cc} | 81 | -- \begin{array}{c|cc} |
55 | -- c_1 & 1.0 & 0.0 \\ | 82 | -- 1.0 & 1.0 & 0.0 \\ |
56 | -- c_2 & -1.0 & 1.0 \\ | 83 | -- 0.0 & -1.0 & 1.0 \\ |
84 | -- \hline | ||
85 | -- & 0.5 & 0.5 \\ | ||
86 | -- & 1.0 & 0.0 \\ | ||
57 | -- \end{array} | 87 | -- \end{array} |
58 | -- \] | 88 | -- \] |
59 | -- | 89 | -- |
@@ -61,20 +91,22 @@ | |||
61 | -- | 91 | -- |
62 | -- \[ | 92 | -- \[ |
63 | -- \begin{array}{c|ccccc} | 93 | -- \begin{array}{c|ccccc} |
64 | -- c_1 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\ | 94 | -- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\ |
65 | -- c_2 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\ | 95 | -- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\ |
66 | -- c_3 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\ | 96 | -- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\ |
67 | -- c_4 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\ | 97 | -- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\ |
68 | -- c_5 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\ | 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 \\ | ||
69 | -- \end{array} | 102 | -- \end{array} |
70 | -- \] | 103 | -- \] |
71 | ----------------------------------------------------------------------------- | 104 | ----------------------------------------------------------------------------- |
72 | module Numeric.Sundials.ARKode.ODE ( odeSolve | 105 | module Numeric.Sundials.ARKode.ODE ( odeSolve |
73 | , odeSolveV | 106 | , odeSolveV |
74 | , odeSolveVWith | 107 | , odeSolveVWith |
75 | , getButcherTable | 108 | , ButcherTable(..) |
76 | , getBT | 109 | , butcherTable |
77 | , btGet | ||
78 | , ODEMethod(..) | 110 | , ODEMethod(..) |
79 | , StepControl(..) | 111 | , StepControl(..) |
80 | , SundialsDiagnostics(..) | 112 | , SundialsDiagnostics(..) |
@@ -457,10 +489,11 @@ solveOdeC method jacH (absTols, relTol) fun f0 ts = unsafePerformIO $ do | |||
457 | ($vec-ptr:(long int *diagMut))[9] = ncfn; | 489 | ($vec-ptr:(long int *diagMut))[9] = ncfn; |
458 | 490 | ||
459 | /* Clean up and return */ | 491 | /* Clean up and return */ |
460 | N_VDestroy(y); /* Free y vector */ | 492 | N_VDestroy(y); /* Free y vector */ |
493 | N_VDestroy(tv); /* Free tv vector */ | ||
461 | ARKodeFree(&arkode_mem); /* Free integrator memory */ | 494 | ARKodeFree(&arkode_mem); /* Free integrator memory */ |
462 | SUNLinSolFree(LS); /* Free linear solver */ | 495 | SUNLinSolFree(LS); /* Free linear solver */ |
463 | SUNMatDestroy(A); /* Free A matrix */ | 496 | SUNMatDestroy(A); /* Free A matrix */ |
464 | 497 | ||
465 | return flag; | 498 | return flag; |
466 | } |] | 499 | } |] |
@@ -482,22 +515,43 @@ solveOdeC method jacH (absTols, relTol) fun f0 ts = unsafePerformIO $ do | |||
482 | else do | 515 | else do |
483 | return $ Left res | 516 | return $ Left res |
484 | 517 | ||
485 | btGet :: ODEMethod -> (Matrix Double, Vector Double) | 518 | data ButcherTable = ButcherTable { am :: Matrix Double |
486 | btGet method = case getBT method of | 519 | , cv :: Vector Double |
487 | Left c -> error $ show c -- FIXME | 520 | , bv :: Vector Double |
488 | Right ((v, w), sqp) -> ( subMatrix (0, 0) (s, s) $ | 521 | , b2v :: Vector Double |
489 | (B.arkSMax >< B.arkSMax) (V.toList v) | 522 | } |
490 | , subVector 0 s w) | 523 | deriving Show |
491 | where | 524 | |
492 | s = fromIntegral $ sqp V.! 0 | 525 | data ButcherTable' a = ButcherTable' { am' :: V.Vector a |
526 | , cv' :: V.Vector a | ||
527 | , bv' :: V.Vector a | ||
528 | , b2v' :: V.Vector a | ||
529 | } | ||
530 | deriving Show | ||
531 | |||
532 | butcherTable :: ODEMethod -> ButcherTable | ||
533 | butcherTable method = | ||
534 | case getBT method of | ||
535 | Left c -> error $ show c -- FIXME | ||
536 | Right (ButcherTable' v w x y, sqp) -> | ||
537 | ButcherTable { am = subMatrix (0, 0) (s, s) $ (B.arkSMax >< B.arkSMax) (V.toList v) | ||
538 | , cv = subVector 0 s w | ||
539 | , bv = subVector 0 s x | ||
540 | , b2v = subVector 0 s y | ||
541 | } | ||
542 | where | ||
543 | s = fromIntegral $ sqp V.! 0 | ||
493 | 544 | ||
494 | getBT :: ODEMethod -> Either Int ((V.Vector Double, V.Vector Double), V.Vector Int) | 545 | getBT :: ODEMethod -> Either Int (ButcherTable' Double, V.Vector Int) |
495 | getBT method = case getButcherTable method of | 546 | getBT method = case getButcherTable method of |
496 | Left c -> Left $ fromIntegral c | 547 | Left c -> |
497 | Right ((v, w), sqp) -> Right $ ((coerce v, coerce w), V.map fromIntegral sqp) | 548 | Left $ fromIntegral c |
549 | Right (ButcherTable' a b c d, sqp) -> | ||
550 | Right $ ( ButcherTable' (coerce a) (coerce b) (coerce c) (coerce d) | ||
551 | , V.map fromIntegral sqp ) | ||
498 | 552 | ||
499 | getButcherTable :: ODEMethod | 553 | getButcherTable :: ODEMethod |
500 | -> Either CInt ((V.Vector CDouble, V.Vector CDouble), V.Vector CInt) | 554 | -> Either CInt (ButcherTable' CDouble, V.Vector CInt) |
501 | getButcherTable method = unsafePerformIO $ do | 555 | getButcherTable method = unsafePerformIO $ do |
502 | -- ARKode seems to want an ODE in order to set and then get the | 556 | -- ARKode seems to want an ODE in order to set and then get the |
503 | -- Butcher tableau so here's one to keep it happy | 557 | -- Butcher tableau so here's one to keep it happy |
@@ -515,8 +569,12 @@ getButcherTable method = unsafePerformIO $ do | |||
515 | btSQPMut <- V.thaw btSQP | 569 | btSQPMut <- V.thaw btSQP |
516 | btAs :: V.Vector CDouble <- createVector (B.arkSMax * B.arkSMax) | 570 | btAs :: V.Vector CDouble <- createVector (B.arkSMax * B.arkSMax) |
517 | btAsMut <- V.thaw btAs | 571 | btAsMut <- V.thaw btAs |
518 | btCs :: V.Vector CDouble <- createVector B.arkSMax | 572 | btCs :: V.Vector CDouble <- createVector B.arkSMax |
519 | btCsMut <- V.thaw btCs | 573 | btBs :: V.Vector CDouble <- createVector B.arkSMax |
574 | btB2s :: V.Vector CDouble <- createVector B.arkSMax | ||
575 | btCsMut <- V.thaw btCs | ||
576 | btBsMut <- V.thaw btBs | ||
577 | btB2sMut <- V.thaw btB2s | ||
520 | let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt | 578 | let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt |
521 | funIO x y f _ptr = do | 579 | funIO x y f _ptr = do |
522 | fImm <- fun x <$> getDataFromContents dim y | 580 | fImm <- fun x <$> getDataFromContents dim y |
@@ -576,7 +634,9 @@ getButcherTable method = unsafePerformIO $ do | |||
576 | } | 634 | } |
577 | 635 | ||
578 | for (i = 0; i < s; i++) { | 636 | for (i = 0; i < s; i++) { |
579 | ($vec-ptr:(double *btCsMut))[i] = ci[i]; | 637 | ($vec-ptr:(double *btCsMut))[i] = ci[i]; |
638 | ($vec-ptr:(double *btBsMut))[i] = bi[i]; | ||
639 | ($vec-ptr:(double *btB2sMut))[i] = b2i[i]; | ||
580 | } | 640 | } |
581 | 641 | ||
582 | /* Clean up and return */ | 642 | /* Clean up and return */ |
@@ -590,7 +650,9 @@ getButcherTable method = unsafePerformIO $ do | |||
590 | x <- V.freeze btAsMut | 650 | x <- V.freeze btAsMut |
591 | y <- V.freeze btSQPMut | 651 | y <- V.freeze btSQPMut |
592 | z <- V.freeze btCsMut | 652 | z <- V.freeze btCsMut |
593 | return $ Right ((x, z), y) | 653 | u <- V.freeze btBsMut |
654 | v <- V.freeze btB2sMut | ||
655 | return $ Right (ButcherTable' { am' = x, cv' = z, bv' = u, b2v' = v }, y) | ||
594 | else do | 656 | else do |
595 | return $ Left res | 657 | return $ Left res |
596 | 658 | ||