From cb09d0e99ae4f10cd2b3f3ac667df83946a9744d Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Sun, 30 Jun 2019 08:58:16 +0100 Subject: 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 --- examples/examples.cabal | 18 +- packages/base/hmatrix.cabal | 4 +- packages/base/src/Numeric/LinearAlgebra/Static.hs | 2 +- packages/base/stack.yaml | 3 +- packages/sundials/ChangeLog.md | 5 - packages/sundials/LICENSE | 30 - packages/sundials/README.md | 8 - packages/sundials/Setup.hs | 2 - packages/sundials/diagrams/brusselator.png | Bin 27362 -> 0 bytes packages/sundials/hmatrix-sundials.cabal | 61 -- packages/sundials/src/Main.hs | 186 ----- .../sundials/src/Numeric/Sundials/ARKode/ODE.hs | 896 --------------------- packages/sundials/src/Numeric/Sundials/Arkode.hsc | 204 ----- .../sundials/src/Numeric/Sundials/CVode/ODE.hs | 471 ----------- packages/sundials/src/Numeric/Sundials/ODEOpts.hs | 32 - packages/sundials/src/helpers.c | 44 - packages/sundials/src/helpers.h | 9 - stack.yaml | 1 - 18 files changed, 6 insertions(+), 1970 deletions(-) delete mode 100644 packages/sundials/ChangeLog.md delete mode 100644 packages/sundials/LICENSE delete mode 100644 packages/sundials/README.md delete mode 100644 packages/sundials/Setup.hs delete mode 100644 packages/sundials/diagrams/brusselator.png delete mode 100644 packages/sundials/hmatrix-sundials.cabal delete mode 100644 packages/sundials/src/Main.hs delete mode 100644 packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs delete mode 100644 packages/sundials/src/Numeric/Sundials/Arkode.hsc delete mode 100644 packages/sundials/src/Numeric/Sundials/CVode/ODE.hs delete mode 100644 packages/sundials/src/Numeric/Sundials/ODEOpts.hs delete mode 100644 packages/sundials/src/helpers.c delete mode 100644 packages/sundials/src/helpers.h diff --git a/examples/examples.cabal b/examples/examples.cabal index 286b2db..e647221 100644 --- a/examples/examples.cabal +++ b/examples/examples.cabal @@ -1,5 +1,5 @@ name: examples -version: 0.19.0.0 +version: 0.20.0.0 synopsis: Example usage of the various hmatrix packages homepage: https://github.com/albertoruiz/hmatrix license: BSD3 @@ -12,22 +12,6 @@ build-type: Simple extra-source-files: ChangeLog.md cabal-version: >=1.10 -executable sundials - main-is: sundials.hs - build-depends: base >=4.10 && <4.11, - hmatrix, - hmatrix-sundials, - hmatrix-gsl - default-language: Haskell2010 - -executable butcherTableau - main-is: ButcherTableau.hs - build-depends: base >=4.10 && <4.11, - hmatrix, - hmatrix-sundials, - pretty - default-language: Haskell2010 - executable vectorShow main-is: VectorShow.hs build-depends: base >=4.10 && <4.11, 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 @@ Name: hmatrix -Version: 0.19.0.0 +Version: 0.20.0.0 License: BSD3 License-file: LICENSE Author: Alberto Ruiz @@ -38,7 +38,7 @@ flag disable-default-paths library - Build-Depends: base >= 4.8 && < 5, + Build-Depends: base >= 4.9 && < 5, binary, array, 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 infixl 4 & -(&) :: forall n . (KnownNat n, 1 <= n) +(&) :: forall n . KnownNat n => R n -> ℝ -> R (n+1) u & x = u # (konst x :: R 1) 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: packages: - '.' extra-deps: [] -resolver: lts-3.3 +resolver: lts-13.21 + 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 @@ -# Revision history for hmatrix-sundials - -## 0.1.0.0 -- 2018-04-21 - -* 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 @@ -Copyright (c) 2018, Dominic Steinitz, Novadiscovery - -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above - copyright notice, this list of conditions and the following - disclaimer in the documentation and/or other materials provided - with the distribution. - - * Neither the name of Dominic Steinitz nor the names of other - contributors may be used to endorse or promote products derived - from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT -OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF 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 @@ -Currently only an interface to the Runge-Kutta methods: -[ARKode](https://computation.llnl.gov/projects/sundials/arkode) - -The interface is almost certainly going to change. Sundials gives a -rich set of "combinators" for controlling the solution of your problem -and reporting on how it performed. The idea is to initially mimic -hmatrix-gsl and add extra, richer functions but ultimately upgrade the -whole 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 @@ -import Distribution.Simple -main = defaultMain diff --git a/packages/sundials/diagrams/brusselator.png b/packages/sundials/diagrams/brusselator.png deleted file mode 100644 index 740cacb..0000000 Binary files a/packages/sundials/diagrams/brusselator.png and /dev/null 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 @@ -name: hmatrix-sundials -version: 0.19.0.0 -synopsis: hmatrix interface to sundials -description: An interface to the solving suite SUNDIALS. Currently, it - mimics the solving interface in hmstrix-gsl but - provides more diagnostic information and the - Butcher Tableaux (for Runge-Kutta methods). -homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials -license: BSD3 -license-file: LICENSE -author: Dominic Steinitz -maintainer: dominic@steinitz.org -copyright: Dominic Steinitz 2018, Novadiscovery 2018 -category: Math -build-type: Simple -extra-source-files: ChangeLog.md, README.md, diagrams/*.png -extra-doc-files: diagrams/*.png -cabal-version: >=1.18 - - -library - build-depends: base >=4.10 && <4.11, - inline-c >=0.6 && <0.7, - vector >=0.12 && <0.13, - template-haskell >=2.12 && <2.13, - containers >=0.5 && <0.6, - hmatrix>=0.18 - extra-libraries: sundials_arkode, - sundials_cvode - other-extensions: QuasiQuotes - hs-source-dirs: src - exposed-modules: Numeric.Sundials.ODEOpts, - Numeric.Sundials.ARKode.ODE, - Numeric.Sundials.CVode.ODE - other-modules: Numeric.Sundials.Arkode - c-sources: src/helpers.c src/helpers.h - default-language: Haskell2010 - -test-suite hmatrix-sundials-testsuite - type: exitcode-stdio-1.0 - main-is: Main.hs - other-modules: Numeric.Sundials.ODEOpts, - Numeric.Sundials.ARKode.ODE, - Numeric.Sundials.CVode.ODE, - Numeric.Sundials.Arkode - build-depends: base >=4.10 && <4.11, - inline-c >=0.6 && <0.7, - vector >=0.12 && <0.13, - template-haskell >=2.12 && <2.13, - containers >=0.5 && <0.6, - hmatrix>=0.18, - plots, - diagrams-lib, - diagrams-rasterific, - lens, - hspec - hs-source-dirs: src - extra-libraries: sundials_arkode, - sundials_cvode - c-sources: src/helpers.c src/helpers.h - 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 @@ -{-# OPTIONS_GHC -Wall #-} - -import qualified Numeric.Sundials.ARKode.ODE as ARK -import qualified Numeric.Sundials.CVode.ODE as CV -import Numeric.LinearAlgebra - -import Plots as P -import qualified Diagrams.Prelude as D -import Diagrams.Backend.Rasterific - -import Control.Lens - -import Test.Hspec - - -lorenz :: Double -> [Double] -> [Double] -lorenz _t u = [ sigma * (y - x) - , x * (rho - z) - y - , x * y - beta * z - ] - where - rho = 28.0 - sigma = 10.0 - beta = 8.0 / 3.0 - x = u !! 0 - y = u !! 1 - z = u !! 2 - -_lorenzJac :: Double -> Vector Double -> Matrix Double -_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y - , sigma , -1.0 , x - , 0.0 , (-x) , (-beta) - ] - where - rho = 28.0 - sigma = 10.0 - beta = 8.0 / 3.0 - x = u ! 0 - y = u ! 1 - z = u ! 2 - -brusselator :: Double -> [Double] -> [Double] -brusselator _t x = [ a - (w + 1) * u + v * u * u - , w * u - v * u * u - , (b - w) / eps - w * u - ] - where - a = 1.0 - b = 3.5 - eps = 5.0e-6 - u = x !! 0 - v = x !! 1 - w = x !! 2 - -_brussJac :: Double -> Vector Double -> Matrix Double -_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w) - , u * u , (-(u * u)) , 0.0 - , (-u) , u , (-1.0) / eps - u - ] - where - y = toList x - u = y !! 0 - v = y !! 1 - w = y !! 2 - eps = 5.0e-6 - -stiffish :: Double -> [Double] -> [Double] -stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] - where - lamda = -100.0 - u = v !! 0 - -stiffishV :: Double -> Vector Double -> Vector Double -stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] - where - lamda = -100.0 - u = v ! 0 - -_stiffJac :: Double -> Vector Double -> Matrix Double -_stiffJac _t _v = (1><1) [ lamda ] - where - lamda = -100.0 - -predatorPrey :: Double -> [Double] -> [Double] -predatorPrey _t v = [ x * a - b * x * y - , d * x * y - c * y - e * y * z - , (-f) * z + g * y * z - ] - where - x = v!!0 - y = v!!1 - z = v!!2 - a = 1.0 - b = 1.0 - c = 1.0 - d = 1.0 - e = 1.0 - f = 1.0 - g = 1.0 - -lSaxis :: [[Double]] -> P.Axis B D.V2 Double -lSaxis xs = P.r2Axis &~ do - let ts = xs!!0 - us = xs!!1 - vs = xs!!2 - ws = xs!!3 - P.linePlot' $ zip ts us - P.linePlot' $ zip ts vs - P.linePlot' $ zip ts ws - -kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double -kSaxis xs = P.r2Axis &~ do - P.linePlot' xs - -main :: IO () -main = do - - let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) - renderRasterific "diagrams/brusselator.png" - (D.dims2D 500.0 500.0) - (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) - - let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) - renderRasterific "diagrams/brusselatorA.png" - (D.dims2D 500.0 500.0) - (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) - - let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) - renderRasterific "diagrams/stiffish.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) - - 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]) - - 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]) - - let maxDiffA = maximum $ map abs $ - zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) - - let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) - - let maxDiffB = maximum $ map abs $ - zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0) - - let maxDiffC = maximum $ map abs $ - zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0) - - let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) - - renderRasterific "diagrams/lorenz.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1)) - - renderRasterific "diagrams/lorenz1.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2)) - - renderRasterific "diagrams/lorenz2.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2)) - - let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0]) - - renderRasterific "diagrams/predatorPrey.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1)) - - renderRasterific "diagrams/predatorPrey1.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2)) - - renderRasterific "diagrams/predatorPrey2.png" - (D.dims2D 500.0 500.0) - (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2)) - - let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0]) - - let maxDiffPpA = maximum $ map abs $ - zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0) - - hspec $ describe "Compare results" $ do - it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6 - it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6 - it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6 - it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3 - 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 @@ -{-# LANGUAGE QuasiQuotes #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE MultiWayIf #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Sundials.ARKode.ODE --- Copyright : Dominic Steinitz 2018, --- Novadiscovery 2018 --- License : BSD --- Maintainer : Dominic Steinitz --- Stability : provisional --- --- Solution of ordinary differential equation (ODE) initial value problems. --- See for more detail. --- --- A simple example: --- --- <> --- --- @ --- import Numeric.Sundials.ARKode.ODE --- import Numeric.LinearAlgebra --- --- import Plots as P --- import qualified Diagrams.Prelude as D --- import Diagrams.Backend.Rasterific --- --- brusselator :: Double -> [Double] -> [Double] --- brusselator _t x = [ a - (w + 1) * u + v * u * u --- , w * u - v * u * u --- , (b - w) / eps - w * u --- ] --- where --- a = 1.0 --- b = 3.5 --- eps = 5.0e-6 --- u = x !! 0 --- v = x !! 1 --- w = x !! 2 --- --- lSaxis :: [[Double]] -> P.Axis B D.V2 Double --- lSaxis xs = P.r2Axis &~ do --- let ts = xs!!0 --- us = xs!!1 --- vs = xs!!2 --- ws = xs!!3 --- P.linePlot' $ zip ts us --- P.linePlot' $ zip ts vs --- P.linePlot' $ zip ts ws --- --- main = do --- let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) --- renderRasterific "diagrams/brusselator.png" --- (D.dims2D 500.0 500.0) --- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) --- @ --- --- With Sundials ARKode, it is possible to retrieve the Butcher tableau for the solver. --- --- @ --- import Numeric.Sundials.ARKode.ODE --- import Numeric.LinearAlgebra --- --- import Data.List (intercalate) --- --- import Text.PrettyPrint.HughesPJClass --- --- --- butcherTableauTex :: ButcherTable -> String --- butcherTableauTex (ButcherTable m c b b2) = --- render $ --- vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}") --- , us --- , text "\\hline" --- , text bs <+> text "\\\\" --- , text b2s <+> text "\\\\" --- , text "\\end{array}" --- ] --- where --- n = rows m --- rs = toLists m --- ss = map (\r -> intercalate " & " $ map show r) rs --- ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss --- us = vcat $ map (\r -> text r <+> text "\\\\") ts --- bs = " & " ++ (intercalate " & " $ map show $ toList b) --- b2s = " & " ++ (intercalate " & " $ map show $ toList b2) --- --- main :: IO () --- main = do --- --- let res = butcherTable (SDIRK_2_1_2 undefined) --- putStrLn $ show res --- putStrLn $ butcherTableauTex res --- --- let resA = butcherTable (KVAERNO_4_2_3 undefined) --- putStrLn $ show resA --- putStrLn $ butcherTableauTex resA --- --- let resB = butcherTable (SDIRK_5_3_4 undefined) --- putStrLn $ show resB --- putStrLn $ butcherTableauTex resB --- @ --- --- Using the code above from the examples gives --- --- KVAERNO_4_2_3 --- --- \[ --- \begin{array}{c|cccc} --- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ --- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\ --- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\ --- 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\ --- \hline --- & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\ --- & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\ --- \end{array} --- \] --- --- SDIRK_2_1_2 --- --- \[ --- \begin{array}{c|cc} --- 1.0 & 1.0 & 0.0 \\ --- 0.0 & -1.0 & 1.0 \\ --- \hline --- & 0.5 & 0.5 \\ --- & 1.0 & 0.0 \\ --- \end{array} --- \] --- --- SDIRK_5_3_4 --- --- \[ --- \begin{array}{c|ccccc} --- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\ --- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\ --- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\ --- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\ --- 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\ --- \hline --- & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\ --- & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\ --- \end{array} --- \] ------------------------------------------------------------------------------ -module Numeric.Sundials.ARKode.ODE ( odeSolve - , odeSolveV - , odeSolveVWith - , odeSolveVWith' - , ButcherTable(..) - , butcherTable - , ODEMethod(..) - , StepControl(..) - ) where - -import qualified Language.C.Inline as C -import qualified Language.C.Inline.Unsafe as CU - -import Data.Monoid ((<>)) -import Data.Maybe (isJust) - -import Foreign.C.Types (CDouble, CInt, CLong) -import Foreign.Ptr (Ptr) -import Foreign.Storable (poke) - -import qualified Data.Vector.Storable as V - -import Data.Coerce (coerce) -import System.IO.Unsafe (unsafePerformIO) -import GHC.Generics (C1, Constructor, (:+:)(..), D1, Rep, Generic, M1(..), - from, conName) - -import Numeric.LinearAlgebra.Devel (createVector) - -import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows, - cols, toLists, size, reshape, - subVector, subMatrix, (><)) - -import Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..)) -import qualified Numeric.Sundials.Arkode as T -import Numeric.Sundials.Arkode (getDataFromContents, putDataInContents, arkSMax, - sDIRK_2_1_2, - bILLINGTON_3_3_2, - tRBDF2_3_3_2, - kVAERNO_4_2_3, - aRK324L2SA_DIRK_4_2_3, - cASH_5_2_4, - cASH_5_3_4, - sDIRK_5_3_4, - kVAERNO_5_3_4, - aRK436L2SA_DIRK_6_3_4, - kVAERNO_7_4_5, - aRK548L2SA_DIRK_8_4_5, - hEUN_EULER_2_1_2, - bOGACKI_SHAMPINE_4_2_3, - aRK324L2SA_ERK_4_2_3, - zONNEVELD_5_3_4, - aRK436L2SA_ERK_6_3_4, - sAYFY_ABURUB_6_3_4, - cASH_KARP_6_4_5, - fEHLBERG_6_4_5, - dORMAND_PRINCE_7_4_5, - aRK548L2SA_ERK_8_4_5, - vERNER_8_5_6, - fEHLBERG_13_7_8) - - -C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx) - -C.include "" -C.include "" -C.include "" -C.include "" -- prototypes for ARKODE fcts., consts. -C.include "" -- serial N_Vector types, fcts., macros -C.include "" -- access to dense SUNMatrix -C.include "" -- access to dense SUNLinearSolver -C.include "" -- access to ARKDls interface -C.include "" -- definition of type realtype -C.include "" -C.include "../../../helpers.h" -C.include "Numeric/Sundials/Arkode_hsc.h" - - --- | Stepping functions -data ODEMethod = SDIRK_2_1_2 Jacobian - | SDIRK_2_1_2' - | BILLINGTON_3_3_2 Jacobian - | BILLINGTON_3_3_2' - | TRBDF2_3_3_2 Jacobian - | TRBDF2_3_3_2' - | KVAERNO_4_2_3 Jacobian - | KVAERNO_4_2_3' - | ARK324L2SA_DIRK_4_2_3 Jacobian - | ARK324L2SA_DIRK_4_2_3' - | CASH_5_2_4 Jacobian - | CASH_5_2_4' - | CASH_5_3_4 Jacobian - | CASH_5_3_4' - | SDIRK_5_3_4 Jacobian - | SDIRK_5_3_4' - | KVAERNO_5_3_4 Jacobian - | KVAERNO_5_3_4' - | ARK436L2SA_DIRK_6_3_4 Jacobian - | ARK436L2SA_DIRK_6_3_4' - | KVAERNO_7_4_5 Jacobian - | KVAERNO_7_4_5' - | ARK548L2SA_DIRK_8_4_5 Jacobian - | ARK548L2SA_DIRK_8_4_5' - | HEUN_EULER_2_1_2 Jacobian - | HEUN_EULER_2_1_2' - | BOGACKI_SHAMPINE_4_2_3 Jacobian - | BOGACKI_SHAMPINE_4_2_3' - | ARK324L2SA_ERK_4_2_3 Jacobian - | ARK324L2SA_ERK_4_2_3' - | ZONNEVELD_5_3_4 Jacobian - | ZONNEVELD_5_3_4' - | ARK436L2SA_ERK_6_3_4 Jacobian - | ARK436L2SA_ERK_6_3_4' - | SAYFY_ABURUB_6_3_4 Jacobian - | SAYFY_ABURUB_6_3_4' - | CASH_KARP_6_4_5 Jacobian - | CASH_KARP_6_4_5' - | FEHLBERG_6_4_5 Jacobian - | FEHLBERG_6_4_5' - | DORMAND_PRINCE_7_4_5 Jacobian - | DORMAND_PRINCE_7_4_5' - | ARK548L2SA_ERK_8_4_5 Jacobian - | ARK548L2SA_ERK_8_4_5' - | VERNER_8_5_6 Jacobian - | VERNER_8_5_6' - | FEHLBERG_13_7_8 Jacobian - | FEHLBERG_13_7_8' - deriving Generic - -constrName :: (HasConstructor (Rep a), Generic a)=> a -> String -constrName = genericConstrName . from - -class HasConstructor (f :: * -> *) where - genericConstrName :: f x -> String - -instance HasConstructor f => HasConstructor (D1 c f) where - genericConstrName (M1 x) = genericConstrName x - -instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where - genericConstrName (L1 l) = genericConstrName l - genericConstrName (R1 r) = genericConstrName r - -instance Constructor c => HasConstructor (C1 c f) where - genericConstrName x = conName x - -instance Show ODEMethod where - show x = constrName x - --- FIXME: We can probably do better here with generics -getMethod :: ODEMethod -> Int -getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2 -getMethod (SDIRK_2_1_2') = sDIRK_2_1_2 -getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2 -getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2 -getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2 -getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2 -getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3 -getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3 -getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3 -getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3 -getMethod (CASH_5_2_4 _) = cASH_5_2_4 -getMethod (CASH_5_2_4') = cASH_5_2_4 -getMethod (CASH_5_3_4 _) = cASH_5_3_4 -getMethod (CASH_5_3_4') = cASH_5_3_4 -getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4 -getMethod (SDIRK_5_3_4') = sDIRK_5_3_4 -getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4 -getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4 -getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4 -getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4 -getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5 -getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5 -getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5 -getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5 -getMethod (HEUN_EULER_2_1_2 _) = hEUN_EULER_2_1_2 -getMethod (HEUN_EULER_2_1_2') = hEUN_EULER_2_1_2 -getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3 -getMethod (BOGACKI_SHAMPINE_4_2_3') = bOGACKI_SHAMPINE_4_2_3 -getMethod (ARK324L2SA_ERK_4_2_3 _) = aRK324L2SA_ERK_4_2_3 -getMethod (ARK324L2SA_ERK_4_2_3') = aRK324L2SA_ERK_4_2_3 -getMethod (ZONNEVELD_5_3_4 _) = zONNEVELD_5_3_4 -getMethod (ZONNEVELD_5_3_4') = zONNEVELD_5_3_4 -getMethod (ARK436L2SA_ERK_6_3_4 _) = aRK436L2SA_ERK_6_3_4 -getMethod (ARK436L2SA_ERK_6_3_4') = aRK436L2SA_ERK_6_3_4 -getMethod (SAYFY_ABURUB_6_3_4 _) = sAYFY_ABURUB_6_3_4 -getMethod (SAYFY_ABURUB_6_3_4') = sAYFY_ABURUB_6_3_4 -getMethod (CASH_KARP_6_4_5 _) = cASH_KARP_6_4_5 -getMethod (CASH_KARP_6_4_5') = cASH_KARP_6_4_5 -getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5 -getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5 -getMethod (DORMAND_PRINCE_7_4_5 _) = dORMAND_PRINCE_7_4_5 -getMethod (DORMAND_PRINCE_7_4_5') = dORMAND_PRINCE_7_4_5 -getMethod (ARK548L2SA_ERK_8_4_5 _) = aRK548L2SA_ERK_8_4_5 -getMethod (ARK548L2SA_ERK_8_4_5') = aRK548L2SA_ERK_8_4_5 -getMethod (VERNER_8_5_6 _) = vERNER_8_5_6 -getMethod (VERNER_8_5_6') = vERNER_8_5_6 -getMethod (FEHLBERG_13_7_8 _) = fEHLBERG_13_7_8 -getMethod (FEHLBERG_13_7_8') = fEHLBERG_13_7_8 - -getJacobian :: ODEMethod -> Maybe Jacobian -getJacobian (SDIRK_2_1_2 j) = Just j -getJacobian (BILLINGTON_3_3_2 j) = Just j -getJacobian (TRBDF2_3_3_2 j) = Just j -getJacobian (KVAERNO_4_2_3 j) = Just j -getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j -getJacobian (CASH_5_2_4 j) = Just j -getJacobian (CASH_5_3_4 j) = Just j -getJacobian (SDIRK_5_3_4 j) = Just j -getJacobian (KVAERNO_5_3_4 j) = Just j -getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j -getJacobian (KVAERNO_7_4_5 j) = Just j -getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j -getJacobian (HEUN_EULER_2_1_2 j) = Just j -getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j -getJacobian (ARK324L2SA_ERK_4_2_3 j) = Just j -getJacobian (ZONNEVELD_5_3_4 j) = Just j -getJacobian (ARK436L2SA_ERK_6_3_4 j) = Just j -getJacobian (SAYFY_ABURUB_6_3_4 j) = Just j -getJacobian (CASH_KARP_6_4_5 j) = Just j -getJacobian (FEHLBERG_6_4_5 j) = Just j -getJacobian (DORMAND_PRINCE_7_4_5 j) = Just j -getJacobian (ARK548L2SA_ERK_8_4_5 j) = Just j -getJacobian (VERNER_8_5_6 j) = Just j -getJacobian (FEHLBERG_13_7_8 j) = Just j -getJacobian _ = Nothing - --- | A version of 'odeSolveVWith' with reasonable default step control. -odeSolveV - :: ODEMethod - -> Maybe Double -- ^ initial step size - by default, ARKode - -- estimates the initial step size to be the - -- solution \(h\) of the equation - -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where - -- \(\ddot{y}\) is an estimated value of the - -- second derivative of the solution at \(t_0\) - -> Double -- ^ absolute tolerance for the state vector - -> Double -- ^ relative tolerance for the state vector - -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> Vector Double -- ^ initial conditions - -> Vector Double -- ^ desired solution times - -> Matrix Double -- ^ solution -odeSolveV meth hi epsAbs epsRel f y0 ts = - odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts - where - g t x0 = coerce $ f t x0 - --- | A version of 'odeSolveV' with reasonable default parameters and --- system of equations defined using lists. FIXME: we should say --- something about the fact we could use the Jacobian but don't for --- compatibility with hmatrix-gsl. -odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> [Double] -- ^ initial conditions - -> Vector Double -- ^ desired solution times - -> Matrix Double -- ^ solution -odeSolve f y0 ts = - -- FIXME: These tolerances are different from the ones in GSL - odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts) - where - g t x0 = V.fromList $ f t (V.toList x0) - -odeSolveVWith :: - ODEMethod - -> StepControl - -> Maybe Double -- ^ initial step size - by default, ARKode - -- estimates the initial step size to be the - -- solution \(h\) of the equation - -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where - -- \(\ddot{y}\) is an estimated value of the second - -- derivative of the solution at \(t_0\) - -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> V.Vector Double -- ^ Initial conditions - -> V.Vector Double -- ^ Desired solution times - -> Matrix Double -- ^ Error code or solution -odeSolveVWith method control initStepSize f y0 tt = - case odeSolveVWith' opts method control initStepSize f y0 tt of - Left (c, _v) -> error $ show c -- FIXME - Right (v, _d) -> v - where - opts = ODEOpts { maxNumSteps = 10000 - , minStep = 1.0e-12 - , relTol = error "relTol" - , absTols = error "absTol" - , initStep = error "initStep" - , maxFail = 10 - } - -odeSolveVWith' :: - ODEOpts - -> ODEMethod - -> StepControl - -> Maybe Double -- ^ initial step size - by default, ARKode - -- estimates the initial step size to be the - -- solution \(h\) of the equation - -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where - -- \(\ddot{y}\) is an estimated value of the second - -- derivative of the solution at \(t_0\) - -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> V.Vector Double -- ^ Initial conditions - -> V.Vector Double -- ^ Desired solution times - -> Either (Matrix Double, Int) (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution -odeSolveVWith' opts method control initStepSize f y0 tt = - case solveOdeC (fromIntegral $ maxFail opts) - (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) - (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) - (coerce f) (coerce y0) (coerce tt) of - Left (v, c) -> Left (reshape l (coerce v), fromIntegral c) - Right (v, d) -> Right (reshape l (coerce v), d) - where - l = size y0 - scise (X aTol rTol) = coerce (V.replicate l aTol, rTol) - scise (X' aTol rTol) = coerce (V.replicate l aTol, rTol) - scise (XX' aTol rTol yScale _yDotScale) = coerce (V.replicate l aTol, yScale * rTol) - -- FIXME; Should we check that the length of ss is correct? - scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol) - jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $ - getJacobian method - matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs } - where - nr = fromIntegral $ rows m - nc = fromIntegral $ cols m - -- FIXME: efficiency - vs = V.fromList $ map coerce $ concat $ toLists m - -solveOdeC :: - CInt -> - CLong -> - CDouble -> - CInt -> - Maybe CDouble -> - (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) -> - (V.Vector CDouble, CDouble) -> - (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> V.Vector CDouble -- ^ Initial conditions - -> V.Vector CDouble -- ^ Desired solution times - -> Either (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or - -- solution and diagnostics -solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize - jacH (aTols, rTol) fun f0 ts = unsafePerformIO $ do - - let isInitStepSize :: CInt - isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize - ss :: CDouble - ss = case initStepSize of - -- It would be better to put an error message here but - -- inline-c seems to evaluate this even if it is never - -- used :( - Nothing -> 0.0 - Just x -> x - - let dim = V.length f0 - nEq :: CLong - nEq = fromIntegral dim - nTs :: CInt - nTs = fromIntegral $ V.length ts - quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs)) - qMatMut <- V.thaw quasiMatrixRes - diagnostics :: V.Vector CLong <- createVector 10 -- FIXME - diagMut <- V.thaw diagnostics - -- We need the types that sundials expects. These are tied together - -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty! - let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt - funIO x y f _ptr = do - -- Convert the pointer we get from C (y) to a vector, and then - -- apply the user-supplied function. - fImm <- fun x <$> getDataFromContents dim y - -- Fill in the provided pointer with the resulting vector. - putDataInContents fImm dim f - -- FIXME: I don't understand what this comment means - -- Unsafe since the function will be called many times. - [CU.exp| int{ 0 } |] - let isJac :: CInt - isJac = fromIntegral $ fromEnum $ isJust jacH - jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix -> - Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector -> - IO CInt - jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do - case jacH of - Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined" - Just jacI -> do j <- jacI t <$> getDataFromContents dim y - poke jacS j - -- FIXME: I don't understand what this comment means - -- Unsafe since the function will be called many times. - [CU.exp| int{ 0 } |] - - res <- [C.block| int { - /* general problem variables */ - - int flag; /* reusable error-checking flag */ - int i, j; /* reusable loop indices */ - N_Vector y = NULL; /* empty vector for storing solution */ - N_Vector tv = NULL; /* empty vector for storing absolute tolerances */ - SUNMatrix A = NULL; /* empty matrix for linear solver */ - SUNLinearSolver LS = NULL; /* empty linear solver object */ - void *arkode_mem = NULL; /* empty ARKode memory structure */ - realtype t; - long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf; - - /* general problem parameters */ - - realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */ - sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */ - - /* Initialize data structures */ - - y = N_VNew_Serial(NEQ); /* Create serial vector for solution */ - if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1; - /* Specify initial condition */ - for (i = 0; i < NEQ; i++) { - NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i]; - }; - - tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */ - if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1; - /* Specify tolerances */ - for (i = 0; i < NEQ; i++) { - NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i]; - }; - - arkode_mem = ARKodeCreate(); /* Create the solver memory */ - if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1; - - /* Call ARKodeInit to initialize the integrator memory and specify the */ - /* right-hand side function in y'=f(t,y), the inital time T0, and */ - /* the initial dependent variable vector y. Note: we treat the */ - /* problem as fully implicit and set f_E to NULL and f_I to f. */ - - /* Here we use the C types defined in helpers.h which tie up with */ - /* the Haskell types defined in CLangToHaskellTypes */ - if ($(int method) < MIN_DIRK_NUM) { - flag = ARKodeInit(arkode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y); - if (check_flag(&flag, "ARKodeInit", 1)) return 1; - } else { - flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y); - if (check_flag(&flag, "ARKodeInit", 1)) return 1; - } - - flag = ARKodeSetMinStep(arkode_mem, $(double minStep_)); - if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1; - flag = ARKodeSetMaxNumSteps(arkode_mem, $(long int maxNumSteps_)); - if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1; - flag = ARKodeSetMaxErrTestFails(arkode_mem, $(int maxErrTestFails)); - if (check_flag(&flag, "ARKodeSetMaxErrTestFails", 1)) return 1; - - /* Set routines */ - flag = ARKodeSVtolerances(arkode_mem, $(double rTol), tv); - if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1; - - /* Initialize dense matrix data structure and solver */ - A = SUNDenseMatrix(NEQ, NEQ); - if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1; - LS = SUNDenseLinearSolver(y, A); - if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1; - - /* Attach matrix and linear solver */ - flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); - if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1; - - /* Set the initial step size if there is one */ - if ($(int isInitStepSize)) { - /* FIXME: We could check if the initial step size is 0 */ - /* or even NaN and then throw an error */ - flag = ARKodeSetInitStep(arkode_mem, $(double ss)); - if (check_flag(&flag, "ARKodeSetInitStep", 1)) return 1; - } - - /* Set the Jacobian if there is one */ - if ($(int isJac)) { - flag = ARKDlsSetJacFn(arkode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[]))); - if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1; - } - - /* Store initial conditions */ - for (j = 0; j < NEQ; j++) { - ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j); - } - - /* Explicitly set the method */ - if ($(int method) >= MIN_DIRK_NUM) { - flag = ARKodeSetIRKTableNum(arkode_mem, $(int method)); - if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1; - } else { - flag = ARKodeSetERKTableNum(arkode_mem, $(int method)); - if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1; - } - - /* Main time-stepping loop: calls ARKode to perform the integration */ - /* Stops when the final time has been reached */ - for (i = 1; i < $(int nTs); i++) { - - flag = ARKode(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */ - if (check_flag(&flag, "ARKode solver failure, stopping integration", 1)) return 1; - - /* Store the results for Haskell */ - for (j = 0; j < NEQ; j++) { - ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j); - } - } - - /* Get some final statistics on how the solve progressed */ - - flag = ARKodeGetNumSteps(arkode_mem, &nst); - check_flag(&flag, "ARKodeGetNumSteps", 1); - ($vec-ptr:(long int *diagMut))[0] = nst; - - flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a); - check_flag(&flag, "ARKodeGetNumStepAttempts", 1); - ($vec-ptr:(long int *diagMut))[1] = nst_a; - - flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi); - check_flag(&flag, "ARKodeGetNumRhsEvals", 1); - ($vec-ptr:(long int *diagMut))[2] = nfe; - ($vec-ptr:(long int *diagMut))[3] = nfi; - - flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups); - check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1); - ($vec-ptr:(long int *diagMut))[4] = nsetups; - - flag = ARKodeGetNumErrTestFails(arkode_mem, &netf); - check_flag(&flag, "ARKodeGetNumErrTestFails", 1); - ($vec-ptr:(long int *diagMut))[5] = netf; - - flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni); - check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1); - ($vec-ptr:(long int *diagMut))[6] = nni; - - flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn); - check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1); - ($vec-ptr:(long int *diagMut))[7] = ncfn; - - flag = ARKDlsGetNumJacEvals(arkode_mem, &nje); - check_flag(&flag, "ARKDlsGetNumJacEvals", 1); - ($vec-ptr:(long int *diagMut))[8] = ncfn; - - flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS); - check_flag(&flag, "ARKDlsGetNumRhsEvals", 1); - ($vec-ptr:(long int *diagMut))[9] = ncfn; - - /* Clean up and return */ - N_VDestroy(y); /* Free y vector */ - N_VDestroy(tv); /* Free tv vector */ - ARKodeFree(&arkode_mem); /* Free integrator memory */ - SUNLinSolFree(LS); /* Free linear solver */ - SUNMatDestroy(A); /* Free A matrix */ - - return flag; - } |] - preD <- V.freeze diagMut - let d = SundialsDiagnostics (fromIntegral $ preD V.!0) - (fromIntegral $ preD V.!1) - (fromIntegral $ preD V.!2) - (fromIntegral $ preD V.!3) - (fromIntegral $ preD V.!4) - (fromIntegral $ preD V.!5) - (fromIntegral $ preD V.!6) - (fromIntegral $ preD V.!7) - (fromIntegral $ preD V.!8) - (fromIntegral $ preD V.!9) - m <- V.freeze qMatMut - if res == 0 - then do - return $ Right (m, d) - else do - return $ Left (m, res) - -data ButcherTable = ButcherTable { am :: Matrix Double - , cv :: Vector Double - , bv :: Vector Double - , b2v :: Vector Double - } - deriving Show - -data ButcherTable' a = ButcherTable' { am' :: V.Vector a - , cv' :: V.Vector a - , bv' :: V.Vector a - , b2v' :: V.Vector a - } - deriving Show - -butcherTable :: ODEMethod -> ButcherTable -butcherTable method = - case getBT method of - Left c -> error $ show c -- FIXME - Right (ButcherTable' v w x y, sqp) -> - ButcherTable { am = subMatrix (0, 0) (s, s) $ (arkSMax >< arkSMax) (V.toList v) - , cv = subVector 0 s w - , bv = subVector 0 s x - , b2v = subVector 0 s y - } - where - s = fromIntegral $ sqp V.! 0 - -getBT :: ODEMethod -> Either Int (ButcherTable' Double, V.Vector Int) -getBT method = case getButcherTable method of - Left c -> - Left $ fromIntegral c - Right (ButcherTable' a b c d, sqp) -> - Right $ ( ButcherTable' (coerce a) (coerce b) (coerce c) (coerce d) - , V.map fromIntegral sqp ) - -getButcherTable :: ODEMethod - -> Either CInt (ButcherTable' CDouble, V.Vector CInt) -getButcherTable method = unsafePerformIO $ do - -- ARKode seems to want an ODE in order to set and then get the - -- Butcher tableau so here's one to keep it happy - let funI :: CDouble -> V.Vector CDouble -> V.Vector CDouble - funI _t ys = V.fromList [ ys V.! 0 ] - let funE :: CDouble -> V.Vector CDouble -> V.Vector CDouble - funE _t ys = V.fromList [ ys V.! 0 ] - f0 = V.fromList [ 1.0 ] - ts = V.fromList [ 0.0 ] - dim = V.length f0 - nEq :: CLong - nEq = fromIntegral dim - mN :: CInt - mN = fromIntegral $ getMethod method - - btSQP :: V.Vector CInt <- createVector 3 - btSQPMut <- V.thaw btSQP - btAs :: V.Vector CDouble <- createVector (arkSMax * arkSMax) - btAsMut <- V.thaw btAs - btCs :: V.Vector CDouble <- createVector arkSMax - btBs :: V.Vector CDouble <- createVector arkSMax - btB2s :: V.Vector CDouble <- createVector arkSMax - btCsMut <- V.thaw btCs - btBsMut <- V.thaw btBs - btB2sMut <- V.thaw btB2s - let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt - funIOI x y f _ptr = do - fImm <- funI x <$> getDataFromContents dim y - putDataInContents fImm dim f - -- FIXME: I don't understand what this comment means - -- Unsafe since the function will be called many times. - [CU.exp| int{ 0 } |] - let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt - funIOE x y f _ptr = do - fImm <- funE x <$> getDataFromContents dim y - putDataInContents fImm dim f - -- FIXME: I don't understand what this comment means - -- Unsafe since the function will be called many times. - [CU.exp| int{ 0 } |] - res <- [C.block| int { - /* general problem variables */ - - int flag; /* reusable error-checking flag */ - N_Vector y = NULL; /* empty vector for storing solution */ - void *arkode_mem = NULL; /* empty ARKode memory structure */ - int i, j; /* reusable loop indices */ - - /* general problem parameters */ - - realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */ - sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars */ - - /* Initialize data structures */ - - y = N_VNew_Serial(NEQ); /* Create serial vector for solution */ - if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1; - /* Specify initial condition */ - for (i = 0; i < NEQ; i++) { - NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i]; - }; - arkode_mem = ARKodeCreate(); /* Create the solver memory */ - if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1; - - 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); - if (check_flag(&flag, "ARKodeInit", 1)) return 1; - - if ($(int mN) >= MIN_DIRK_NUM) { - flag = ARKodeSetIRKTableNum(arkode_mem, $(int mN)); - if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1; - } else { - flag = ARKodeSetERKTableNum(arkode_mem, $(int mN)); - if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1; - } - - int s, q, p; - realtype *ai = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype)); - realtype *ae = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype)); - realtype *ci = (realtype *)malloc(ARK_S_MAX * sizeof(realtype)); - realtype *ce = (realtype *)malloc(ARK_S_MAX * sizeof(realtype)); - realtype *bi = (realtype *)malloc(ARK_S_MAX * sizeof(realtype)); - realtype *be = (realtype *)malloc(ARK_S_MAX * sizeof(realtype)); - realtype *b2i = (realtype *)malloc(ARK_S_MAX * sizeof(realtype)); - realtype *b2e = (realtype *)malloc(ARK_S_MAX * sizeof(realtype)); - flag = ARKodeGetCurrentButcherTables(arkode_mem, &s, &q, &p, ai, ae, ci, ce, bi, be, b2i, b2e); - if (check_flag(&flag, "ARKode", 1)) return 1; - $vec-ptr:(int *btSQPMut)[0] = s; - $vec-ptr:(int *btSQPMut)[1] = q; - $vec-ptr:(int *btSQPMut)[2] = p; - for (i = 0; i < s; i++) { - for (j = 0; j < s; j++) { - /* FIXME: double should be realtype */ - ($vec-ptr:(double *btAsMut))[i * ARK_S_MAX + j] = ai[i * ARK_S_MAX + j]; - } - } - - for (i = 0; i < s; i++) { - ($vec-ptr:(double *btCsMut))[i] = ci[i]; - ($vec-ptr:(double *btBsMut))[i] = bi[i]; - ($vec-ptr:(double *btB2sMut))[i] = b2i[i]; - } - - /* Clean up and return */ - N_VDestroy(y); /* Free y vector */ - ARKodeFree(&arkode_mem); /* Free integrator memory */ - - return flag; - } |] - if res == 0 - then do - x <- V.freeze btAsMut - y <- V.freeze btSQPMut - z <- V.freeze btCsMut - u <- V.freeze btBsMut - v <- V.freeze btB2sMut - return $ Right (ButcherTable' { am' = x, cv' = z, bv' = u, b2v' = v }, y) - else do - return $ Left res - --- | Adaptive step-size control --- functions. --- --- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control) --- allows the user to control the step size adjustment using --- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where --- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\) --- is the required relative error, \(s_i\) is a vector of scaling --- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and --- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\). --- --- [ARKode](https://computation.llnl.gov/projects/sundials/arkode) --- allows the user to control the step size adjustment using --- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with --- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl), --- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no --- effect. -data 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 - | 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 - | XX' Double Double Double Double -- ^ include both via relative tolerance - -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\) - | 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 @@ -{-# LANGUAGE QuasiQuotes #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE EmptyDataDecls #-} - -module Numeric.Sundials.Arkode where - -import Foreign -import Foreign.C.Types - -import Language.C.Types as CT - -import qualified Data.Vector.Storable as VS -import qualified Data.Vector.Storable.Mutable as VM - -import qualified Language.Haskell.TH as TH -import qualified Data.Map as Map -import Language.C.Inline.Context - -import qualified Data.Vector.Storable as V - - -#include -#include -#include -#include -#include -#include -#include - - -data SunVector -data SunMatrix = SunMatrix { rows :: CInt - , cols :: CInt - , vals :: V.Vector CDouble - } - --- | This is true only if configured/ built as 64 bits -type SunIndexType = CLong - -sunTypesTable :: Map.Map TypeSpecifier TH.TypeQ -sunTypesTable = Map.fromList - [ - (TypeName "sunindextype", [t| SunIndexType |] ) - , (TypeName "SunVector", [t| SunVector |] ) - , (TypeName "SunMatrix", [t| SunMatrix |] ) - ] - -sunCtx :: Context -sunCtx = mempty {ctxTypesTable = sunTypesTable} - -getMatrixDataFromContents :: Ptr SunMatrix -> IO SunMatrix -getMatrixDataFromContents ptr = do - qtr <- getContentMatrixPtr ptr - rs <- getNRows qtr - cs <- getNCols qtr - rtr <- getMatrixData qtr - vs <- vectorFromC (fromIntegral $ rs * cs) rtr - return $ SunMatrix { rows = rs, cols = cs, vals = vs } - -putMatrixDataFromContents :: SunMatrix -> Ptr SunMatrix -> IO () -putMatrixDataFromContents mat ptr = do - let rs = rows mat - cs = cols mat - vs = vals mat - qtr <- getContentMatrixPtr ptr - putNRows rs qtr - putNCols cs qtr - rtr <- getMatrixData qtr - vectorToC vs (fromIntegral $ rs * cs) rtr - -instance Storable SunMatrix where - poke = flip putMatrixDataFromContents - peek = getMatrixDataFromContents - sizeOf _ = error "sizeOf not supported for SunMatrix" - alignment _ = error "alignment not supported for SunMatrix" - -vectorFromC :: Storable a => Int -> Ptr a -> IO (VS.Vector a) -vectorFromC len ptr = do - ptr' <- newForeignPtr_ ptr - VS.freeze $ VM.unsafeFromForeignPtr0 ptr' len - -vectorToC :: Storable a => VS.Vector a -> Int -> Ptr a -> IO () -vectorToC vec len ptr = do - ptr' <- newForeignPtr_ ptr - VS.copy (VM.unsafeFromForeignPtr0 ptr' len) vec - -getDataFromContents :: Int -> Ptr SunVector -> IO (VS.Vector CDouble) -getDataFromContents len ptr = do - qtr <- getContentPtr ptr - rtr <- getData qtr - vectorFromC len rtr - -putDataInContents :: Storable a => VS.Vector a -> Int -> Ptr b -> IO () -putDataInContents vec len ptr = do - qtr <- getContentPtr ptr - rtr <- getData qtr - vectorToC vec len rtr - -#def typedef struct _generic_N_Vector SunVector; -#def typedef struct _N_VectorContent_Serial SunContent; - -#def typedef struct _generic_SUNMatrix SunMatrix; -#def typedef struct _SUNMatrixContent_Dense SunMatrixContent; - -getContentMatrixPtr :: Storable a => Ptr b -> IO a -getContentMatrixPtr ptr = (#peek SunMatrix, content) ptr - -getNRows :: Ptr b -> IO CInt -getNRows ptr = (#peek SunMatrixContent, M) ptr -putNRows :: CInt -> Ptr b -> IO () -putNRows nr ptr = (#poke SunMatrixContent, M) ptr nr - -getNCols :: Ptr b -> IO CInt -getNCols ptr = (#peek SunMatrixContent, N) ptr -putNCols :: CInt -> Ptr b -> IO () -putNCols nc ptr = (#poke SunMatrixContent, N) ptr nc - -getMatrixData :: Storable a => Ptr b -> IO a -getMatrixData ptr = (#peek SunMatrixContent, data) ptr - -getContentPtr :: Storable a => Ptr b -> IO a -getContentPtr ptr = (#peek SunVector, content) ptr - -getData :: Storable a => Ptr b -> IO a -getData ptr = (#peek SunContent, data) ptr - -cV_ADAMS :: Int -cV_ADAMS = #const CV_ADAMS -cV_BDF :: Int -cV_BDF = #const CV_BDF - -arkSMax :: Int -arkSMax = #const ARK_S_MAX - -mIN_DIRK_NUM, mAX_DIRK_NUM :: Int -mIN_DIRK_NUM = #const MIN_DIRK_NUM -mAX_DIRK_NUM = #const MAX_DIRK_NUM - --- FIXME: We could just use inline-c instead - --- Butcher table accessors -- implicit -sDIRK_2_1_2 :: Int -sDIRK_2_1_2 = #const SDIRK_2_1_2 -bILLINGTON_3_3_2 :: Int -bILLINGTON_3_3_2 = #const BILLINGTON_3_3_2 -tRBDF2_3_3_2 :: Int -tRBDF2_3_3_2 = #const TRBDF2_3_3_2 -kVAERNO_4_2_3 :: Int -kVAERNO_4_2_3 = #const KVAERNO_4_2_3 -aRK324L2SA_DIRK_4_2_3 :: Int -aRK324L2SA_DIRK_4_2_3 = #const ARK324L2SA_DIRK_4_2_3 -cASH_5_2_4 :: Int -cASH_5_2_4 = #const CASH_5_2_4 -cASH_5_3_4 :: Int -cASH_5_3_4 = #const CASH_5_3_4 -sDIRK_5_3_4 :: Int -sDIRK_5_3_4 = #const SDIRK_5_3_4 -kVAERNO_5_3_4 :: Int -kVAERNO_5_3_4 = #const KVAERNO_5_3_4 -aRK436L2SA_DIRK_6_3_4 :: Int -aRK436L2SA_DIRK_6_3_4 = #const ARK436L2SA_DIRK_6_3_4 -kVAERNO_7_4_5 :: Int -kVAERNO_7_4_5 = #const KVAERNO_7_4_5 -aRK548L2SA_DIRK_8_4_5 :: Int -aRK548L2SA_DIRK_8_4_5 = #const ARK548L2SA_DIRK_8_4_5 - --- #define DEFAULT_DIRK_2 SDIRK_2_1_2 --- #define DEFAULT_DIRK_3 ARK324L2SA_DIRK_4_2_3 --- #define DEFAULT_DIRK_4 SDIRK_5_3_4 --- #define DEFAULT_DIRK_5 ARK548L2SA_DIRK_8_4_5 - --- Butcher table accessors -- explicit -hEUN_EULER_2_1_2 :: Int -hEUN_EULER_2_1_2 = #const HEUN_EULER_2_1_2 -bOGACKI_SHAMPINE_4_2_3 :: Int -bOGACKI_SHAMPINE_4_2_3 = #const BOGACKI_SHAMPINE_4_2_3 -aRK324L2SA_ERK_4_2_3 :: Int -aRK324L2SA_ERK_4_2_3 = #const ARK324L2SA_ERK_4_2_3 -zONNEVELD_5_3_4 :: Int -zONNEVELD_5_3_4 = #const ZONNEVELD_5_3_4 -aRK436L2SA_ERK_6_3_4 :: Int -aRK436L2SA_ERK_6_3_4 = #const ARK436L2SA_ERK_6_3_4 -sAYFY_ABURUB_6_3_4 :: Int -sAYFY_ABURUB_6_3_4 = #const SAYFY_ABURUB_6_3_4 -cASH_KARP_6_4_5 :: Int -cASH_KARP_6_4_5 = #const CASH_KARP_6_4_5 -fEHLBERG_6_4_5 :: Int -fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5 -dORMAND_PRINCE_7_4_5 :: Int -dORMAND_PRINCE_7_4_5 = #const DORMAND_PRINCE_7_4_5 -aRK548L2SA_ERK_8_4_5 :: Int -aRK548L2SA_ERK_8_4_5 = #const ARK548L2SA_ERK_8_4_5 -vERNER_8_5_6 :: Int -vERNER_8_5_6 = #const VERNER_8_5_6 -fEHLBERG_13_7_8 :: Int -fEHLBERG_13_7_8 = #const FEHLBERG_13_7_8 - --- #define DEFAULT_ERK_2 HEUN_EULER_2_1_2 --- #define DEFAULT_ERK_3 BOGACKI_SHAMPINE_4_2_3 --- #define DEFAULT_ERK_4 ZONNEVELD_5_3_4 --- #define DEFAULT_ERK_5 CASH_KARP_6_4_5 --- #define DEFAULT_ERK_6 VERNER_8_5_6 --- #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 @@ -{-# OPTIONS_GHC -Wall #-} - -{-# LANGUAGE QuasiQuotes #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE MultiWayIf #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE ScopedTypeVariables #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Sundials.CVode.ODE --- Copyright : Dominic Steinitz 2018, --- Novadiscovery 2018 --- License : BSD --- Maintainer : Dominic Steinitz --- Stability : provisional --- --- Solution of ordinary differential equation (ODE) initial value problems. --- --- --- --- A simple example: --- --- <> --- --- @ --- import Numeric.Sundials.CVode.ODE --- import Numeric.LinearAlgebra --- --- import Plots as P --- import qualified Diagrams.Prelude as D --- import Diagrams.Backend.Rasterific --- --- brusselator :: Double -> [Double] -> [Double] --- brusselator _t x = [ a - (w + 1) * u + v * u * u --- , w * u - v * u * u --- , (b - w) / eps - w * u --- ] --- where --- a = 1.0 --- b = 3.5 --- eps = 5.0e-6 --- u = x !! 0 --- v = x !! 1 --- w = x !! 2 --- --- lSaxis :: [[Double]] -> P.Axis B D.V2 Double --- lSaxis xs = P.r2Axis &~ do --- let ts = xs!!0 --- us = xs!!1 --- vs = xs!!2 --- ws = xs!!3 --- P.linePlot' $ zip ts us --- P.linePlot' $ zip ts vs --- P.linePlot' $ zip ts ws --- --- main = do --- let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) --- renderRasterific "diagrams/brusselator.png" --- (D.dims2D 500.0 500.0) --- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) --- @ --- ------------------------------------------------------------------------------ -module Numeric.Sundials.CVode.ODE ( odeSolve - , odeSolveV - , odeSolveVWith - , odeSolveVWith' - , ODEMethod(..) - , StepControl(..) - ) where - -import qualified Language.C.Inline as C -import qualified Language.C.Inline.Unsafe as CU - -import Data.Monoid ((<>)) -import Data.Maybe (isJust) - -import Foreign.C.Types (CDouble, CInt, CLong) -import Foreign.Ptr (Ptr) -import Foreign.Storable (poke) - -import qualified Data.Vector.Storable as V - -import Data.Coerce (coerce) -import System.IO.Unsafe (unsafePerformIO) - -import Numeric.LinearAlgebra.Devel (createVector) - -import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows, - cols, toLists, size, reshape) - -import Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF, - getDataFromContents, putDataInContents) -import qualified Numeric.Sundials.Arkode as T -import Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..)) - - -C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx) - -C.include "" -C.include "" -C.include "" -C.include "" -- prototypes for CVODE fcts., consts. -C.include "" -- serial N_Vector types, fcts., macros -C.include "" -- access to dense SUNMatrix -C.include "" -- access to dense SUNLinearSolver -C.include "" -- access to CVDls interface -C.include "" -- definition of type realtype -C.include "" -C.include "../../../helpers.h" -C.include "Numeric/Sundials/Arkode_hsc.h" - - --- | Stepping functions -data ODEMethod = ADAMS - | BDF - -getMethod :: ODEMethod -> Int -getMethod (ADAMS) = cV_ADAMS -getMethod (BDF) = cV_BDF - -getJacobian :: ODEMethod -> Maybe Jacobian -getJacobian _ = Nothing - --- | A version of 'odeSolveVWith' with reasonable default step control. -odeSolveV - :: ODEMethod - -> Maybe Double -- ^ initial step size - by default, CVode - -- estimates the initial step size to be the - -- solution \(h\) of the equation - -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where - -- \(\ddot{y}\) is an estimated value of the - -- second derivative of the solution at \(t_0\) - -> Double -- ^ absolute tolerance for the state vector - -> Double -- ^ relative tolerance for the state vector - -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> Vector Double -- ^ initial conditions - -> Vector Double -- ^ desired solution times - -> Matrix Double -- ^ solution -odeSolveV meth hi epsAbs epsRel f y0 ts = - odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts - where - g t x0 = coerce $ f t x0 - --- | A version of 'odeSolveV' with reasonable default parameters and --- system of equations defined using lists. FIXME: we should say --- something about the fact we could use the Jacobian but don't for --- compatibility with hmatrix-gsl. -odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> [Double] -- ^ initial conditions - -> Vector Double -- ^ desired solution times - -> Matrix Double -- ^ solution -odeSolve f y0 ts = - -- FIXME: These tolerances are different from the ones in GSL - odeSolveVWith BDF (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts) - where - g t x0 = V.fromList $ f t (V.toList x0) - -odeSolveVWith :: - ODEMethod - -> StepControl - -> Maybe Double -- ^ initial step size - by default, CVode - -- estimates the initial step size to be the - -- solution \(h\) of the equation - -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where - -- \(\ddot{y}\) is an estimated value of the second - -- derivative of the solution at \(t_0\) - -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> V.Vector Double -- ^ Initial conditions - -> V.Vector Double -- ^ Desired solution times - -> Matrix Double -- ^ Error code or solution -odeSolveVWith method control initStepSize f y0 tt = - case odeSolveVWith' opts method control initStepSize f y0 tt of - Left (c, _v) -> error $ show c -- FIXME - Right (v, _d) -> v - where - opts = ODEOpts { maxNumSteps = 10000 - , minStep = 1.0e-12 - , relTol = error "relTol" - , absTols = error "absTol" - , initStep = error "initStep" - , maxFail = 10 - } - -odeSolveVWith' :: - ODEOpts - -> ODEMethod - -> StepControl - -> Maybe Double -- ^ initial step size - by default, CVode - -- estimates the initial step size to be the - -- solution \(h\) of the equation - -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where - -- \(\ddot{y}\) is an estimated value of the second - -- derivative of the solution at \(t_0\) - -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> V.Vector Double -- ^ Initial conditions - -> V.Vector Double -- ^ Desired solution times - -> Either (Matrix Double, Int) (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution -odeSolveVWith' opts method control initStepSize f y0 tt = - case solveOdeC (fromIntegral $ maxFail opts) - (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) - (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) - (coerce f) (coerce y0) (coerce tt) of - Left (v, c) -> Left (reshape l (coerce v), fromIntegral c) - Right (v, d) -> Right (reshape l (coerce v), d) - where - l = size y0 - scise (X aTol rTol) = coerce (V.replicate l aTol, rTol) - scise (X' aTol rTol) = coerce (V.replicate l aTol, rTol) - scise (XX' aTol rTol yScale _yDotScale) = coerce (V.replicate l aTol, yScale * rTol) - -- FIXME; Should we check that the length of ss is correct? - scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol) - jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $ - getJacobian method - matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs } - where - nr = fromIntegral $ rows m - nc = fromIntegral $ cols m - -- FIXME: efficiency - vs = V.fromList $ map coerce $ concat $ toLists m - -solveOdeC :: - CInt -> - CLong -> - CDouble -> - CInt -> - Maybe CDouble -> - (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) -> - (V.Vector CDouble, CDouble) -> - (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) - -> V.Vector CDouble -- ^ Initial conditions - -> V.Vector CDouble -- ^ Desired solution times - -> Either (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or - -- solution and diagnostics -solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize - jacH (aTols, rTol) fun f0 ts = - unsafePerformIO $ do - - let isInitStepSize :: CInt - isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize - ss :: CDouble - ss = case initStepSize of - -- It would be better to put an error message here but - -- inline-c seems to evaluate this even if it is never - -- used :( - Nothing -> 0.0 - Just x -> x - - let dim = V.length f0 - nEq :: CLong - nEq = fromIntegral dim - nTs :: CInt - nTs = fromIntegral $ V.length ts - quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs)) - qMatMut <- V.thaw quasiMatrixRes - diagnostics :: V.Vector CLong <- createVector 10 -- FIXME - diagMut <- V.thaw diagnostics - -- We need the types that sundials expects. These are tied together - -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty! - let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt - funIO x y f _ptr = do - -- Convert the pointer we get from C (y) to a vector, and then - -- apply the user-supplied function. - fImm <- fun x <$> getDataFromContents dim y - -- Fill in the provided pointer with the resulting vector. - putDataInContents fImm dim f - -- FIXME: I don't understand what this comment means - -- Unsafe since the function will be called many times. - [CU.exp| int{ 0 } |] - let isJac :: CInt - isJac = fromIntegral $ fromEnum $ isJust jacH - jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix -> - Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector -> - IO CInt - jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do - case jacH of - Nothing -> error "Numeric.Sundials.CVode.ODE: Jacobian not defined" - Just jacI -> do j <- jacI t <$> getDataFromContents dim y - poke jacS j - -- FIXME: I don't understand what this comment means - -- Unsafe since the function will be called many times. - [CU.exp| int{ 0 } |] - - res <- [C.block| int { - /* general problem variables */ - - int flag; /* reusable error-checking flag */ - int i, j; /* reusable loop indices */ - N_Vector y = NULL; /* empty vector for storing solution */ - N_Vector tv = NULL; /* empty vector for storing absolute tolerances */ - - SUNMatrix A = NULL; /* empty matrix for linear solver */ - SUNLinearSolver LS = NULL; /* empty linear solver object */ - void *cvode_mem = NULL; /* empty CVODE memory structure */ - realtype t; - long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge; - - /* general problem parameters */ - - realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */ - sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */ - - /* Initialize data structures */ - - y = N_VNew_Serial(NEQ); /* Create serial vector for solution */ - if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1; - /* Specify initial condition */ - for (i = 0; i < NEQ; i++) { - NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i]; - }; - - cvode_mem = CVodeCreate($(int method), CV_NEWTON); - if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1); - - /* Call CVodeInit to initialize the integrator memory and specify the - * user's right hand side function in y'=f(t,y), the inital time T0, and - * the initial dependent variable vector y. */ - flag = CVodeInit(cvode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y); - if (check_flag(&flag, "CVodeInit", 1)) return(1); - - tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */ - if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1; - /* Specify tolerances */ - for (i = 0; i < NEQ; i++) { - NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i]; - }; - - flag = CVodeSetMinStep(cvode_mem, $(double minStep_)); - if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1; - flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_)); - if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1; - flag = CVodeSetMaxErrTestFails(cvode_mem, $(int maxErrTestFails)); - if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return 1; - - /* Call CVodeSVtolerances to specify the scalar relative tolerance - * and vector absolute tolerances */ - flag = CVodeSVtolerances(cvode_mem, $(double rTol), tv); - if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1); - - /* Initialize dense matrix data structure and solver */ - A = SUNDenseMatrix(NEQ, NEQ); - if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1; - LS = SUNDenseLinearSolver(y, A); - if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1; - - /* Attach matrix and linear solver */ - flag = CVDlsSetLinearSolver(cvode_mem, LS, A); - if (check_flag(&flag, "CVDlsSetLinearSolver", 1)) return 1; - - /* Set the initial step size if there is one */ - if ($(int isInitStepSize)) { - /* FIXME: We could check if the initial step size is 0 */ - /* or even NaN and then throw an error */ - flag = CVodeSetInitStep(cvode_mem, $(double ss)); - if (check_flag(&flag, "CVodeSetInitStep", 1)) return 1; - } - - /* Set the Jacobian if there is one */ - if ($(int isJac)) { - flag = CVDlsSetJacFn(cvode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[]))); - if (check_flag(&flag, "CVDlsSetJacFn", 1)) return 1; - } - - /* Store initial conditions */ - for (j = 0; j < NEQ; j++) { - ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j); - } - - /* Main time-stepping loop: calls CVode to perform the integration */ - /* Stops when the final time has been reached */ - for (i = 1; i < $(int nTs); i++) { - - flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[i], y, &t, CV_NORMAL); /* call integrator */ - if (check_flag(&flag, "CVode solver failure, stopping integration", 1)) return 1; - - /* Store the results for Haskell */ - for (j = 0; j < NEQ; j++) { - ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j); - } - } - - /* Get some final statistics on how the solve progressed */ - - flag = CVodeGetNumSteps(cvode_mem, &nst); - check_flag(&flag, "CVodeGetNumSteps", 1); - ($vec-ptr:(long int *diagMut))[0] = nst; - - /* FIXME */ - ($vec-ptr:(long int *diagMut))[1] = 0; - - flag = CVodeGetNumRhsEvals(cvode_mem, &nfe); - check_flag(&flag, "CVodeGetNumRhsEvals", 1); - ($vec-ptr:(long int *diagMut))[2] = nfe; - /* FIXME */ - ($vec-ptr:(long int *diagMut))[3] = 0; - - flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups); - check_flag(&flag, "CVodeGetNumLinSolvSetups", 1); - ($vec-ptr:(long int *diagMut))[4] = nsetups; - - flag = CVodeGetNumErrTestFails(cvode_mem, &netf); - check_flag(&flag, "CVodeGetNumErrTestFails", 1); - ($vec-ptr:(long int *diagMut))[5] = netf; - - flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni); - check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1); - ($vec-ptr:(long int *diagMut))[6] = nni; - - flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn); - check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1); - ($vec-ptr:(long int *diagMut))[7] = ncfn; - - flag = CVDlsGetNumJacEvals(cvode_mem, &nje); - check_flag(&flag, "CVDlsGetNumJacEvals", 1); - ($vec-ptr:(long int *diagMut))[8] = ncfn; - - flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS); - check_flag(&flag, "CVDlsGetNumRhsEvals", 1); - ($vec-ptr:(long int *diagMut))[9] = ncfn; - - /* Clean up and return */ - - N_VDestroy(y); /* Free y vector */ - N_VDestroy(tv); /* Free tv vector */ - CVodeFree(&cvode_mem); /* Free integrator memory */ - SUNLinSolFree(LS); /* Free linear solver */ - SUNMatDestroy(A); /* Free A matrix */ - - return flag; - } |] - preD <- V.freeze diagMut - let d = SundialsDiagnostics (fromIntegral $ preD V.!0) - (fromIntegral $ preD V.!1) - (fromIntegral $ preD V.!2) - (fromIntegral $ preD V.!3) - (fromIntegral $ preD V.!4) - (fromIntegral $ preD V.!5) - (fromIntegral $ preD V.!6) - (fromIntegral $ preD V.!7) - (fromIntegral $ preD V.!8) - (fromIntegral $ preD V.!9) - m <- V.freeze qMatMut - if res == 0 - then do - return $ Right (m, d) - else do - return $ Left (m, res) - --- | Adaptive step-size control --- functions. --- --- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control) --- allows the user to control the step size adjustment using --- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where --- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\) --- is the required relative error, \(s_i\) is a vector of scaling --- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and --- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\). --- --- [ARKode](https://computation.llnl.gov/projects/sundials/arkode) --- allows the user to control the step size adjustment using --- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with --- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl), --- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no --- effect. -data 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 - | 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 - | XX' Double Double Double Double -- ^ include both via relative tolerance - -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\) - | 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 @@ -module Numeric.Sundials.ODEOpts where - -import Data.Word (Word32) -import qualified Data.Vector.Storable as VS - -import Numeric.LinearAlgebra.HMatrix (Vector, Matrix) - - -type Jacobian = Double -> Vector Double -> Matrix Double - -data ODEOpts = ODEOpts { - maxNumSteps :: Word32 - , minStep :: Double - , relTol :: Double - , absTols :: VS.Vector Double - , initStep :: Maybe Double - , maxFail :: Word32 - } deriving (Read, Show, Eq, Ord) - -data SundialsDiagnostics = SundialsDiagnostics { - aRKodeGetNumSteps :: Int - , aRKodeGetNumStepAttempts :: Int - , aRKodeGetNumRhsEvals_fe :: Int - , aRKodeGetNumRhsEvals_fi :: Int - , aRKodeGetNumLinSolvSetups :: Int - , aRKodeGetNumErrTestFails :: Int - , aRKodeGetNumNonlinSolvIters :: Int - , aRKodeGetNumNonlinSolvConvFails :: Int - , aRKDlsGetNumJacEvals :: Int - , aRKDlsGetNumRhsEvals :: Int - } deriving Show - 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 @@ -#include -#include -#include /* prototypes for ARKODE fcts., consts. */ -#include /* serial N_Vector types, fcts., macros */ -#include /* access to dense SUNMatrix */ -#include /* access to dense SUNLinearSolver */ -#include /* access to ARKDls interface */ -#include /* definition of type realtype */ -#include - -/* Check function return value... - opt == 0 means SUNDIALS function allocates memory so check if - returned NULL pointer - opt == 1 means SUNDIALS function returns a flag so check if - flag >= 0 - opt == 2 means function allocates memory so check if returned - NULL pointer -*/ -int check_flag(void *flagvalue, const char *funcname, int opt) -{ - int *errflag; - - /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ - if (opt == 0 && flagvalue == NULL) { - fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", - funcname); - return 1; } - - /* Check if flag < 0 */ - else if (opt == 1) { - errflag = (int *) flagvalue; - if (*errflag < 0) { - fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", - funcname, *errflag); - return 1; }} - - /* Check if function returned NULL pointer - no memory allocated */ - else if (opt == 2 && flagvalue == NULL) { - fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", - funcname); - return 1; } - - return 0; -} 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 @@ -/* Check function return value... - opt == 0 means SUNDIALS function allocates memory so check if - returned NULL pointer - opt == 1 means SUNDIALS function returns a flag so check if - flag >= 0 - opt == 2 means function allocates memory so check if returned - NULL pointer -*/ -int check_flag(void *flagvalue, const char *funcname, int opt); diff --git a/stack.yaml b/stack.yaml index 9b858c0..619860b 100644 --- a/stack.yaml +++ b/stack.yaml @@ -13,7 +13,6 @@ packages: - packages/gsl/ - packages/glpk/ - packages/base/ -- packages/sundials/ - examples/ extra-deps: - diagrams-rasterific-1.4 -- cgit v1.2.3