summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2019-06-30 08:58:16 +0100
committerDominic Steinitz <dominic@steinitz.org>2019-06-30 08:58:16 +0100
commitcb09d0e99ae4f10cd2b3f3ac667df83946a9744d (patch)
tree7469898bbde43392e2a4306450139cd434edba0d /packages
parentc4b80ef9951b533d6bbbb34df8109f3290546296 (diff)
Remove sundials as it has its own repo now. Fix
* https://github.com/haskell-numerics/hmatrix/issues/304 * https://github.com/haskell-numerics/hmatrix/issues/302 * https://github.com/haskell-numerics/hmatrix/issues/290
Diffstat (limited to 'packages')
-rw-r--r--packages/base/hmatrix.cabal4
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Static.hs2
-rw-r--r--packages/base/stack.yaml3
-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.pngbin27362 -> 0 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.hs896
-rw-r--r--packages/sundials/src/Numeric/Sundials/Arkode.hsc204
-rw-r--r--packages/sundials/src/Numeric/Sundials/CVode/ODE.hs471
-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
16 files changed, 5 insertions, 1952 deletions
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 1380a36..76e3825 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -1,5 +1,5 @@
1Name: hmatrix 1Name: hmatrix
2Version: 0.19.0.0 2Version: 0.20.0.0
3License: BSD3 3License: BSD3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
@@ -38,7 +38,7 @@ flag disable-default-paths
38 38
39library 39library
40 40
41 Build-Depends: base >= 4.8 && < 5, 41 Build-Depends: base >= 4.9 && < 5,
42 binary, 42 binary,
43 array, 43 array,
44 deepseq, 44 deepseq,
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 2e05c90..e5ce4e2 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -92,7 +92,7 @@ ud1 (R (Dim v)) = v
92 92
93 93
94infixl 4 & 94infixl 4 &
95(&) :: forall n . (KnownNat n, 1 <= n) 95(&) :: forall n . KnownNat n
96 => R n -> ℝ -> R (n+1) 96 => R n -> ℝ -> R (n+1)
97u & x = u # (konst x :: R 1) 97u & x = u # (konst x :: R 1)
98 98
diff --git a/packages/base/stack.yaml b/packages/base/stack.yaml
index f4001c6..28fc91d 100644
--- a/packages/base/stack.yaml
+++ b/packages/base/stack.yaml
@@ -4,4 +4,5 @@ flags:
4packages: 4packages:
5- '.' 5- '.'
6extra-deps: [] 6extra-deps: []
7resolver: lts-3.3 7resolver: lts-13.21
8
diff --git a/packages/sundials/ChangeLog.md b/packages/sundials/ChangeLog.md
deleted file mode 100644
index 7b15777..0000000
--- a/packages/sundials/ChangeLog.md
+++ /dev/null
@@ -1,5 +0,0 @@
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
deleted file mode 100644
index a162e98..0000000
--- a/packages/sundials/LICENSE
+++ /dev/null
@@ -1,30 +0,0 @@
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
deleted file mode 100644
index 2fac5c2..0000000
--- a/packages/sundials/README.md
+++ /dev/null
@@ -1,8 +0,0 @@
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
deleted file mode 100644
index 9a994af..0000000
--- a/packages/sundials/Setup.hs
+++ /dev/null
@@ -1,2 +0,0 @@
1import Distribution.Simple
2main = defaultMain
diff --git a/packages/sundials/diagrams/brusselator.png b/packages/sundials/diagrams/brusselator.png
deleted file mode 100644
index 740cacb..0000000
--- a/packages/sundials/diagrams/brusselator.png
+++ /dev/null
Binary files differ
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
deleted file mode 100644
index cd2be4e..0000000
--- a/packages/sundials/hmatrix-sundials.cabal
+++ /dev/null
@@ -1,61 +0,0 @@
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
deleted file mode 100644
index 16c21c5..0000000
--- a/packages/sundials/src/Main.hs
+++ /dev/null
@@ -1,186 +0,0 @@
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
deleted file mode 100644
index 48ac887..0000000
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ /dev/null
@@ -1,896 +0,0 @@
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, _v) -> 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 (Matrix Double, 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 (v, c) -> Left (reshape l (coerce v), 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 (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or
491 -- solution and diagnostics
492solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
493 jacH (aTols, rTol) fun f0 ts = unsafePerformIO $ do
494
495 let isInitStepSize :: CInt
496 isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
497 ss :: CDouble
498 ss = case initStepSize of
499 -- It would be better to put an error message here but
500 -- inline-c seems to evaluate this even if it is never
501 -- used :(
502 Nothing -> 0.0
503 Just x -> x
504
505 let dim = V.length f0
506 nEq :: CLong
507 nEq = fromIntegral dim
508 nTs :: CInt
509 nTs = fromIntegral $ V.length ts
510 quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
511 qMatMut <- V.thaw quasiMatrixRes
512 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
513 diagMut <- V.thaw diagnostics
514 -- We need the types that sundials expects. These are tied together
515 -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
516 let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
517 funIO x y f _ptr = do
518 -- Convert the pointer we get from C (y) to a vector, and then
519 -- apply the user-supplied function.
520 fImm <- fun x <$> getDataFromContents dim y
521 -- Fill in the provided pointer with the resulting vector.
522 putDataInContents fImm dim f
523 -- FIXME: I don't understand what this comment means
524 -- Unsafe since the function will be called many times.
525 [CU.exp| int{ 0 } |]
526 let isJac :: CInt
527 isJac = fromIntegral $ fromEnum $ isJust jacH
528 jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
529 Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
530 IO CInt
531 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
532 case jacH of
533 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
534 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
535 poke jacS j
536 -- FIXME: I don't understand what this comment means
537 -- Unsafe since the function will be called many times.
538 [CU.exp| int{ 0 } |]
539
540 res <- [C.block| int {
541 /* general problem variables */
542
543 int flag; /* reusable error-checking flag */
544 int i, j; /* reusable loop indices */
545 N_Vector y = NULL; /* empty vector for storing solution */
546 N_Vector tv = NULL; /* empty vector for storing absolute tolerances */
547 SUNMatrix A = NULL; /* empty matrix for linear solver */
548 SUNLinearSolver LS = NULL; /* empty linear solver object */
549 void *arkode_mem = NULL; /* empty ARKode memory structure */
550 realtype t;
551 long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
552
553 /* general problem parameters */
554
555 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
556 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */
557
558 /* Initialize data structures */
559
560 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
561 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
562 /* Specify initial condition */
563 for (i = 0; i < NEQ; i++) {
564 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
565 };
566
567 tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
568 if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
569 /* Specify tolerances */
570 for (i = 0; i < NEQ; i++) {
571 NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
572 };
573
574 arkode_mem = ARKodeCreate(); /* Create the solver memory */
575 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
576
577 /* Call ARKodeInit to initialize the integrator memory and specify the */
578 /* right-hand side function in y'=f(t,y), the inital time T0, and */
579 /* the initial dependent variable vector y. Note: we treat the */
580 /* problem as fully implicit and set f_E to NULL and f_I to f. */
581
582 /* Here we use the C types defined in helpers.h which tie up with */
583 /* the Haskell types defined in CLangToHaskellTypes */
584 if ($(int method) < MIN_DIRK_NUM) {
585 flag = ARKodeInit(arkode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y);
586 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
587 } else {
588 flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
589 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
590 }
591
592 flag = ARKodeSetMinStep(arkode_mem, $(double minStep_));
593 if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1;
594 flag = ARKodeSetMaxNumSteps(arkode_mem, $(long int maxNumSteps_));
595 if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1;
596 flag = ARKodeSetMaxErrTestFails(arkode_mem, $(int maxErrTestFails));
597 if (check_flag(&flag, "ARKodeSetMaxErrTestFails", 1)) return 1;
598
599 /* Set routines */
600 flag = ARKodeSVtolerances(arkode_mem, $(double rTol), tv);
601 if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1;
602
603 /* Initialize dense matrix data structure and solver */
604 A = SUNDenseMatrix(NEQ, NEQ);
605 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
606 LS = SUNDenseLinearSolver(y, A);
607 if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
608
609 /* Attach matrix and linear solver */
610 flag = ARKDlsSetLinearSolver(arkode_mem, LS, A);
611 if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
612
613 /* Set the initial step size if there is one */
614 if ($(int isInitStepSize)) {
615 /* FIXME: We could check if the initial step size is 0 */
616 /* or even NaN and then throw an error */
617 flag = ARKodeSetInitStep(arkode_mem, $(double ss));
618 if (check_flag(&flag, "ARKodeSetInitStep", 1)) return 1;
619 }
620
621 /* Set the Jacobian if there is one */
622 if ($(int isJac)) {
623 flag = ARKDlsSetJacFn(arkode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
624 if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
625 }
626
627 /* Store initial conditions */
628 for (j = 0; j < NEQ; j++) {
629 ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
630 }
631
632 /* Explicitly set the method */
633 if ($(int method) >= MIN_DIRK_NUM) {
634 flag = ARKodeSetIRKTableNum(arkode_mem, $(int method));
635 if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
636 } else {
637 flag = ARKodeSetERKTableNum(arkode_mem, $(int method));
638 if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
639 }
640
641 /* Main time-stepping loop: calls ARKode to perform the integration */
642 /* Stops when the final time has been reached */
643 for (i = 1; i < $(int nTs); i++) {
644
645 flag = ARKode(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */
646 if (check_flag(&flag, "ARKode solver failure, stopping integration", 1)) return 1;
647
648 /* Store the results for Haskell */
649 for (j = 0; j < NEQ; j++) {
650 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
651 }
652 }
653
654 /* Get some final statistics on how the solve progressed */
655
656 flag = ARKodeGetNumSteps(arkode_mem, &nst);
657 check_flag(&flag, "ARKodeGetNumSteps", 1);
658 ($vec-ptr:(long int *diagMut))[0] = nst;
659
660 flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
661 check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
662 ($vec-ptr:(long int *diagMut))[1] = nst_a;
663
664 flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
665 check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
666 ($vec-ptr:(long int *diagMut))[2] = nfe;
667 ($vec-ptr:(long int *diagMut))[3] = nfi;
668
669 flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
670 check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
671 ($vec-ptr:(long int *diagMut))[4] = nsetups;
672
673 flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
674 check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
675 ($vec-ptr:(long int *diagMut))[5] = netf;
676
677 flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
678 check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
679 ($vec-ptr:(long int *diagMut))[6] = nni;
680
681 flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
682 check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
683 ($vec-ptr:(long int *diagMut))[7] = ncfn;
684
685 flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
686 check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
687 ($vec-ptr:(long int *diagMut))[8] = ncfn;
688
689 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
690 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
691 ($vec-ptr:(long int *diagMut))[9] = ncfn;
692
693 /* Clean up and return */
694 N_VDestroy(y); /* Free y vector */
695 N_VDestroy(tv); /* Free tv vector */
696 ARKodeFree(&arkode_mem); /* Free integrator memory */
697 SUNLinSolFree(LS); /* Free linear solver */
698 SUNMatDestroy(A); /* Free A matrix */
699
700 return flag;
701 } |]
702 preD <- V.freeze diagMut
703 let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
704 (fromIntegral $ preD V.!1)
705 (fromIntegral $ preD V.!2)
706 (fromIntegral $ preD V.!3)
707 (fromIntegral $ preD V.!4)
708 (fromIntegral $ preD V.!5)
709 (fromIntegral $ preD V.!6)
710 (fromIntegral $ preD V.!7)
711 (fromIntegral $ preD V.!8)
712 (fromIntegral $ preD V.!9)
713 m <- V.freeze qMatMut
714 if res == 0
715 then do
716 return $ Right (m, d)
717 else do
718 return $ Left (m, res)
719
720data ButcherTable = ButcherTable { am :: Matrix Double
721 , cv :: Vector Double
722 , bv :: Vector Double
723 , b2v :: Vector Double
724 }
725 deriving Show
726
727data ButcherTable' a = ButcherTable' { am' :: V.Vector a
728 , cv' :: V.Vector a
729 , bv' :: V.Vector a
730 , b2v' :: V.Vector a
731 }
732 deriving Show
733
734butcherTable :: ODEMethod -> ButcherTable
735butcherTable method =
736 case getBT method of
737 Left c -> error $ show c -- FIXME
738 Right (ButcherTable' v w x y, sqp) ->
739 ButcherTable { am = subMatrix (0, 0) (s, s) $ (arkSMax >< arkSMax) (V.toList v)
740 , cv = subVector 0 s w
741 , bv = subVector 0 s x
742 , b2v = subVector 0 s y
743 }
744 where
745 s = fromIntegral $ sqp V.! 0
746
747getBT :: ODEMethod -> Either Int (ButcherTable' Double, V.Vector Int)
748getBT method = case getButcherTable method of
749 Left c ->
750 Left $ fromIntegral c
751 Right (ButcherTable' a b c d, sqp) ->
752 Right $ ( ButcherTable' (coerce a) (coerce b) (coerce c) (coerce d)
753 , V.map fromIntegral sqp )
754
755getButcherTable :: ODEMethod
756 -> Either CInt (ButcherTable' CDouble, V.Vector CInt)
757getButcherTable method = unsafePerformIO $ do
758 -- ARKode seems to want an ODE in order to set and then get the
759 -- Butcher tableau so here's one to keep it happy
760 let funI :: CDouble -> V.Vector CDouble -> V.Vector CDouble
761 funI _t ys = V.fromList [ ys V.! 0 ]
762 let funE :: CDouble -> V.Vector CDouble -> V.Vector CDouble
763 funE _t ys = V.fromList [ ys V.! 0 ]
764 f0 = V.fromList [ 1.0 ]
765 ts = V.fromList [ 0.0 ]
766 dim = V.length f0
767 nEq :: CLong
768 nEq = fromIntegral dim
769 mN :: CInt
770 mN = fromIntegral $ getMethod method
771
772 btSQP :: V.Vector CInt <- createVector 3
773 btSQPMut <- V.thaw btSQP
774 btAs :: V.Vector CDouble <- createVector (arkSMax * arkSMax)
775 btAsMut <- V.thaw btAs
776 btCs :: V.Vector CDouble <- createVector arkSMax
777 btBs :: V.Vector CDouble <- createVector arkSMax
778 btB2s :: V.Vector CDouble <- createVector arkSMax
779 btCsMut <- V.thaw btCs
780 btBsMut <- V.thaw btBs
781 btB2sMut <- V.thaw btB2s
782 let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
783 funIOI x y f _ptr = do
784 fImm <- funI x <$> getDataFromContents dim y
785 putDataInContents fImm dim f
786 -- FIXME: I don't understand what this comment means
787 -- Unsafe since the function will be called many times.
788 [CU.exp| int{ 0 } |]
789 let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
790 funIOE x y f _ptr = do
791 fImm <- funE 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 res <- [C.block| int {
797 /* general problem variables */
798
799 int flag; /* reusable error-checking flag */
800 N_Vector y = NULL; /* empty vector for storing solution */
801 void *arkode_mem = NULL; /* empty ARKode memory structure */
802 int i, j; /* reusable loop indices */
803
804 /* general problem parameters */
805
806 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
807 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars */
808
809 /* Initialize data structures */
810
811 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
812 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
813 /* Specify initial condition */
814 for (i = 0; i < NEQ; i++) {
815 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
816 };
817 arkode_mem = ARKodeCreate(); /* Create the solver memory */
818 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
819
820 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);
821 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
822
823 if ($(int mN) >= MIN_DIRK_NUM) {
824 flag = ARKodeSetIRKTableNum(arkode_mem, $(int mN));
825 if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
826 } else {
827 flag = ARKodeSetERKTableNum(arkode_mem, $(int mN));
828 if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
829 }
830
831 int s, q, p;
832 realtype *ai = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
833 realtype *ae = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
834 realtype *ci = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
835 realtype *ce = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
836 realtype *bi = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
837 realtype *be = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
838 realtype *b2i = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
839 realtype *b2e = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
840 flag = ARKodeGetCurrentButcherTables(arkode_mem, &s, &q, &p, ai, ae, ci, ce, bi, be, b2i, b2e);
841 if (check_flag(&flag, "ARKode", 1)) return 1;
842 $vec-ptr:(int *btSQPMut)[0] = s;
843 $vec-ptr:(int *btSQPMut)[1] = q;
844 $vec-ptr:(int *btSQPMut)[2] = p;
845 for (i = 0; i < s; i++) {
846 for (j = 0; j < s; j++) {
847 /* FIXME: double should be realtype */
848 ($vec-ptr:(double *btAsMut))[i * ARK_S_MAX + j] = ai[i * ARK_S_MAX + j];
849 }
850 }
851
852 for (i = 0; i < s; i++) {
853 ($vec-ptr:(double *btCsMut))[i] = ci[i];
854 ($vec-ptr:(double *btBsMut))[i] = bi[i];
855 ($vec-ptr:(double *btB2sMut))[i] = b2i[i];
856 }
857
858 /* Clean up and return */
859 N_VDestroy(y); /* Free y vector */
860 ARKodeFree(&arkode_mem); /* Free integrator memory */
861
862 return flag;
863 } |]
864 if res == 0
865 then do
866 x <- V.freeze btAsMut
867 y <- V.freeze btSQPMut
868 z <- V.freeze btCsMut
869 u <- V.freeze btBsMut
870 v <- V.freeze btB2sMut
871 return $ Right (ButcherTable' { am' = x, cv' = z, bv' = u, b2v' = v }, y)
872 else do
873 return $ Left res
874
875-- | Adaptive step-size control
876-- functions.
877--
878-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
879-- allows the user to control the step size adjustment using
880-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
881-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
882-- is the required relative error, \(s_i\) is a vector of scaling
883-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
884-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
885--
886-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
887-- allows the user to control the step size adjustment using
888-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
889-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
890-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
891-- effect.
892data 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
893 | 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
894 | XX' Double Double Double Double -- ^ include both via relative tolerance
895 -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
896 | 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
deleted file mode 100644
index 0850258..0000000
--- a/packages/sundials/src/Numeric/Sundials/Arkode.hsc
+++ /dev/null
@@ -1,204 +0,0 @@
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
deleted file mode 100644
index ad7cf51..0000000
--- a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
+++ /dev/null
@@ -1,471 +0,0 @@
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, _v) -> 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 (Matrix Double, 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 (v, c) -> Left (reshape l (coerce v), 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 (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or
235 -- solution and diagnostics
236solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
237 jacH (aTols, rTol) fun f0 ts =
238 unsafePerformIO $ do
239
240 let isInitStepSize :: CInt
241 isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
242 ss :: CDouble
243 ss = case initStepSize of
244 -- It would be better to put an error message here but
245 -- inline-c seems to evaluate this even if it is never
246 -- used :(
247 Nothing -> 0.0
248 Just x -> x
249
250 let dim = V.length f0
251 nEq :: CLong
252 nEq = fromIntegral dim
253 nTs :: CInt
254 nTs = fromIntegral $ V.length ts
255 quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
256 qMatMut <- V.thaw quasiMatrixRes
257 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
258 diagMut <- V.thaw diagnostics
259 -- We need the types that sundials expects. These are tied together
260 -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
261 let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
262 funIO x y f _ptr = do
263 -- Convert the pointer we get from C (y) to a vector, and then
264 -- apply the user-supplied function.
265 fImm <- fun x <$> getDataFromContents dim y
266 -- Fill in the provided pointer with the resulting vector.
267 putDataInContents fImm dim f
268 -- FIXME: I don't understand what this comment means
269 -- Unsafe since the function will be called many times.
270 [CU.exp| int{ 0 } |]
271 let isJac :: CInt
272 isJac = fromIntegral $ fromEnum $ isJust jacH
273 jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
274 Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
275 IO CInt
276 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
277 case jacH of
278 Nothing -> error "Numeric.Sundials.CVode.ODE: Jacobian not defined"
279 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
280 poke jacS j
281 -- FIXME: I don't understand what this comment means
282 -- Unsafe since the function will be called many times.
283 [CU.exp| int{ 0 } |]
284
285 res <- [C.block| int {
286 /* general problem variables */
287
288 int flag; /* reusable error-checking flag */
289 int i, j; /* reusable loop indices */
290 N_Vector y = NULL; /* empty vector for storing solution */
291 N_Vector tv = NULL; /* empty vector for storing absolute tolerances */
292
293 SUNMatrix A = NULL; /* empty matrix for linear solver */
294 SUNLinearSolver LS = NULL; /* empty linear solver object */
295 void *cvode_mem = NULL; /* empty CVODE memory structure */
296 realtype t;
297 long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
298
299 /* general problem parameters */
300
301 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
302 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */
303
304 /* Initialize data structures */
305
306 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
307 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
308 /* Specify initial condition */
309 for (i = 0; i < NEQ; i++) {
310 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
311 };
312
313 cvode_mem = CVodeCreate($(int method), CV_NEWTON);
314 if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
315
316 /* Call CVodeInit to initialize the integrator memory and specify the
317 * user's right hand side function in y'=f(t,y), the inital time T0, and
318 * the initial dependent variable vector y. */
319 flag = CVodeInit(cvode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
320 if (check_flag(&flag, "CVodeInit", 1)) return(1);
321
322 tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
323 if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
324 /* Specify tolerances */
325 for (i = 0; i < NEQ; i++) {
326 NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
327 };
328
329 flag = CVodeSetMinStep(cvode_mem, $(double minStep_));
330 if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1;
331 flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_));
332 if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1;
333 flag = CVodeSetMaxErrTestFails(cvode_mem, $(int maxErrTestFails));
334 if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return 1;
335
336 /* Call CVodeSVtolerances to specify the scalar relative tolerance
337 * and vector absolute tolerances */
338 flag = CVodeSVtolerances(cvode_mem, $(double rTol), tv);
339 if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
340
341 /* Initialize dense matrix data structure and solver */
342 A = SUNDenseMatrix(NEQ, NEQ);
343 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
344 LS = SUNDenseLinearSolver(y, A);
345 if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
346
347 /* Attach matrix and linear solver */
348 flag = CVDlsSetLinearSolver(cvode_mem, LS, A);
349 if (check_flag(&flag, "CVDlsSetLinearSolver", 1)) return 1;
350
351 /* Set the initial step size if there is one */
352 if ($(int isInitStepSize)) {
353 /* FIXME: We could check if the initial step size is 0 */
354 /* or even NaN and then throw an error */
355 flag = CVodeSetInitStep(cvode_mem, $(double ss));
356 if (check_flag(&flag, "CVodeSetInitStep", 1)) return 1;
357 }
358
359 /* Set the Jacobian if there is one */
360 if ($(int isJac)) {
361 flag = CVDlsSetJacFn(cvode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
362 if (check_flag(&flag, "CVDlsSetJacFn", 1)) return 1;
363 }
364
365 /* Store initial conditions */
366 for (j = 0; j < NEQ; j++) {
367 ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
368 }
369
370 /* Main time-stepping loop: calls CVode to perform the integration */
371 /* Stops when the final time has been reached */
372 for (i = 1; i < $(int nTs); i++) {
373
374 flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[i], y, &t, CV_NORMAL); /* call integrator */
375 if (check_flag(&flag, "CVode solver failure, stopping integration", 1)) return 1;
376
377 /* Store the results for Haskell */
378 for (j = 0; j < NEQ; j++) {
379 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
380 }
381 }
382
383 /* Get some final statistics on how the solve progressed */
384
385 flag = CVodeGetNumSteps(cvode_mem, &nst);
386 check_flag(&flag, "CVodeGetNumSteps", 1);
387 ($vec-ptr:(long int *diagMut))[0] = nst;
388
389 /* FIXME */
390 ($vec-ptr:(long int *diagMut))[1] = 0;
391
392 flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
393 check_flag(&flag, "CVodeGetNumRhsEvals", 1);
394 ($vec-ptr:(long int *diagMut))[2] = nfe;
395 /* FIXME */
396 ($vec-ptr:(long int *diagMut))[3] = 0;
397
398 flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
399 check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
400 ($vec-ptr:(long int *diagMut))[4] = nsetups;
401
402 flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
403 check_flag(&flag, "CVodeGetNumErrTestFails", 1);
404 ($vec-ptr:(long int *diagMut))[5] = netf;
405
406 flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
407 check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
408 ($vec-ptr:(long int *diagMut))[6] = nni;
409
410 flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
411 check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
412 ($vec-ptr:(long int *diagMut))[7] = ncfn;
413
414 flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
415 check_flag(&flag, "CVDlsGetNumJacEvals", 1);
416 ($vec-ptr:(long int *diagMut))[8] = ncfn;
417
418 flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
419 check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
420 ($vec-ptr:(long int *diagMut))[9] = ncfn;
421
422 /* Clean up and return */
423
424 N_VDestroy(y); /* Free y vector */
425 N_VDestroy(tv); /* Free tv vector */
426 CVodeFree(&cvode_mem); /* Free integrator memory */
427 SUNLinSolFree(LS); /* Free linear solver */
428 SUNMatDestroy(A); /* Free A matrix */
429
430 return flag;
431 } |]
432 preD <- V.freeze diagMut
433 let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
434 (fromIntegral $ preD V.!1)
435 (fromIntegral $ preD V.!2)
436 (fromIntegral $ preD V.!3)
437 (fromIntegral $ preD V.!4)
438 (fromIntegral $ preD V.!5)
439 (fromIntegral $ preD V.!6)
440 (fromIntegral $ preD V.!7)
441 (fromIntegral $ preD V.!8)
442 (fromIntegral $ preD V.!9)
443 m <- V.freeze qMatMut
444 if res == 0
445 then do
446 return $ Right (m, d)
447 else do
448 return $ Left (m, res)
449
450-- | Adaptive step-size control
451-- functions.
452--
453-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
454-- allows the user to control the step size adjustment using
455-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
456-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
457-- is the required relative error, \(s_i\) is a vector of scaling
458-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
459-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
460--
461-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
462-- allows the user to control the step size adjustment using
463-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
464-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
465-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
466-- effect.
467data 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
468 | 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
469 | XX' Double Double Double Double -- ^ include both via relative tolerance
470 -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
471 | 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
deleted file mode 100644
index 027d99a..0000000
--- a/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
+++ /dev/null
@@ -1,32 +0,0 @@
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
deleted file mode 100644
index f0ca592..0000000
--- a/packages/sundials/src/helpers.c
+++ /dev/null
@@ -1,44 +0,0 @@
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
deleted file mode 100644
index 3d8fbc0..0000000
--- a/packages/sundials/src/helpers.h
+++ /dev/null
@@ -1,9 +0,0 @@
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-11-11 10:15:32 +0000