summaryrefslogtreecommitdiff
path: root/packages/sundials
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials')
-rw-r--r--packages/sundials/ChangeLog.md5
-rw-r--r--packages/sundials/LICENSE30
-rw-r--r--packages/sundials/README.md8
-rw-r--r--packages/sundials/Setup.hs2
-rw-r--r--packages/sundials/diagrams/brusselator.pngbin0 -> 27362 bytes
-rw-r--r--packages/sundials/hmatrix-sundials.cabal61
-rw-r--r--packages/sundials/src/Main.hs186
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs903
-rw-r--r--packages/sundials/src/Numeric/Sundials/Arkode.hsc204
-rw-r--r--packages/sundials/src/Numeric/Sundials/CVode/ODE.hs476
-rw-r--r--packages/sundials/src/Numeric/Sundials/ODEOpts.hs32
-rw-r--r--packages/sundials/src/helpers.c44
-rw-r--r--packages/sundials/src/helpers.h9
13 files changed, 1960 insertions, 0 deletions
diff --git a/packages/sundials/ChangeLog.md b/packages/sundials/ChangeLog.md
new file mode 100644
index 0000000..7b15777
--- /dev/null
+++ b/packages/sundials/ChangeLog.md
@@ -0,0 +1,5 @@
1# Revision history for hmatrix-sundials
2
3## 0.1.0.0 -- 2018-04-21
4
5* First version. Released on an unsuspecting world. Just Runge-Kutta methods to start with.
diff --git a/packages/sundials/LICENSE b/packages/sundials/LICENSE
new file mode 100644
index 0000000..a162e98
--- /dev/null
+++ b/packages/sundials/LICENSE
@@ -0,0 +1,30 @@
1Copyright (c) 2018, Dominic Steinitz, Novadiscovery
2
3All rights reserved.
4
5Redistribution and use in source and binary forms, with or without
6modification, are permitted provided that the following conditions are met:
7
8 * Redistributions of source code must retain the above copyright
9 notice, this list of conditions and the following disclaimer.
10
11 * Redistributions in binary form must reproduce the above
12 copyright notice, this list of conditions and the following
13 disclaimer in the documentation and/or other materials provided
14 with the distribution.
15
16 * Neither the name of Dominic Steinitz nor the names of other
17 contributors may be used to endorse or promote products derived
18 from this software without specific prior written permission.
19
20THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/packages/sundials/README.md b/packages/sundials/README.md
new file mode 100644
index 0000000..2fac5c2
--- /dev/null
+++ b/packages/sundials/README.md
@@ -0,0 +1,8 @@
1Currently only an interface to the Runge-Kutta methods:
2[ARKode](https://computation.llnl.gov/projects/sundials/arkode)
3
4The interface is almost certainly going to change. Sundials gives a
5rich set of "combinators" for controlling the solution of your problem
6and reporting on how it performed. The idea is to initially mimic
7hmatrix-gsl and add extra, richer functions but ultimately upgrade the
8whole interface both for sundials and for gsl.
diff --git a/packages/sundials/Setup.hs b/packages/sundials/Setup.hs
new file mode 100644
index 0000000..9a994af
--- /dev/null
+++ b/packages/sundials/Setup.hs
@@ -0,0 +1,2 @@
1import Distribution.Simple
2main = defaultMain
diff --git a/packages/sundials/diagrams/brusselator.png b/packages/sundials/diagrams/brusselator.png
new file mode 100644
index 0000000..740cacb
--- /dev/null
+++ b/packages/sundials/diagrams/brusselator.png
Binary files differ
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
new file mode 100644
index 0000000..cd2be4e
--- /dev/null
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -0,0 +1,61 @@
1name: hmatrix-sundials
2version: 0.19.0.0
3synopsis: hmatrix interface to sundials
4description: An interface to the solving suite SUNDIALS. Currently, it
5 mimics the solving interface in hmstrix-gsl but
6 provides more diagnostic information and the
7 Butcher Tableaux (for Runge-Kutta methods).
8homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials
9license: BSD3
10license-file: LICENSE
11author: Dominic Steinitz
12maintainer: dominic@steinitz.org
13copyright: Dominic Steinitz 2018, Novadiscovery 2018
14category: Math
15build-type: Simple
16extra-source-files: ChangeLog.md, README.md, diagrams/*.png
17extra-doc-files: diagrams/*.png
18cabal-version: >=1.18
19
20
21library
22 build-depends: base >=4.10 && <4.11,
23 inline-c >=0.6 && <0.7,
24 vector >=0.12 && <0.13,
25 template-haskell >=2.12 && <2.13,
26 containers >=0.5 && <0.6,
27 hmatrix>=0.18
28 extra-libraries: sundials_arkode,
29 sundials_cvode
30 other-extensions: QuasiQuotes
31 hs-source-dirs: src
32 exposed-modules: Numeric.Sundials.ODEOpts,
33 Numeric.Sundials.ARKode.ODE,
34 Numeric.Sundials.CVode.ODE
35 other-modules: Numeric.Sundials.Arkode
36 c-sources: src/helpers.c src/helpers.h
37 default-language: Haskell2010
38
39test-suite hmatrix-sundials-testsuite
40 type: exitcode-stdio-1.0
41 main-is: Main.hs
42 other-modules: Numeric.Sundials.ODEOpts,
43 Numeric.Sundials.ARKode.ODE,
44 Numeric.Sundials.CVode.ODE,
45 Numeric.Sundials.Arkode
46 build-depends: base >=4.10 && <4.11,
47 inline-c >=0.6 && <0.7,
48 vector >=0.12 && <0.13,
49 template-haskell >=2.12 && <2.13,
50 containers >=0.5 && <0.6,
51 hmatrix>=0.18,
52 plots,
53 diagrams-lib,
54 diagrams-rasterific,
55 lens,
56 hspec
57 hs-source-dirs: src
58 extra-libraries: sundials_arkode,
59 sundials_cvode
60 c-sources: src/helpers.c src/helpers.h
61 default-language: Haskell2010
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
new file mode 100644
index 0000000..16c21c5
--- /dev/null
+++ b/packages/sundials/src/Main.hs
@@ -0,0 +1,186 @@
1{-# OPTIONS_GHC -Wall #-}
2
3import qualified Numeric.Sundials.ARKode.ODE as ARK
4import qualified Numeric.Sundials.CVode.ODE as CV
5import Numeric.LinearAlgebra
6
7import Plots as P
8import qualified Diagrams.Prelude as D
9import Diagrams.Backend.Rasterific
10
11import Control.Lens
12
13import Test.Hspec
14
15
16lorenz :: Double -> [Double] -> [Double]
17lorenz _t u = [ sigma * (y - x)
18 , x * (rho - z) - y
19 , x * y - beta * z
20 ]
21 where
22 rho = 28.0
23 sigma = 10.0
24 beta = 8.0 / 3.0
25 x = u !! 0
26 y = u !! 1
27 z = u !! 2
28
29_lorenzJac :: Double -> Vector Double -> Matrix Double
30_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y
31 , sigma , -1.0 , x
32 , 0.0 , (-x) , (-beta)
33 ]
34 where
35 rho = 28.0
36 sigma = 10.0
37 beta = 8.0 / 3.0
38 x = u ! 0
39 y = u ! 1
40 z = u ! 2
41
42brusselator :: Double -> [Double] -> [Double]
43brusselator _t x = [ a - (w + 1) * u + v * u * u
44 , w * u - v * u * u
45 , (b - w) / eps - w * u
46 ]
47 where
48 a = 1.0
49 b = 3.5
50 eps = 5.0e-6
51 u = x !! 0
52 v = x !! 1
53 w = x !! 2
54
55_brussJac :: Double -> Vector Double -> Matrix Double
56_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
57 , u * u , (-(u * u)) , 0.0
58 , (-u) , u , (-1.0) / eps - u
59 ]
60 where
61 y = toList x
62 u = y !! 0
63 v = y !! 1
64 w = y !! 2
65 eps = 5.0e-6
66
67stiffish :: Double -> [Double] -> [Double]
68stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
69 where
70 lamda = -100.0
71 u = v !! 0
72
73stiffishV :: Double -> Vector Double -> Vector Double
74stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
75 where
76 lamda = -100.0
77 u = v ! 0
78
79_stiffJac :: Double -> Vector Double -> Matrix Double
80_stiffJac _t _v = (1><1) [ lamda ]
81 where
82 lamda = -100.0
83
84predatorPrey :: Double -> [Double] -> [Double]
85predatorPrey _t v = [ x * a - b * x * y
86 , d * x * y - c * y - e * y * z
87 , (-f) * z + g * y * z
88 ]
89 where
90 x = v!!0
91 y = v!!1
92 z = v!!2
93 a = 1.0
94 b = 1.0
95 c = 1.0
96 d = 1.0
97 e = 1.0
98 f = 1.0
99 g = 1.0
100
101lSaxis :: [[Double]] -> P.Axis B D.V2 Double
102lSaxis xs = P.r2Axis &~ do
103 let ts = xs!!0
104 us = xs!!1
105 vs = xs!!2
106 ws = xs!!3
107 P.linePlot' $ zip ts us
108 P.linePlot' $ zip ts vs
109 P.linePlot' $ zip ts ws
110
111kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
112kSaxis xs = P.r2Axis &~ do
113 P.linePlot' xs
114
115main :: IO ()
116main = do
117
118 let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
119 renderRasterific "diagrams/brusselator.png"
120 (D.dims2D 500.0 500.0)
121 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
122
123 let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
124 renderRasterific "diagrams/brusselatorA.png"
125 (D.dims2D 500.0 500.0)
126 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))
127
128 let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
129 renderRasterific "diagrams/stiffish.png"
130 (D.dims2D 500.0 500.0)
131 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2))
132
133 let res2a = ARK.odeSolveV (ARK.SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
134
135 let res2b = ARK.odeSolveV (ARK.TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
136
137 let maxDiffA = maximum $ map abs $
138 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)
139
140 let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
141
142 let maxDiffB = maximum $ map abs $
143 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0)
144
145 let maxDiffC = maximum $ map abs $
146 zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0)
147
148 let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
149
150 renderRasterific "diagrams/lorenz.png"
151 (D.dims2D 500.0 500.0)
152 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))
153
154 renderRasterific "diagrams/lorenz1.png"
155 (D.dims2D 500.0 500.0)
156 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2))
157
158 renderRasterific "diagrams/lorenz2.png"
159 (D.dims2D 500.0 500.0)
160 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2))
161
162 let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
163
164 renderRasterific "diagrams/predatorPrey.png"
165 (D.dims2D 500.0 500.0)
166 (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1))
167
168 renderRasterific "diagrams/predatorPrey1.png"
169 (D.dims2D 500.0 500.0)
170 (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2))
171
172 renderRasterific "diagrams/predatorPrey2.png"
173 (D.dims2D 500.0 500.0)
174 (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2))
175
176 let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
177
178 let maxDiffPpA = maximum $ map abs $
179 zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0)
180
181 hspec $ describe "Compare results" $ do
182 it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6
183 it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6
184 it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6
185 it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3
186
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
new file mode 100644
index 0000000..fafc237
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -0,0 +1,903 @@
1{-# LANGUAGE QuasiQuotes #-}
2{-# LANGUAGE TemplateHaskell #-}
3{-# LANGUAGE MultiWayIf #-}
4{-# LANGUAGE OverloadedStrings #-}
5{-# LANGUAGE ScopedTypeVariables #-}
6{-# LANGUAGE DeriveGeneric #-}
7{-# LANGUAGE TypeOperators #-}
8{-# LANGUAGE KindSignatures #-}
9{-# LANGUAGE TypeSynonymInstances #-}
10{-# LANGUAGE FlexibleInstances #-}
11{-# LANGUAGE FlexibleContexts #-}
12
13-----------------------------------------------------------------------------
14-- |
15-- Module : Numeric.Sundials.ARKode.ODE
16-- Copyright : Dominic Steinitz 2018,
17-- Novadiscovery 2018
18-- License : BSD
19-- Maintainer : Dominic Steinitz
20-- Stability : provisional
21--
22-- Solution of ordinary differential equation (ODE) initial value problems.
23-- See <https://computation.llnl.gov/projects/sundials/sundials-software> for more detail.
24--
25-- A simple example:
26--
27-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
28--
29-- @
30-- import Numeric.Sundials.ARKode.ODE
31-- import Numeric.LinearAlgebra
32--
33-- import Plots as P
34-- import qualified Diagrams.Prelude as D
35-- import Diagrams.Backend.Rasterific
36--
37-- brusselator :: Double -> [Double] -> [Double]
38-- brusselator _t x = [ a - (w + 1) * u + v * u * u
39-- , w * u - v * u * u
40-- , (b - w) / eps - w * u
41-- ]
42-- where
43-- a = 1.0
44-- b = 3.5
45-- eps = 5.0e-6
46-- u = x !! 0
47-- v = x !! 1
48-- w = x !! 2
49--
50-- lSaxis :: [[Double]] -> P.Axis B D.V2 Double
51-- lSaxis xs = P.r2Axis &~ do
52-- let ts = xs!!0
53-- us = xs!!1
54-- vs = xs!!2
55-- ws = xs!!3
56-- P.linePlot' $ zip ts us
57-- P.linePlot' $ zip ts vs
58-- P.linePlot' $ zip ts ws
59--
60-- main = do
61-- let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
62-- renderRasterific "diagrams/brusselator.png"
63-- (D.dims2D 500.0 500.0)
64-- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
65-- @
66--
67-- With Sundials ARKode, it is possible to retrieve the Butcher tableau for the solver.
68--
69-- @
70-- import Numeric.Sundials.ARKode.ODE
71-- import Numeric.LinearAlgebra
72--
73-- import Data.List (intercalate)
74--
75-- import Text.PrettyPrint.HughesPJClass
76--
77--
78-- butcherTableauTex :: ButcherTable -> String
79-- butcherTableauTex (ButcherTable m c b b2) =
80-- render $
81-- vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}")
82-- , us
83-- , text "\\hline"
84-- , text bs <+> text "\\\\"
85-- , text b2s <+> text "\\\\"
86-- , text "\\end{array}"
87-- ]
88-- where
89-- n = rows m
90-- rs = toLists m
91-- ss = map (\r -> intercalate " & " $ map show r) rs
92-- ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss
93-- us = vcat $ map (\r -> text r <+> text "\\\\") ts
94-- bs = " & " ++ (intercalate " & " $ map show $ toList b)
95-- b2s = " & " ++ (intercalate " & " $ map show $ toList b2)
96--
97-- main :: IO ()
98-- main = do
99--
100-- let res = butcherTable (SDIRK_2_1_2 undefined)
101-- putStrLn $ show res
102-- putStrLn $ butcherTableauTex res
103--
104-- let resA = butcherTable (KVAERNO_4_2_3 undefined)
105-- putStrLn $ show resA
106-- putStrLn $ butcherTableauTex resA
107--
108-- let resB = butcherTable (SDIRK_5_3_4 undefined)
109-- putStrLn $ show resB
110-- putStrLn $ butcherTableauTex resB
111-- @
112--
113-- Using the code above from the examples gives
114--
115-- KVAERNO_4_2_3
116--
117-- \[
118-- \begin{array}{c|cccc}
119-- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
120-- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\
121-- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
122-- 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
123-- \hline
124-- & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
125-- & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
126-- \end{array}
127-- \]
128--
129-- SDIRK_2_1_2
130--
131-- \[
132-- \begin{array}{c|cc}
133-- 1.0 & 1.0 & 0.0 \\
134-- 0.0 & -1.0 & 1.0 \\
135-- \hline
136-- & 0.5 & 0.5 \\
137-- & 1.0 & 0.0 \\
138-- \end{array}
139-- \]
140--
141-- SDIRK_5_3_4
142--
143-- \[
144-- \begin{array}{c|ccccc}
145-- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\
146-- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\
147-- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\
148-- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\
149-- 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
150-- \hline
151-- & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
152-- & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\
153-- \end{array}
154-- \]
155-----------------------------------------------------------------------------
156module Numeric.Sundials.ARKode.ODE ( odeSolve
157 , odeSolveV
158 , odeSolveVWith
159 , odeSolveVWith'
160 , ButcherTable(..)
161 , butcherTable
162 , ODEMethod(..)
163 , StepControl(..)
164 ) where
165
166import qualified Language.C.Inline as C
167import qualified Language.C.Inline.Unsafe as CU
168
169import Data.Monoid ((<>))
170import Data.Maybe (isJust)
171
172import Foreign.C.Types (CDouble, CInt, CLong)
173import Foreign.Ptr (Ptr)
174import Foreign.Storable (poke)
175
176import qualified Data.Vector.Storable as V
177
178import Data.Coerce (coerce)
179import System.IO.Unsafe (unsafePerformIO)
180import GHC.Generics (C1, Constructor, (:+:)(..), D1, Rep, Generic, M1(..),
181 from, conName)
182
183import Numeric.LinearAlgebra.Devel (createVector)
184
185import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
186 cols, toLists, size, reshape,
187 subVector, subMatrix, (><))
188
189import Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..))
190import qualified Numeric.Sundials.Arkode as T
191import Numeric.Sundials.Arkode (getDataFromContents, putDataInContents, arkSMax,
192 sDIRK_2_1_2,
193 bILLINGTON_3_3_2,
194 tRBDF2_3_3_2,
195 kVAERNO_4_2_3,
196 aRK324L2SA_DIRK_4_2_3,
197 cASH_5_2_4,
198 cASH_5_3_4,
199 sDIRK_5_3_4,
200 kVAERNO_5_3_4,
201 aRK436L2SA_DIRK_6_3_4,
202 kVAERNO_7_4_5,
203 aRK548L2SA_DIRK_8_4_5,
204 hEUN_EULER_2_1_2,
205 bOGACKI_SHAMPINE_4_2_3,
206 aRK324L2SA_ERK_4_2_3,
207 zONNEVELD_5_3_4,
208 aRK436L2SA_ERK_6_3_4,
209 sAYFY_ABURUB_6_3_4,
210 cASH_KARP_6_4_5,
211 fEHLBERG_6_4_5,
212 dORMAND_PRINCE_7_4_5,
213 aRK548L2SA_ERK_8_4_5,
214 vERNER_8_5_6,
215 fEHLBERG_13_7_8)
216
217
218C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
219
220C.include "<stdlib.h>"
221C.include "<stdio.h>"
222C.include "<math.h>"
223C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts.
224C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros
225C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
226C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
227C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface
228C.include "<sundials/sundials_types.h>" -- definition of type realtype
229C.include "<sundials/sundials_math.h>"
230C.include "../../../helpers.h"
231C.include "Numeric/Sundials/Arkode_hsc.h"
232
233
234-- | Stepping functions
235data ODEMethod = SDIRK_2_1_2 Jacobian
236 | SDIRK_2_1_2'
237 | BILLINGTON_3_3_2 Jacobian
238 | BILLINGTON_3_3_2'
239 | TRBDF2_3_3_2 Jacobian
240 | TRBDF2_3_3_2'
241 | KVAERNO_4_2_3 Jacobian
242 | KVAERNO_4_2_3'
243 | ARK324L2SA_DIRK_4_2_3 Jacobian
244 | ARK324L2SA_DIRK_4_2_3'
245 | CASH_5_2_4 Jacobian
246 | CASH_5_2_4'
247 | CASH_5_3_4 Jacobian
248 | CASH_5_3_4'
249 | SDIRK_5_3_4 Jacobian
250 | SDIRK_5_3_4'
251 | KVAERNO_5_3_4 Jacobian
252 | KVAERNO_5_3_4'
253 | ARK436L2SA_DIRK_6_3_4 Jacobian
254 | ARK436L2SA_DIRK_6_3_4'
255 | KVAERNO_7_4_5 Jacobian
256 | KVAERNO_7_4_5'
257 | ARK548L2SA_DIRK_8_4_5 Jacobian
258 | ARK548L2SA_DIRK_8_4_5'
259 | HEUN_EULER_2_1_2 Jacobian
260 | HEUN_EULER_2_1_2'
261 | BOGACKI_SHAMPINE_4_2_3 Jacobian
262 | BOGACKI_SHAMPINE_4_2_3'
263 | ARK324L2SA_ERK_4_2_3 Jacobian
264 | ARK324L2SA_ERK_4_2_3'
265 | ZONNEVELD_5_3_4 Jacobian
266 | ZONNEVELD_5_3_4'
267 | ARK436L2SA_ERK_6_3_4 Jacobian
268 | ARK436L2SA_ERK_6_3_4'
269 | SAYFY_ABURUB_6_3_4 Jacobian
270 | SAYFY_ABURUB_6_3_4'
271 | CASH_KARP_6_4_5 Jacobian
272 | CASH_KARP_6_4_5'
273 | FEHLBERG_6_4_5 Jacobian
274 | FEHLBERG_6_4_5'
275 | DORMAND_PRINCE_7_4_5 Jacobian
276 | DORMAND_PRINCE_7_4_5'
277 | ARK548L2SA_ERK_8_4_5 Jacobian
278 | ARK548L2SA_ERK_8_4_5'
279 | VERNER_8_5_6 Jacobian
280 | VERNER_8_5_6'
281 | FEHLBERG_13_7_8 Jacobian
282 | FEHLBERG_13_7_8'
283 deriving Generic
284
285constrName :: (HasConstructor (Rep a), Generic a)=> a -> String
286constrName = genericConstrName . from
287
288class HasConstructor (f :: * -> *) where
289 genericConstrName :: f x -> String
290
291instance HasConstructor f => HasConstructor (D1 c f) where
292 genericConstrName (M1 x) = genericConstrName x
293
294instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where
295 genericConstrName (L1 l) = genericConstrName l
296 genericConstrName (R1 r) = genericConstrName r
297
298instance Constructor c => HasConstructor (C1 c f) where
299 genericConstrName x = conName x
300
301instance Show ODEMethod where
302 show x = constrName x
303
304-- FIXME: We can probably do better here with generics
305getMethod :: ODEMethod -> Int
306getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2
307getMethod (SDIRK_2_1_2') = sDIRK_2_1_2
308getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2
309getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2
310getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2
311getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2
312getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3
313getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3
314getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3
315getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3
316getMethod (CASH_5_2_4 _) = cASH_5_2_4
317getMethod (CASH_5_2_4') = cASH_5_2_4
318getMethod (CASH_5_3_4 _) = cASH_5_3_4
319getMethod (CASH_5_3_4') = cASH_5_3_4
320getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4
321getMethod (SDIRK_5_3_4') = sDIRK_5_3_4
322getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4
323getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4
324getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4
325getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4
326getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5
327getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5
328getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5
329getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5
330getMethod (HEUN_EULER_2_1_2 _) = hEUN_EULER_2_1_2
331getMethod (HEUN_EULER_2_1_2') = hEUN_EULER_2_1_2
332getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3
333getMethod (BOGACKI_SHAMPINE_4_2_3') = bOGACKI_SHAMPINE_4_2_3
334getMethod (ARK324L2SA_ERK_4_2_3 _) = aRK324L2SA_ERK_4_2_3
335getMethod (ARK324L2SA_ERK_4_2_3') = aRK324L2SA_ERK_4_2_3
336getMethod (ZONNEVELD_5_3_4 _) = zONNEVELD_5_3_4
337getMethod (ZONNEVELD_5_3_4') = zONNEVELD_5_3_4
338getMethod (ARK436L2SA_ERK_6_3_4 _) = aRK436L2SA_ERK_6_3_4
339getMethod (ARK436L2SA_ERK_6_3_4') = aRK436L2SA_ERK_6_3_4
340getMethod (SAYFY_ABURUB_6_3_4 _) = sAYFY_ABURUB_6_3_4
341getMethod (SAYFY_ABURUB_6_3_4') = sAYFY_ABURUB_6_3_4
342getMethod (CASH_KARP_6_4_5 _) = cASH_KARP_6_4_5
343getMethod (CASH_KARP_6_4_5') = cASH_KARP_6_4_5
344getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5
345getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5
346getMethod (DORMAND_PRINCE_7_4_5 _) = dORMAND_PRINCE_7_4_5
347getMethod (DORMAND_PRINCE_7_4_5') = dORMAND_PRINCE_7_4_5
348getMethod (ARK548L2SA_ERK_8_4_5 _) = aRK548L2SA_ERK_8_4_5
349getMethod (ARK548L2SA_ERK_8_4_5') = aRK548L2SA_ERK_8_4_5
350getMethod (VERNER_8_5_6 _) = vERNER_8_5_6
351getMethod (VERNER_8_5_6') = vERNER_8_5_6
352getMethod (FEHLBERG_13_7_8 _) = fEHLBERG_13_7_8
353getMethod (FEHLBERG_13_7_8') = fEHLBERG_13_7_8
354
355getJacobian :: ODEMethod -> Maybe Jacobian
356getJacobian (SDIRK_2_1_2 j) = Just j
357getJacobian (BILLINGTON_3_3_2 j) = Just j
358getJacobian (TRBDF2_3_3_2 j) = Just j
359getJacobian (KVAERNO_4_2_3 j) = Just j
360getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j
361getJacobian (CASH_5_2_4 j) = Just j
362getJacobian (CASH_5_3_4 j) = Just j
363getJacobian (SDIRK_5_3_4 j) = Just j
364getJacobian (KVAERNO_5_3_4 j) = Just j
365getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j
366getJacobian (KVAERNO_7_4_5 j) = Just j
367getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j
368getJacobian (HEUN_EULER_2_1_2 j) = Just j
369getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j
370getJacobian (ARK324L2SA_ERK_4_2_3 j) = Just j
371getJacobian (ZONNEVELD_5_3_4 j) = Just j
372getJacobian (ARK436L2SA_ERK_6_3_4 j) = Just j
373getJacobian (SAYFY_ABURUB_6_3_4 j) = Just j
374getJacobian (CASH_KARP_6_4_5 j) = Just j
375getJacobian (FEHLBERG_6_4_5 j) = Just j
376getJacobian (DORMAND_PRINCE_7_4_5 j) = Just j
377getJacobian (ARK548L2SA_ERK_8_4_5 j) = Just j
378getJacobian (VERNER_8_5_6 j) = Just j
379getJacobian (FEHLBERG_13_7_8 j) = Just j
380getJacobian _ = Nothing
381
382-- | A version of 'odeSolveVWith' with reasonable default step control.
383odeSolveV
384 :: ODEMethod
385 -> Maybe Double -- ^ initial step size - by default, ARKode
386 -- estimates the initial step size to be the
387 -- solution \(h\) of the equation
388 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
389 -- \(\ddot{y}\) is an estimated value of the
390 -- second derivative of the solution at \(t_0\)
391 -> Double -- ^ absolute tolerance for the state vector
392 -> Double -- ^ relative tolerance for the state vector
393 -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
394 -> Vector Double -- ^ initial conditions
395 -> Vector Double -- ^ desired solution times
396 -> Matrix Double -- ^ solution
397odeSolveV meth hi epsAbs epsRel f y0 ts =
398 odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts
399 where
400 g t x0 = coerce $ f t x0
401
402-- | A version of 'odeSolveV' with reasonable default parameters and
403-- system of equations defined using lists. FIXME: we should say
404-- something about the fact we could use the Jacobian but don't for
405-- compatibility with hmatrix-gsl.
406odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
407 -> [Double] -- ^ initial conditions
408 -> Vector Double -- ^ desired solution times
409 -> Matrix Double -- ^ solution
410odeSolve f y0 ts =
411 -- FIXME: These tolerances are different from the ones in GSL
412 odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts)
413 where
414 g t x0 = V.fromList $ f t (V.toList x0)
415
416odeSolveVWith ::
417 ODEMethod
418 -> StepControl
419 -> Maybe Double -- ^ initial step size - by default, ARKode
420 -- estimates the initial step size to be the
421 -- solution \(h\) of the equation
422 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
423 -- \(\ddot{y}\) is an estimated value of the second
424 -- derivative of the solution at \(t_0\)
425 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
426 -> V.Vector Double -- ^ Initial conditions
427 -> V.Vector Double -- ^ Desired solution times
428 -> Matrix Double -- ^ Error code or solution
429odeSolveVWith method control initStepSize f y0 tt =
430 case odeSolveVWith' opts method control initStepSize f y0 tt of
431 Left c -> error $ show c -- FIXME
432 Right (v, _d) -> v
433 where
434 opts = ODEOpts { maxNumSteps = 10000
435 , minStep = 1.0e-12
436 , relTol = error "relTol"
437 , absTols = error "absTol"
438 , initStep = error "initStep"
439 , maxFail = 10
440 }
441
442odeSolveVWith' ::
443 ODEOpts
444 -> ODEMethod
445 -> StepControl
446 -> Maybe Double -- ^ initial step size - by default, ARKode
447 -- estimates the initial step size to be the
448 -- solution \(h\) of the equation
449 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
450 -- \(\ddot{y}\) is an estimated value of the second
451 -- derivative of the solution at \(t_0\)
452 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
453 -> V.Vector Double -- ^ Initial conditions
454 -> V.Vector Double -- ^ Desired solution times
455 -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
456odeSolveVWith' opts method control initStepSize f y0 tt =
457 case solveOdeC (fromIntegral $ maxFail opts)
458 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
459 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
460 (coerce f) (coerce y0) (coerce tt) of
461 Left c -> Left $ fromIntegral c
462 Right (v, d) -> Right (reshape l (coerce v), d)
463 where
464 l = size y0
465 scise (X aTol rTol) = coerce (V.replicate l aTol, rTol)
466 scise (X' aTol rTol) = coerce (V.replicate l aTol, rTol)
467 scise (XX' aTol rTol yScale _yDotScale) = coerce (V.replicate l aTol, yScale * rTol)
468 -- FIXME; Should we check that the length of ss is correct?
469 scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
470 jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
471 getJacobian method
472 matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
473 where
474 nr = fromIntegral $ rows m
475 nc = fromIntegral $ cols m
476 -- FIXME: efficiency
477 vs = V.fromList $ map coerce $ concat $ toLists m
478
479solveOdeC ::
480 CInt ->
481 CLong ->
482 CDouble ->
483 CInt ->
484 Maybe CDouble ->
485 (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
486 (V.Vector CDouble, CDouble) ->
487 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
488 -> V.Vector CDouble -- ^ Initial conditions
489 -> V.Vector CDouble -- ^ Desired solution times
490 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
491solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
492 jacH (aTols, rTol) fun f0 ts = unsafePerformIO $ do
493
494 let isInitStepSize :: CInt
495 isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
496 ss :: CDouble
497 ss = case initStepSize of
498 -- It would be better to put an error message here but
499 -- inline-c seems to evaluate this even if it is never
500 -- used :(
501 Nothing -> 0.0
502 Just x -> x
503
504 let dim = V.length f0
505 nEq :: CLong
506 nEq = fromIntegral dim
507 nTs :: CInt
508 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))
512 qMatMut <- V.thaw quasiMatrixRes
513 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
514 diagMut <- V.thaw diagnostics
515 -- We need the types that sundials expects. These are tied together
516 -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
517 let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
518 funIO x y f _ptr = do
519 -- Convert the pointer we get from C (y) to a vector, and then
520 -- apply the user-supplied function.
521 fImm <- fun x <$> getDataFromContents dim y
522 -- Fill in the provided pointer with the resulting vector.
523 putDataInContents fImm dim f
524 -- FIXME: I don't understand what this comment means
525 -- Unsafe since the function will be called many times.
526 [CU.exp| int{ 0 } |]
527 let isJac :: CInt
528 isJac = fromIntegral $ fromEnum $ isJust jacH
529 jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
530 Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
531 IO CInt
532 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
533 case jacH of
534 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
535 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
536 poke jacS j
537 -- FIXME: I don't understand what this comment means
538 -- Unsafe since the function will be called many times.
539 [CU.exp| int{ 0 } |]
540
541 res <- [C.block| int {
542 /* general problem variables */
543
544 int flag; /* reusable error-checking flag */
545 int i, j; /* reusable loop indices */
546 N_Vector y = NULL; /* empty vector for storing solution */
547 N_Vector tv = NULL; /* empty vector for storing absolute tolerances */
548 SUNMatrix A = NULL; /* empty matrix for linear solver */
549 SUNLinearSolver LS = NULL; /* empty linear solver object */
550 void *arkode_mem = NULL; /* empty ARKode memory structure */
551 realtype t;
552 long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
553
554 /* general problem parameters */
555
556 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
557 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */
558
559 /* Initialize data structures */
560
561 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
562 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
563 /* Specify initial condition */
564 for (i = 0; i < NEQ; i++) {
565 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
566 };
567
568 tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
569 if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
570 /* Specify tolerances */
571 for (i = 0; i < NEQ; i++) {
572 NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
573 };
574
575 arkode_mem = ARKodeCreate(); /* Create the solver memory */
576 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
577
578 /* Call ARKodeInit to initialize the integrator memory and specify the */
579 /* right-hand side function in y'=f(t,y), the inital time T0, and */
580 /* the initial dependent variable vector y. Note: we treat the */
581 /* problem as fully implicit and set f_E to NULL and f_I to f. */
582
583 /* Here we use the C types defined in helpers.h which tie up with */
584 /* the Haskell types defined in CLangToHaskellTypes */
585 if ($(int method) < MIN_DIRK_NUM) {
586 flag = ARKodeInit(arkode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y);
587 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
588 } else {
589 flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
590 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
591 }
592
593 flag = ARKodeSetMinStep(arkode_mem, $(double minStep_));
594 if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1;
595 flag = ARKodeSetMaxNumSteps(arkode_mem, $(long int maxNumSteps_));
596 if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1;
597 flag = ARKodeSetMaxErrTestFails(arkode_mem, $(int maxErrTestFails));
598 if (check_flag(&flag, "ARKodeSetMaxErrTestFails", 1)) return 1;
599
600 /* Set routines */
601 flag = ARKodeSVtolerances(arkode_mem, $(double rTol), tv);
602 if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1;
603
604 /* Initialize dense matrix data structure and solver */
605 A = SUNDenseMatrix(NEQ, NEQ);
606 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
607 LS = SUNDenseLinearSolver(y, A);
608 if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
609
610 /* Attach matrix and linear solver */
611 flag = ARKDlsSetLinearSolver(arkode_mem, LS, A);
612 if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
613
614 /* Set the initial step size if there is one */
615 if ($(int isInitStepSize)) {
616 /* FIXME: We could check if the initial step size is 0 */
617 /* or even NaN and then throw an error */
618 flag = ARKodeSetInitStep(arkode_mem, $(double ss));
619 if (check_flag(&flag, "ARKodeSetInitStep", 1)) return 1;
620 }
621
622 /* Set the Jacobian if there is one */
623 if ($(int isJac)) {
624 flag = ARKDlsSetJacFn(arkode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
625 if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
626 }
627
628 /* Store initial conditions */
629 for (j = 0; j < NEQ; j++) {
630 ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
631 }
632
633 /* Explicitly set the method */
634 if ($(int method) >= MIN_DIRK_NUM) {
635 flag = ARKodeSetIRKTableNum(arkode_mem, $(int method));
636 if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
637 } else {
638 flag = ARKodeSetERKTableNum(arkode_mem, $(int method));
639 if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
640 }
641
642 /* Main time-stepping loop: calls ARKode to perform the integration */
643 /* Stops when the final time has been reached */
644 for (i = 1; i < $(int nTs); i++) {
645
646 flag = ARKode(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */
647 if (check_flag(&flag, "ARKode", 1)) break;
648
649 /* Store the results for Haskell */
650 for (j = 0; j < NEQ; j++) {
651 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
652 }
653
654 /* unsuccessful solve: break */
655 if (flag < 0) {
656 fprintf(stderr,"Solver failure, stopping integration\n");
657 break;
658 }
659 }
660
661 /* Get some final statistics on how the solve progressed */
662
663 flag = ARKodeGetNumSteps(arkode_mem, &nst);
664 check_flag(&flag, "ARKodeGetNumSteps", 1);
665 ($vec-ptr:(long int *diagMut))[0] = nst;
666
667 flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
668 check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
669 ($vec-ptr:(long int *diagMut))[1] = nst_a;
670
671 flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
672 check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
673 ($vec-ptr:(long int *diagMut))[2] = nfe;
674 ($vec-ptr:(long int *diagMut))[3] = nfi;
675
676 flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
677 check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
678 ($vec-ptr:(long int *diagMut))[4] = nsetups;
679
680 flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
681 check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
682 ($vec-ptr:(long int *diagMut))[5] = netf;
683
684 flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
685 check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
686 ($vec-ptr:(long int *diagMut))[6] = nni;
687
688 flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
689 check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
690 ($vec-ptr:(long int *diagMut))[7] = ncfn;
691
692 flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
693 check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
694 ($vec-ptr:(long int *diagMut))[8] = ncfn;
695
696 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
697 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
698 ($vec-ptr:(long int *diagMut))[9] = ncfn;
699
700 /* Clean up and return */
701 N_VDestroy(y); /* Free y vector */
702 N_VDestroy(tv); /* Free tv vector */
703 ARKodeFree(&arkode_mem); /* Free integrator memory */
704 SUNLinSolFree(LS); /* Free linear solver */
705 SUNMatDestroy(A); /* Free A matrix */
706
707 return flag;
708 } |]
709 if res == 0
710 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)
724 else do
725 return $ Left res
726
727data ButcherTable = ButcherTable { am :: Matrix Double
728 , cv :: Vector Double
729 , bv :: Vector Double
730 , b2v :: Vector Double
731 }
732 deriving Show
733
734data ButcherTable' a = ButcherTable' { am' :: V.Vector a
735 , cv' :: V.Vector a
736 , bv' :: V.Vector a
737 , b2v' :: V.Vector a
738 }
739 deriving Show
740
741butcherTable :: ODEMethod -> ButcherTable
742butcherTable method =
743 case getBT method of
744 Left c -> error $ show c -- FIXME
745 Right (ButcherTable' v w x y, sqp) ->
746 ButcherTable { am = subMatrix (0, 0) (s, s) $ (arkSMax >< arkSMax) (V.toList v)
747 , cv = subVector 0 s w
748 , bv = subVector 0 s x
749 , b2v = subVector 0 s y
750 }
751 where
752 s = fromIntegral $ sqp V.! 0
753
754getBT :: ODEMethod -> Either Int (ButcherTable' Double, V.Vector Int)
755getBT method = case getButcherTable method of
756 Left c ->
757 Left $ fromIntegral c
758 Right (ButcherTable' a b c d, sqp) ->
759 Right $ ( ButcherTable' (coerce a) (coerce b) (coerce c) (coerce d)
760 , V.map fromIntegral sqp )
761
762getButcherTable :: ODEMethod
763 -> Either CInt (ButcherTable' CDouble, V.Vector CInt)
764getButcherTable method = unsafePerformIO $ do
765 -- ARKode seems to want an ODE in order to set and then get the
766 -- Butcher tableau so here's one to keep it happy
767 let funI :: CDouble -> V.Vector CDouble -> V.Vector CDouble
768 funI _t ys = V.fromList [ ys V.! 0 ]
769 let funE :: CDouble -> V.Vector CDouble -> V.Vector CDouble
770 funE _t ys = V.fromList [ ys V.! 0 ]
771 f0 = V.fromList [ 1.0 ]
772 ts = V.fromList [ 0.0 ]
773 dim = V.length f0
774 nEq :: CLong
775 nEq = fromIntegral dim
776 mN :: CInt
777 mN = fromIntegral $ getMethod method
778
779 btSQP :: V.Vector CInt <- createVector 3
780 btSQPMut <- V.thaw btSQP
781 btAs :: V.Vector CDouble <- createVector (arkSMax * arkSMax)
782 btAsMut <- V.thaw btAs
783 btCs :: V.Vector CDouble <- createVector arkSMax
784 btBs :: V.Vector CDouble <- createVector arkSMax
785 btB2s :: V.Vector CDouble <- createVector arkSMax
786 btCsMut <- V.thaw btCs
787 btBsMut <- V.thaw btBs
788 btB2sMut <- V.thaw btB2s
789 let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
790 funIOI x y f _ptr = do
791 fImm <- funI x <$> getDataFromContents dim y
792 putDataInContents fImm dim f
793 -- FIXME: I don't understand what this comment means
794 -- Unsafe since the function will be called many times.
795 [CU.exp| int{ 0 } |]
796 let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
797 funIOE x y f _ptr = do
798 fImm <- funE x <$> getDataFromContents dim y
799 putDataInContents fImm dim f
800 -- FIXME: I don't understand what this comment means
801 -- Unsafe since the function will be called many times.
802 [CU.exp| int{ 0 } |]
803 res <- [C.block| int {
804 /* general problem variables */
805
806 int flag; /* reusable error-checking flag */
807 N_Vector y = NULL; /* empty vector for storing solution */
808 void *arkode_mem = NULL; /* empty ARKode memory structure */
809 int i, j; /* reusable loop indices */
810
811 /* general problem parameters */
812
813 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
814 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars */
815
816 /* Initialize data structures */
817
818 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
819 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
820 /* Specify initial condition */
821 for (i = 0; i < NEQ; i++) {
822 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
823 };
824 arkode_mem = ARKodeCreate(); /* Create the solver memory */
825 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
826
827 flag = ARKodeInit(arkode_mem, $fun:(int (* funIOE) (double t, SunVector y[], SunVector dydt[], void * params)), $fun:(int (* funIOI) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
828 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
829
830 if ($(int mN) >= MIN_DIRK_NUM) {
831 flag = ARKodeSetIRKTableNum(arkode_mem, $(int mN));
832 if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
833 } else {
834 flag = ARKodeSetERKTableNum(arkode_mem, $(int mN));
835 if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
836 }
837
838 int s, q, p;
839 realtype *ai = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
840 realtype *ae = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
841 realtype *ci = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
842 realtype *ce = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
843 realtype *bi = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
844 realtype *be = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
845 realtype *b2i = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
846 realtype *b2e = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
847 flag = ARKodeGetCurrentButcherTables(arkode_mem, &s, &q, &p, ai, ae, ci, ce, bi, be, b2i, b2e);
848 if (check_flag(&flag, "ARKode", 1)) return 1;
849 $vec-ptr:(int *btSQPMut)[0] = s;
850 $vec-ptr:(int *btSQPMut)[1] = q;
851 $vec-ptr:(int *btSQPMut)[2] = p;
852 for (i = 0; i < s; i++) {
853 for (j = 0; j < s; j++) {
854 /* FIXME: double should be realtype */
855 ($vec-ptr:(double *btAsMut))[i * ARK_S_MAX + j] = ai[i * ARK_S_MAX + j];
856 }
857 }
858
859 for (i = 0; i < s; i++) {
860 ($vec-ptr:(double *btCsMut))[i] = ci[i];
861 ($vec-ptr:(double *btBsMut))[i] = bi[i];
862 ($vec-ptr:(double *btB2sMut))[i] = b2i[i];
863 }
864
865 /* Clean up and return */
866 N_VDestroy(y); /* Free y vector */
867 ARKodeFree(&arkode_mem); /* Free integrator memory */
868
869 return flag;
870 } |]
871 if res == 0
872 then do
873 x <- V.freeze btAsMut
874 y <- V.freeze btSQPMut
875 z <- V.freeze btCsMut
876 u <- V.freeze btBsMut
877 v <- V.freeze btB2sMut
878 return $ Right (ButcherTable' { am' = x, cv' = z, bv' = u, b2v' = v }, y)
879 else do
880 return $ Left res
881
882-- | Adaptive step-size control
883-- functions.
884--
885-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
886-- allows the user to control the step size adjustment using
887-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
888-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
889-- is the required relative error, \(s_i\) is a vector of scaling
890-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
891-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
892--
893-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
894-- allows the user to control the step size adjustment using
895-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
896-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
897-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
898-- effect.
899data StepControl = X Double Double -- ^ absolute and relative tolerance for \(y\); in GSL terms, \(a_{y} = 1\) and \(a_{dy/dt} = 0\); in ARKode terms, the \(\eta^{abs}_i\) are identical
900 | X' Double Double -- ^ absolute and relative tolerance for \(\dot{y}\); in GSL terms, \(a_{y} = 0\) and \(a_{dy/dt} = 1\); in ARKode terms, the latter is treated as the relative tolerance for \(y\) so this is the same as specifying 'X' which may be entirely incorrect for the given problem
901 | XX' Double Double Double Double -- ^ include both via relative tolerance
902 -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
903 | ScXX' Double Double Double Double (Vector Double) -- ^ scale absolute tolerance of \(y_i\); in ARKode terms, \(a_{{dy}/{dt}}\) is ignored, \(\eta^{abs}_i = s_i \epsilon^{abs}\) and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
diff --git a/packages/sundials/src/Numeric/Sundials/Arkode.hsc b/packages/sundials/src/Numeric/Sundials/Arkode.hsc
new file mode 100644
index 0000000..0850258
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/Arkode.hsc
@@ -0,0 +1,204 @@
1{-# LANGUAGE QuasiQuotes #-}
2{-# LANGUAGE TemplateHaskell #-}
3{-# LANGUAGE OverloadedStrings #-}
4{-# LANGUAGE EmptyDataDecls #-}
5
6module Numeric.Sundials.Arkode where
7
8import Foreign
9import Foreign.C.Types
10
11import Language.C.Types as CT
12
13import qualified Data.Vector.Storable as VS
14import qualified Data.Vector.Storable.Mutable as VM
15
16import qualified Language.Haskell.TH as TH
17import qualified Data.Map as Map
18import Language.C.Inline.Context
19
20import qualified Data.Vector.Storable as V
21
22
23#include <stdio.h>
24#include <sundials/sundials_nvector.h>
25#include <sundials/sundials_matrix.h>
26#include <nvector/nvector_serial.h>
27#include <sunmatrix/sunmatrix_dense.h>
28#include <arkode/arkode.h>
29#include <cvode/cvode.h>
30
31
32data SunVector
33data SunMatrix = SunMatrix { rows :: CInt
34 , cols :: CInt
35 , vals :: V.Vector CDouble
36 }
37
38-- | This is true only if configured/ built as 64 bits
39type SunIndexType = CLong
40
41sunTypesTable :: Map.Map TypeSpecifier TH.TypeQ
42sunTypesTable = Map.fromList
43 [
44 (TypeName "sunindextype", [t| SunIndexType |] )
45 , (TypeName "SunVector", [t| SunVector |] )
46 , (TypeName "SunMatrix", [t| SunMatrix |] )
47 ]
48
49sunCtx :: Context
50sunCtx = mempty {ctxTypesTable = sunTypesTable}
51
52getMatrixDataFromContents :: Ptr SunMatrix -> IO SunMatrix
53getMatrixDataFromContents ptr = do
54 qtr <- getContentMatrixPtr ptr
55 rs <- getNRows qtr
56 cs <- getNCols qtr
57 rtr <- getMatrixData qtr
58 vs <- vectorFromC (fromIntegral $ rs * cs) rtr
59 return $ SunMatrix { rows = rs, cols = cs, vals = vs }
60
61putMatrixDataFromContents :: SunMatrix -> Ptr SunMatrix -> IO ()
62putMatrixDataFromContents mat ptr = do
63 let rs = rows mat
64 cs = cols mat
65 vs = vals mat
66 qtr <- getContentMatrixPtr ptr
67 putNRows rs qtr
68 putNCols cs qtr
69 rtr <- getMatrixData qtr
70 vectorToC vs (fromIntegral $ rs * cs) rtr
71
72instance Storable SunMatrix where
73 poke = flip putMatrixDataFromContents
74 peek = getMatrixDataFromContents
75 sizeOf _ = error "sizeOf not supported for SunMatrix"
76 alignment _ = error "alignment not supported for SunMatrix"
77
78vectorFromC :: Storable a => Int -> Ptr a -> IO (VS.Vector a)
79vectorFromC len ptr = do
80 ptr' <- newForeignPtr_ ptr
81 VS.freeze $ VM.unsafeFromForeignPtr0 ptr' len
82
83vectorToC :: Storable a => VS.Vector a -> Int -> Ptr a -> IO ()
84vectorToC vec len ptr = do
85 ptr' <- newForeignPtr_ ptr
86 VS.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
87
88getDataFromContents :: Int -> Ptr SunVector -> IO (VS.Vector CDouble)
89getDataFromContents len ptr = do
90 qtr <- getContentPtr ptr
91 rtr <- getData qtr
92 vectorFromC len rtr
93
94putDataInContents :: Storable a => VS.Vector a -> Int -> Ptr b -> IO ()
95putDataInContents vec len ptr = do
96 qtr <- getContentPtr ptr
97 rtr <- getData qtr
98 vectorToC vec len rtr
99
100#def typedef struct _generic_N_Vector SunVector;
101#def typedef struct _N_VectorContent_Serial SunContent;
102
103#def typedef struct _generic_SUNMatrix SunMatrix;
104#def typedef struct _SUNMatrixContent_Dense SunMatrixContent;
105
106getContentMatrixPtr :: Storable a => Ptr b -> IO a
107getContentMatrixPtr ptr = (#peek SunMatrix, content) ptr
108
109getNRows :: Ptr b -> IO CInt
110getNRows ptr = (#peek SunMatrixContent, M) ptr
111putNRows :: CInt -> Ptr b -> IO ()
112putNRows nr ptr = (#poke SunMatrixContent, M) ptr nr
113
114getNCols :: Ptr b -> IO CInt
115getNCols ptr = (#peek SunMatrixContent, N) ptr
116putNCols :: CInt -> Ptr b -> IO ()
117putNCols nc ptr = (#poke SunMatrixContent, N) ptr nc
118
119getMatrixData :: Storable a => Ptr b -> IO a
120getMatrixData ptr = (#peek SunMatrixContent, data) ptr
121
122getContentPtr :: Storable a => Ptr b -> IO a
123getContentPtr ptr = (#peek SunVector, content) ptr
124
125getData :: Storable a => Ptr b -> IO a
126getData ptr = (#peek SunContent, data) ptr
127
128cV_ADAMS :: Int
129cV_ADAMS = #const CV_ADAMS
130cV_BDF :: Int
131cV_BDF = #const CV_BDF
132
133arkSMax :: Int
134arkSMax = #const ARK_S_MAX
135
136mIN_DIRK_NUM, mAX_DIRK_NUM :: Int
137mIN_DIRK_NUM = #const MIN_DIRK_NUM
138mAX_DIRK_NUM = #const MAX_DIRK_NUM
139
140-- FIXME: We could just use inline-c instead
141
142-- Butcher table accessors -- implicit
143sDIRK_2_1_2 :: Int
144sDIRK_2_1_2 = #const SDIRK_2_1_2
145bILLINGTON_3_3_2 :: Int
146bILLINGTON_3_3_2 = #const BILLINGTON_3_3_2
147tRBDF2_3_3_2 :: Int
148tRBDF2_3_3_2 = #const TRBDF2_3_3_2
149kVAERNO_4_2_3 :: Int
150kVAERNO_4_2_3 = #const KVAERNO_4_2_3
151aRK324L2SA_DIRK_4_2_3 :: Int
152aRK324L2SA_DIRK_4_2_3 = #const ARK324L2SA_DIRK_4_2_3
153cASH_5_2_4 :: Int
154cASH_5_2_4 = #const CASH_5_2_4
155cASH_5_3_4 :: Int
156cASH_5_3_4 = #const CASH_5_3_4
157sDIRK_5_3_4 :: Int
158sDIRK_5_3_4 = #const SDIRK_5_3_4
159kVAERNO_5_3_4 :: Int
160kVAERNO_5_3_4 = #const KVAERNO_5_3_4
161aRK436L2SA_DIRK_6_3_4 :: Int
162aRK436L2SA_DIRK_6_3_4 = #const ARK436L2SA_DIRK_6_3_4
163kVAERNO_7_4_5 :: Int
164kVAERNO_7_4_5 = #const KVAERNO_7_4_5
165aRK548L2SA_DIRK_8_4_5 :: Int
166aRK548L2SA_DIRK_8_4_5 = #const ARK548L2SA_DIRK_8_4_5
167
168-- #define DEFAULT_DIRK_2 SDIRK_2_1_2
169-- #define DEFAULT_DIRK_3 ARK324L2SA_DIRK_4_2_3
170-- #define DEFAULT_DIRK_4 SDIRK_5_3_4
171-- #define DEFAULT_DIRK_5 ARK548L2SA_DIRK_8_4_5
172
173-- Butcher table accessors -- explicit
174hEUN_EULER_2_1_2 :: Int
175hEUN_EULER_2_1_2 = #const HEUN_EULER_2_1_2
176bOGACKI_SHAMPINE_4_2_3 :: Int
177bOGACKI_SHAMPINE_4_2_3 = #const BOGACKI_SHAMPINE_4_2_3
178aRK324L2SA_ERK_4_2_3 :: Int
179aRK324L2SA_ERK_4_2_3 = #const ARK324L2SA_ERK_4_2_3
180zONNEVELD_5_3_4 :: Int
181zONNEVELD_5_3_4 = #const ZONNEVELD_5_3_4
182aRK436L2SA_ERK_6_3_4 :: Int
183aRK436L2SA_ERK_6_3_4 = #const ARK436L2SA_ERK_6_3_4
184sAYFY_ABURUB_6_3_4 :: Int
185sAYFY_ABURUB_6_3_4 = #const SAYFY_ABURUB_6_3_4
186cASH_KARP_6_4_5 :: Int
187cASH_KARP_6_4_5 = #const CASH_KARP_6_4_5
188fEHLBERG_6_4_5 :: Int
189fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5
190dORMAND_PRINCE_7_4_5 :: Int
191dORMAND_PRINCE_7_4_5 = #const DORMAND_PRINCE_7_4_5
192aRK548L2SA_ERK_8_4_5 :: Int
193aRK548L2SA_ERK_8_4_5 = #const ARK548L2SA_ERK_8_4_5
194vERNER_8_5_6 :: Int
195vERNER_8_5_6 = #const VERNER_8_5_6
196fEHLBERG_13_7_8 :: Int
197fEHLBERG_13_7_8 = #const FEHLBERG_13_7_8
198
199-- #define DEFAULT_ERK_2 HEUN_EULER_2_1_2
200-- #define DEFAULT_ERK_3 BOGACKI_SHAMPINE_4_2_3
201-- #define DEFAULT_ERK_4 ZONNEVELD_5_3_4
202-- #define DEFAULT_ERK_5 CASH_KARP_6_4_5
203-- #define DEFAULT_ERK_6 VERNER_8_5_6
204-- #define DEFAULT_ERK_8 FEHLBERG_13_7_8
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
new file mode 100644
index 0000000..a6f185e
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
@@ -0,0 +1,476 @@
1{-# OPTIONS_GHC -Wall #-}
2
3{-# LANGUAGE QuasiQuotes #-}
4{-# LANGUAGE TemplateHaskell #-}
5{-# LANGUAGE MultiWayIf #-}
6{-# LANGUAGE OverloadedStrings #-}
7{-# LANGUAGE ScopedTypeVariables #-}
8
9-----------------------------------------------------------------------------
10-- |
11-- Module : Numeric.Sundials.CVode.ODE
12-- Copyright : Dominic Steinitz 2018,
13-- Novadiscovery 2018
14-- License : BSD
15-- Maintainer : Dominic Steinitz
16-- Stability : provisional
17--
18-- Solution of ordinary differential equation (ODE) initial value problems.
19--
20-- <https://computation.llnl.gov/projects/sundials/sundials-software>
21--
22-- A simple example:
23--
24-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
25--
26-- @
27-- import Numeric.Sundials.CVode.ODE
28-- import Numeric.LinearAlgebra
29--
30-- import Plots as P
31-- import qualified Diagrams.Prelude as D
32-- import Diagrams.Backend.Rasterific
33--
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
46--
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
56--
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))
62-- @
63--
64-----------------------------------------------------------------------------
65module Numeric.Sundials.CVode.ODE ( odeSolve
66 , odeSolveV
67 , odeSolveVWith
68 , odeSolveVWith'
69 , ODEMethod(..)
70 , StepControl(..)
71 ) where
72
73import qualified Language.C.Inline as C
74import qualified Language.C.Inline.Unsafe as CU
75
76import Data.Monoid ((<>))
77import Data.Maybe (isJust)
78
79import Foreign.C.Types (CDouble, CInt, CLong)
80import Foreign.Ptr (Ptr)
81import Foreign.Storable (poke)
82
83import qualified Data.Vector.Storable as V
84
85import Data.Coerce (coerce)
86import System.IO.Unsafe (unsafePerformIO)
87
88import Numeric.LinearAlgebra.Devel (createVector)
89
90import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
91 cols, toLists, size, reshape)
92
93import Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF,
94 getDataFromContents, putDataInContents)
95import qualified Numeric.Sundials.Arkode as T
96import Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..))
97
98
99C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
100
101C.include "<stdlib.h>"
102C.include "<stdio.h>"
103C.include "<math.h>"
104C.include "<cvode/cvode.h>" -- prototypes for CVODE fcts., consts.
105C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros
106C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
107C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
108C.include "<cvode/cvode_direct.h>" -- access to CVDls interface
109C.include "<sundials/sundials_types.h>" -- definition of type realtype
110C.include "<sundials/sundials_math.h>"
111C.include "../../../helpers.h"
112C.include "Numeric/Sundials/Arkode_hsc.h"
113
114
115-- | Stepping functions
116data ODEMethod = ADAMS
117 | BDF
118
119getMethod :: ODEMethod -> Int
120getMethod (ADAMS) = cV_ADAMS
121getMethod (BDF) = cV_BDF
122
123getJacobian :: ODEMethod -> Maybe Jacobian
124getJacobian _ = Nothing
125
126-- | A version of 'odeSolveVWith' with reasonable default step control.
127odeSolveV
128 :: ODEMethod
129 -> Maybe Double -- ^ initial step size - by default, CVode
130 -- estimates the initial step size to be the
131 -- solution \(h\) of the equation
132 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
133 -- \(\ddot{y}\) is an estimated value of the
134 -- second derivative of the solution at \(t_0\)
135 -> Double -- ^ absolute tolerance for the state vector
136 -> Double -- ^ relative tolerance for the state vector
137 -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
138 -> Vector Double -- ^ initial conditions
139 -> Vector Double -- ^ desired solution times
140 -> Matrix Double -- ^ solution
141odeSolveV meth hi epsAbs epsRel f y0 ts =
142 odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts
143 where
144 g t x0 = coerce $ f t x0
145
146-- | A version of 'odeSolveV' with reasonable default parameters and
147-- system of equations defined using lists. FIXME: we should say
148-- something about the fact we could use the Jacobian but don't for
149-- compatibility with hmatrix-gsl.
150odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
151 -> [Double] -- ^ initial conditions
152 -> Vector Double -- ^ desired solution times
153 -> Matrix Double -- ^ solution
154odeSolve f y0 ts =
155 -- FIXME: These tolerances are different from the ones in GSL
156 odeSolveVWith BDF (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts)
157 where
158 g t x0 = V.fromList $ f t (V.toList x0)
159
160odeSolveVWith ::
161 ODEMethod
162 -> StepControl
163 -> Maybe Double -- ^ initial step size - by default, CVode
164 -- estimates the initial step size to be the
165 -- solution \(h\) of the equation
166 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
167 -- \(\ddot{y}\) is an estimated value of the second
168 -- derivative of the solution at \(t_0\)
169 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
170 -> V.Vector Double -- ^ Initial conditions
171 -> V.Vector Double -- ^ Desired solution times
172 -> Matrix Double -- ^ Error code or solution
173odeSolveVWith method control initStepSize f y0 tt =
174 case odeSolveVWith' opts method control initStepSize f y0 tt of
175 Left c -> error $ show c -- FIXME
176 Right (v, _d) -> v
177 where
178 opts = ODEOpts { maxNumSteps = 10000
179 , minStep = 1.0e-12
180 , relTol = error "relTol"
181 , absTols = error "absTol"
182 , initStep = error "initStep"
183 , maxFail = 10
184 }
185
186odeSolveVWith' ::
187 ODEOpts
188 -> ODEMethod
189 -> StepControl
190 -> Maybe Double -- ^ initial step size - by default, CVode
191 -- estimates the initial step size to be the
192 -- solution \(h\) of the equation
193 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
194 -- \(\ddot{y}\) is an estimated value of the second
195 -- derivative of the solution at \(t_0\)
196 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
197 -> V.Vector Double -- ^ Initial conditions
198 -> V.Vector Double -- ^ Desired solution times
199 -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
200odeSolveVWith' opts method control initStepSize f y0 tt =
201 case solveOdeC (fromIntegral $ maxFail opts)
202 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
203 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
204 (coerce f) (coerce y0) (coerce tt) of
205 Left c -> Left $ fromIntegral c
206 Right (v, d) -> Right (reshape l (coerce v), d)
207 where
208 l = size y0
209 scise (X aTol rTol) = coerce (V.replicate l aTol, rTol)
210 scise (X' aTol rTol) = coerce (V.replicate l aTol, rTol)
211 scise (XX' aTol rTol yScale _yDotScale) = coerce (V.replicate l aTol, yScale * rTol)
212 -- FIXME; Should we check that the length of ss is correct?
213 scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
214 jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
215 getJacobian method
216 matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
217 where
218 nr = fromIntegral $ rows m
219 nc = fromIntegral $ cols m
220 -- FIXME: efficiency
221 vs = V.fromList $ map coerce $ concat $ toLists m
222
223solveOdeC ::
224 CInt ->
225 CLong ->
226 CDouble ->
227 CInt ->
228 Maybe CDouble ->
229 (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
230 (V.Vector CDouble, CDouble) ->
231 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
232 -> V.Vector CDouble -- ^ Initial conditions
233 -> V.Vector CDouble -- ^ Desired solution times
234 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
235solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
236 jacH (aTols, rTol) fun f0 ts =
237 unsafePerformIO $ do
238
239 let isInitStepSize :: CInt
240 isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
241 ss :: CDouble
242 ss = case initStepSize of
243 -- It would be better to put an error message here but
244 -- inline-c seems to evaluate this even if it is never
245 -- used :(
246 Nothing -> 0.0
247 Just x -> x
248
249 let dim = V.length f0
250 nEq :: CLong
251 nEq = fromIntegral dim
252 nTs :: CInt
253 nTs = fromIntegral $ V.length ts
254 quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
255 qMatMut <- V.thaw quasiMatrixRes
256 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
257 diagMut <- V.thaw diagnostics
258 -- We need the types that sundials expects. These are tied together
259 -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
260 let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
261 funIO x y f _ptr = do
262 -- Convert the pointer we get from C (y) to a vector, and then
263 -- apply the user-supplied function.
264 fImm <- fun x <$> getDataFromContents dim y
265 -- Fill in the provided pointer with the resulting vector.
266 putDataInContents fImm dim f
267 -- FIXME: I don't understand what this comment means
268 -- Unsafe since the function will be called many times.
269 [CU.exp| int{ 0 } |]
270 let isJac :: CInt
271 isJac = fromIntegral $ fromEnum $ isJust jacH
272 jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
273 Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
274 IO CInt
275 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
276 case jacH of
277 Nothing -> error "Numeric.Sundials.CVode.ODE: Jacobian not defined"
278 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
279 poke jacS j
280 -- FIXME: I don't understand what this comment means
281 -- Unsafe since the function will be called many times.
282 [CU.exp| int{ 0 } |]
283
284 res <- [C.block| int {
285 /* general problem variables */
286
287 int flag; /* reusable error-checking flag */
288 int i, j; /* reusable loop indices */
289 N_Vector y = NULL; /* empty vector for storing solution */
290 N_Vector tv = NULL; /* empty vector for storing absolute tolerances */
291
292 SUNMatrix A = NULL; /* empty matrix for linear solver */
293 SUNLinearSolver LS = NULL; /* empty linear solver object */
294 void *cvode_mem = NULL; /* empty CVODE memory structure */
295 realtype t;
296 long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
297
298 /* general problem parameters */
299
300 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
301 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */
302
303 /* Initialize data structures */
304
305 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
306 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
307 /* Specify initial condition */
308 for (i = 0; i < NEQ; i++) {
309 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
310 };
311
312 cvode_mem = CVodeCreate($(int method), CV_NEWTON);
313 if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
314
315 /* Call CVodeInit to initialize the integrator memory and specify the
316 * user's right hand side function in y'=f(t,y), the inital time T0, and
317 * the initial dependent variable vector y. */
318 flag = CVodeInit(cvode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
319 if (check_flag(&flag, "CVodeInit", 1)) return(1);
320
321 tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
322 if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
323 /* Specify tolerances */
324 for (i = 0; i < NEQ; i++) {
325 NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
326 };
327
328 flag = CVodeSetMinStep(cvode_mem, $(double minStep_));
329 if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1;
330 flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_));
331 if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1;
332 flag = CVodeSetMaxErrTestFails(cvode_mem, $(int maxErrTestFails));
333 if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return 1;
334
335 /* Call CVodeSVtolerances to specify the scalar relative tolerance
336 * and vector absolute tolerances */
337 flag = CVodeSVtolerances(cvode_mem, $(double rTol), tv);
338 if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
339
340 /* Initialize dense matrix data structure and solver */
341 A = SUNDenseMatrix(NEQ, NEQ);
342 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
343 LS = SUNDenseLinearSolver(y, A);
344 if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
345
346 /* Attach matrix and linear solver */
347 flag = CVDlsSetLinearSolver(cvode_mem, LS, A);
348 if (check_flag(&flag, "CVDlsSetLinearSolver", 1)) return 1;
349
350 /* Set the initial step size if there is one */
351 if ($(int isInitStepSize)) {
352 /* FIXME: We could check if the initial step size is 0 */
353 /* or even NaN and then throw an error */
354 flag = CVodeSetInitStep(cvode_mem, $(double ss));
355 if (check_flag(&flag, "CVodeSetInitStep", 1)) return 1;
356 }
357
358 /* Set the Jacobian if there is one */
359 if ($(int isJac)) {
360 flag = CVDlsSetJacFn(cvode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
361 if (check_flag(&flag, "CVDlsSetJacFn", 1)) return 1;
362 }
363
364 /* Store initial conditions */
365 for (j = 0; j < NEQ; j++) {
366 ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
367 }
368
369 /* Main time-stepping loop: calls CVode to perform the integration */
370 /* Stops when the final time has been reached */
371 for (i = 1; i < $(int nTs); i++) {
372
373 flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[i], y, &t, CV_NORMAL); /* call integrator */
374 if (check_flag(&flag, "CVode", 1)) break;
375
376 /* Store the results for Haskell */
377 for (j = 0; j < NEQ; j++) {
378 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
379 }
380
381 /* unsuccessful solve: break */
382 if (flag < 0) {
383 fprintf(stderr,"Solver failure, stopping integration\n");
384 break;
385 }
386 }
387
388 /* Get some final statistics on how the solve progressed */
389
390 flag = CVodeGetNumSteps(cvode_mem, &nst);
391 check_flag(&flag, "CVodeGetNumSteps", 1);
392 ($vec-ptr:(long int *diagMut))[0] = nst;
393
394 /* FIXME */
395 ($vec-ptr:(long int *diagMut))[1] = 0;
396
397 flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
398 check_flag(&flag, "CVodeGetNumRhsEvals", 1);
399 ($vec-ptr:(long int *diagMut))[2] = nfe;
400 /* FIXME */
401 ($vec-ptr:(long int *diagMut))[3] = 0;
402
403 flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
404 check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
405 ($vec-ptr:(long int *diagMut))[4] = nsetups;
406
407 flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
408 check_flag(&flag, "CVodeGetNumErrTestFails", 1);
409 ($vec-ptr:(long int *diagMut))[5] = netf;
410
411 flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
412 check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
413 ($vec-ptr:(long int *diagMut))[6] = nni;
414
415 flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
416 check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
417 ($vec-ptr:(long int *diagMut))[7] = ncfn;
418
419 flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
420 check_flag(&flag, "CVDlsGetNumJacEvals", 1);
421 ($vec-ptr:(long int *diagMut))[8] = ncfn;
422
423 flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
424 check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
425 ($vec-ptr:(long int *diagMut))[9] = ncfn;
426
427 /* Clean up and return */
428
429 N_VDestroy(y); /* Free y vector */
430 N_VDestroy(tv); /* Free tv vector */
431 CVodeFree(&cvode_mem); /* Free integrator memory */
432 SUNLinSolFree(LS); /* Free linear solver */
433 SUNMatDestroy(A); /* Free A matrix */
434
435 return flag;
436 } |]
437 if res == 0
438 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)
452 else do
453 return $ Left res
454
455-- | Adaptive step-size control
456-- functions.
457--
458-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
459-- allows the user to control the step size adjustment using
460-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
461-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
462-- is the required relative error, \(s_i\) is a vector of scaling
463-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
464-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
465--
466-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
467-- allows the user to control the step size adjustment using
468-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
469-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
470-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
471-- effect.
472data StepControl = X Double Double -- ^ absolute and relative tolerance for \(y\); in GSL terms, \(a_{y} = 1\) and \(a_{dy/dt} = 0\); in ARKode terms, the \(\eta^{abs}_i\) are identical
473 | X' Double Double -- ^ absolute and relative tolerance for \(\dot{y}\); in GSL terms, \(a_{y} = 0\) and \(a_{dy/dt} = 1\); in ARKode terms, the latter is treated as the relative tolerance for \(y\) so this is the same as specifying 'X' which may be entirely incorrect for the given problem
474 | XX' Double Double Double Double -- ^ include both via relative tolerance
475 -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
476 | ScXX' Double Double Double Double (Vector Double) -- ^ scale absolute tolerance of \(y_i\); in ARKode terms, \(a_{{dy}/{dt}}\) is ignored, \(\eta^{abs}_i = s_i \epsilon^{abs}\) and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
diff --git a/packages/sundials/src/Numeric/Sundials/ODEOpts.hs b/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
new file mode 100644
index 0000000..027d99a
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
@@ -0,0 +1,32 @@
1module Numeric.Sundials.ODEOpts where
2
3import Data.Word (Word32)
4import qualified Data.Vector.Storable as VS
5
6import Numeric.LinearAlgebra.HMatrix (Vector, Matrix)
7
8
9type Jacobian = Double -> Vector Double -> Matrix Double
10
11data ODEOpts = ODEOpts {
12 maxNumSteps :: Word32
13 , minStep :: Double
14 , relTol :: Double
15 , absTols :: VS.Vector Double
16 , initStep :: Maybe Double
17 , maxFail :: Word32
18 } deriving (Read, Show, Eq, Ord)
19
20data SundialsDiagnostics = SundialsDiagnostics {
21 aRKodeGetNumSteps :: Int
22 , aRKodeGetNumStepAttempts :: Int
23 , aRKodeGetNumRhsEvals_fe :: Int
24 , aRKodeGetNumRhsEvals_fi :: Int
25 , aRKodeGetNumLinSolvSetups :: Int
26 , aRKodeGetNumErrTestFails :: Int
27 , aRKodeGetNumNonlinSolvIters :: Int
28 , aRKodeGetNumNonlinSolvConvFails :: Int
29 , aRKDlsGetNumJacEvals :: Int
30 , aRKDlsGetNumRhsEvals :: Int
31 } deriving Show
32
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
new file mode 100644
index 0000000..f0ca592
--- /dev/null
+++ b/packages/sundials/src/helpers.c
@@ -0,0 +1,44 @@
1#include <stdio.h>
2#include <math.h>
3#include <arkode/arkode.h> /* prototypes for ARKODE fcts., consts. */
4#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
5#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
6#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
7#include <arkode/arkode_direct.h> /* access to ARKDls interface */
8#include <sundials/sundials_types.h> /* definition of type realtype */
9#include <sundials/sundials_math.h>
10
11/* Check function return value...
12 opt == 0 means SUNDIALS function allocates memory so check if
13 returned NULL pointer
14 opt == 1 means SUNDIALS function returns a flag so check if
15 flag >= 0
16 opt == 2 means function allocates memory so check if returned
17 NULL pointer
18*/
19int check_flag(void *flagvalue, const char *funcname, int opt)
20{
21 int *errflag;
22
23 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
24 if (opt == 0 && flagvalue == NULL) {
25 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
26 funcname);
27 return 1; }
28
29 /* Check if flag < 0 */
30 else if (opt == 1) {
31 errflag = (int *) flagvalue;
32 if (*errflag < 0) {
33 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
34 funcname, *errflag);
35 return 1; }}
36
37 /* Check if function returned NULL pointer - no memory allocated */
38 else if (opt == 2 && flagvalue == NULL) {
39 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
40 funcname);
41 return 1; }
42
43 return 0;
44}
diff --git a/packages/sundials/src/helpers.h b/packages/sundials/src/helpers.h
new file mode 100644
index 0000000..3d8fbc0
--- /dev/null
+++ b/packages/sundials/src/helpers.h
@@ -0,0 +1,9 @@
1/* Check function return value...
2 opt == 0 means SUNDIALS function allocates memory so check if
3 returned NULL pointer
4 opt == 1 means SUNDIALS function returns a flag so check if
5 flag >= 0
6 opt == 2 means function allocates memory so check if returned
7 NULL pointer
8*/
9int check_flag(void *flagvalue, const char *funcname, int opt);
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-30 03:32:22 +0000