hmatrix-sundials 0.19.0.0 → 0.19.1.0
raw patch · 14 files changed
+1273/−333 lines, 14 filesbinary-added
Files
- ChangeLog.md +10/−2
- README.md +4/−2
- diagrams/lorenzA.png binary
- diagrams/predatorPrey.png binary
- diagrams/predatorPrey1.png binary
- diagrams/predatorPrey2.png binary
- hmatrix-sundials.cabal +12/−8
- src/Arkode.hsc +0/−114
- src/Main.hs +60/−12
- src/Numeric/Sundials/ARKode/ODE.hs +153/−155
- src/Numeric/Sundials/Arkode.hsc +209/−0
- src/Numeric/Sundials/CVode/ODE.hs +793/−0
- src/Numeric/Sundials/ODEOpts.hs +32/−0
- src/Types.hs +0/−40
ChangeLog.md view
@@ -1,5 +1,13 @@ # Revision history for hmatrix-sundials -## 0.1.0.0 -- 2018-04-21+## 0.19.0.0 -- 2018-04-21 -* First version. Released on an unsuspecting world. Just Runge-Kutta methods to start with.+* First version. Released on an unsuspecting world. Just Runge-Kutta+ methods to start with.++## 0.19.1.0 -- 2018-08-14++* Methods in CVODE: Adams-Moulton and Backward Differentiation+ Formulas (BDFs).+* An experimental interface for finding when a variable has crossed a+ certain value.
README.md view
@@ -1,5 +1,7 @@-Currently only an interface to the Runge-Kutta methods:-[ARKode](https://computation.llnl.gov/projects/sundials/arkode)+An interface to the Runge-Kutta methods:+[ARKode](https://computation.llnl.gov/projects/sundials/arkode) and+the methods in+[CVode](https://computation.llnl.gov/projects/sundials/cvode) The interface is almost certainly going to change. Sundials gives a rich set of "combinators" for controlling the solution of your problem
− diagrams/lorenzA.png
binary file changed (19746 → absent bytes)
+ diagrams/predatorPrey.png view
binary file changed (absent → 18107 bytes)
+ diagrams/predatorPrey1.png view
binary file changed (absent → 13377 bytes)
+ diagrams/predatorPrey2.png view
binary file changed (absent → 17022 bytes)
hmatrix-sundials.cabal view
@@ -1,5 +1,5 @@ name: hmatrix-sundials-version: 0.19.0.0+version: 0.19.1.0 synopsis: hmatrix interface to sundials description: An interface to the solving suite SUNDIALS. Currently, it mimics the solving interface in hmstrix-gsl but@@ -25,21 +25,24 @@ template-haskell >=2.12 && <2.13, containers >=0.5 && <0.6, hmatrix>=0.18- extra-libraries: sundials_arkode+ extra-libraries: sundials_arkode,+ sundials_cvode other-extensions: QuasiQuotes hs-source-dirs: src- exposed-modules: Numeric.Sundials.ARKode.ODE- other-modules: Types,- Arkode+ 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: Types,+ other-modules: Numeric.Sundials.ODEOpts, Numeric.Sundials.ARKode.ODE,- Arkode+ 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,@@ -52,6 +55,7 @@ lens, hspec hs-source-dirs: src- extra-libraries: sundials_arkode+ extra-libraries: sundials_arkode,+ sundials_cvode c-sources: src/helpers.c src/helpers.h default-language: Haskell2010
− src/Arkode.hsc
@@ -1,114 +0,0 @@-module Arkode where--import Foreign-import Foreign.C.Types---#include <stdio.h>-#include <sundials/sundials_nvector.h>-#include <sundials/sundials_matrix.h>-#include <nvector/nvector_serial.h>-#include <sunmatrix/sunmatrix_dense.h>-#include <arkode/arkode.h>---#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--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
src/Main.hs view
@@ -1,6 +1,7 @@ {-# OPTIONS_GHC -Wall #-} -import Numeric.Sundials.ARKode.ODE+import qualified Numeric.Sundials.ARKode.ODE as ARK+import qualified Numeric.Sundials.CVode.ODE as CV import Numeric.LinearAlgebra import Plots as P@@ -80,6 +81,23 @@ 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@@ -97,34 +115,38 @@ main :: IO () main = do - let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])+ 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 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])+ 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 = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])+ 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 = odeSolveV (SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])+ 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 = odeSolveV (TRBDF2_3_3_2') 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 maxDiff = maximum $ map abs $- zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)+ let maxDiffA = maximum $ map abs $+ zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) - hspec $ describe "Compare results" $ do- it "for two different RK methods" $- maxDiff < 1.0e-6+ let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) - let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 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))@@ -136,3 +158,29 @@ 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+
src/Numeric/Sundials/ARKode/ODE.hs view
@@ -1,5 +1,3 @@-{-# OPTIONS_GHC -Wall #-}- {-# LANGUAGE QuasiQuotes #-} {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE MultiWayIf #-}@@ -22,8 +20,7 @@ -- Stability : provisional -- -- Solution of ordinary differential equation (ODE) initial value problems.------ <https://computation.llnl.gov/projects/sundials/sundials-software>+-- See <https://computation.llnl.gov/projects/sundials/sundials-software> for more detail. -- -- A simple example: --@@ -67,6 +64,54 @@ -- (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 -- -- \[@@ -116,8 +161,6 @@ , butcherTable , ODEMethod(..) , StepControl(..)- , Jacobian- , SundialsDiagnostics(..) ) where import qualified Language.C.Inline as C@@ -126,27 +169,50 @@ import Data.Monoid ((<>)) import Data.Maybe (isJust) -import Foreign.C.Types+import Foreign.C.Types (CDouble, CInt, CLong) import Foreign.Ptr (Ptr)-import Foreign.ForeignPtr (newForeignPtr_)-import Foreign.Storable (Storable)+import Foreign.Storable (poke) import qualified Data.Vector.Storable as V-import qualified Data.Vector.Storable.Mutable as VM import Data.Coerce (coerce) import System.IO.Unsafe (unsafePerformIO)-import GHC.Generics+import GHC.Generics (C1, Constructor, (:+:)(..), D1, Rep, Generic, M1(..),+ from, conName) import Numeric.LinearAlgebra.Devel (createVector) -import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),- subMatrix, rows, cols, toLists,- size, subVector)+import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,+ cols, toLists, size, reshape,+ subVector, subMatrix, (><)) -import qualified Types as T-import Arkode-import qualified Arkode as B+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)@@ -162,70 +228,9 @@ C.include "<sundials/sundials_types.h>" -- definition of type realtype C.include "<sundials/sundials_math.h>" C.include "../../../helpers.h"-C.include "Arkode_hsc.h"+C.include "Numeric/Sundials/Arkode_hsc.h" -getDataFromContents :: Int -> Ptr T.SunVector -> IO (V.Vector CDouble)-getDataFromContents len ptr = do- qtr <- B.getContentPtr ptr- rtr <- B.getData qtr- vectorFromC len rtr---- FIXME: Potentially an instance of Storable-_getMatrixDataFromContents :: Ptr T.SunMatrix -> IO T.SunMatrix-_getMatrixDataFromContents ptr = do- qtr <- B.getContentMatrixPtr ptr- rs <- B.getNRows qtr- cs <- B.getNCols qtr- rtr <- B.getMatrixData qtr- vs <- vectorFromC (fromIntegral $ rs * cs) rtr- return $ T.SunMatrix { T.rows = rs, T.cols = cs, T.vals = vs }--putMatrixDataFromContents :: T.SunMatrix -> Ptr T.SunMatrix -> IO ()-putMatrixDataFromContents mat ptr = do- let rs = T.rows mat- cs = T.cols mat- vs = T.vals mat- qtr <- B.getContentMatrixPtr ptr- B.putNRows rs qtr- B.putNCols cs qtr- rtr <- B.getMatrixData qtr- vectorToC vs (fromIntegral $ rs * cs) rtr--- FIXME: END--putDataInContents :: Storable a => V.Vector a -> Int -> Ptr b -> IO ()-putDataInContents vec len ptr = do- qtr <- B.getContentPtr ptr- rtr <- B.getData qtr- vectorToC vec len rtr---- Utils--vectorFromC :: Storable a => Int -> Ptr a -> IO (V.Vector a)-vectorFromC len ptr = do- ptr' <- newForeignPtr_ ptr- V.freeze $ VM.unsafeFromForeignPtr0 ptr' len--vectorToC :: Storable a => V.Vector a -> Int -> Ptr a -> IO ()-vectorToC vec len ptr = do- ptr' <- newForeignPtr_ ptr- V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec--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--type Jacobian = Double -> Vector Double -> Matrix Double- -- | Stepping functions data ODEMethod = SDIRK_2_1_2 Jacobian | SDIRK_2_1_2'@@ -390,15 +395,9 @@ -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution odeSolveV meth hi epsAbs epsRel f y0 ts =- case odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts of- Left c -> error $ show c -- FIXME- -- FIXME: Can we do better than using lists?- Right (v, _d) -> (nR >< nC) (V.toList v)- where- us = toList ts- nR = length us- nC = size y0- g t x0 = coerce $ f t x0+ 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@@ -410,16 +409,11 @@ -> Matrix Double -- ^ solution odeSolve f y0 ts = -- FIXME: These tolerances are different from the ones in GSL- case odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts) of- Left c -> error $ show c -- FIXME- Right (v, _d) -> (nR >< nC) (V.toList v)+ odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1) Nothing g (V.fromList y0) (V.fromList $ toList ts) where- us = toList ts- nR = length us- nC = length y0 g t x0 = V.fromList $ f t (V.toList x0) -odeSolveVWith' ::+odeSolveVWith :: ODEMethod -> StepControl -> Maybe Double -- ^ initial step size - by default, ARKode@@ -432,16 +426,22 @@ -> 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 method control initStepSize f y0 tt of- Left c -> error $ show c -- FIXME- Right (v, _d) -> (nR >< nC) (V.toList v)+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- nR = V.length tt- nC = V.length y0+ opts = ODEOpts { maxNumSteps = 10000+ , minStep = 1.0e-12+ , relTol = error "relTol"+ , absTols = error "absTol"+ , initStep = error "initStep"+ , maxFail = 10+ } -odeSolveVWith ::- ODEMethod+odeSolveVWith' ::+ ODEOpts+ -> ODEMethod -> StepControl -> Maybe Double -- ^ initial step size - by default, ARKode -- estimates the initial step size to be the@@ -452,19 +452,21 @@ -> (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 Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution-odeSolveVWith method control initStepSize f y0 tt =- case solveOdeC (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)+ -> 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 c -> Left $ fromIntegral c- Right (v, d) -> Right (coerce v, d)+ 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 absTol relTol) = coerce (V.replicate l absTol, relTol)- scise (X' absTol relTol) = coerce (V.replicate l absTol, relTol)- scise (XX' absTol relTol yScale _yDotScale) = coerce (V.replicate l absTol, yScale * relTol)+ 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' absTol relTol yScale _yDotScale ss) = coerce (V.map (* absTol) ss, yScale * relTol)+ 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 }@@ -476,14 +478,19 @@ 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 CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution-solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO $ do+ -> 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@@ -494,22 +501,18 @@ -- used :( Nothing -> 0.0 Just x -> x+ let dim = V.length f0 nEq :: CLong nEq = fromIntegral dim nTs :: CInt nTs = fromIntegral $ V.length ts- -- FIXME: fMut is not actually mutatated- fMut <- V.thaw f0- tMut <- V.thaw ts- -- FIXME: I believe this gets taken from the ghc heap and so should- -- be subject to garbage collection. 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 'Types'. FIXME: The Haskell type is currently empty!+ -- 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@@ -529,7 +532,7 @@ case jacH of Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined" Just jacI -> do j <- jacI t <$> getDataFromContents dim y- putMatrixDataFromContents j jacS+ 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 } |]@@ -549,7 +552,7 @@ /* general problem parameters */ - realtype T0 = RCONST(($vec-ptr:(double *tMut))[0]); /* initial time */+ realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */ sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */ /* Initialize data structures */@@ -558,14 +561,14 @@ 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 *fMut))[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 *absTols))[i];+ NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i]; }; arkode_mem = ARKodeCreate(); /* Create the solver memory */@@ -577,7 +580,7 @@ /* 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 Types */+ /* 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;@@ -586,14 +589,15 @@ if (check_flag(&flag, "ARKodeInit", 1)) return 1; } - /* FIXME: A hack for initial testing */- flag = ARKodeSetMinStep(arkode_mem, 1.0e-12);+ flag = ARKodeSetMinStep(arkode_mem, $(double minStep_)); if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1;- flag = ARKodeSetMaxNumSteps(arkode_mem, 10000);+ 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 relTol), tv);+ flag = ARKodeSVtolerances(arkode_mem, $(double rTol), tv); if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1; /* Initialize dense matrix data structure and solver */@@ -638,19 +642,13 @@ /* Stops when the final time has been reached */ for (i = 1; i < $(int nTs); i++) { - flag = ARKode(arkode_mem, ($vec-ptr:(double *tMut))[i], y, &t, ARK_NORMAL); /* call integrator */- if (check_flag(&flag, "ARKode", 1)) break;+ 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); }-- /* unsuccessful solve: break */- if (flag < 0) {- fprintf(stderr,"Solver failure, stopping integration\n");- break;- } } /* Get some final statistics on how the solve progressed */@@ -701,23 +699,23 @@ 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- 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 return $ Right (m, d) else do- return $ Left res+ return $ Left (m, res) data ButcherTable = ButcherTable { am :: Matrix Double , cv :: Vector Double@@ -738,7 +736,7 @@ case getBT method of Left c -> error $ show c -- FIXME Right (ButcherTable' v w x y, sqp) ->- ButcherTable { am = subMatrix (0, 0) (s, s) $ (B.arkSMax >< B.arkSMax) (V.toList v)+ 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@@ -773,11 +771,11 @@ btSQP :: V.Vector CInt <- createVector 3 btSQPMut <- V.thaw btSQP- btAs :: V.Vector CDouble <- createVector (B.arkSMax * B.arkSMax)+ btAs :: V.Vector CDouble <- createVector (arkSMax * arkSMax) btAsMut <- V.thaw btAs- btCs :: V.Vector CDouble <- createVector B.arkSMax- btBs :: V.Vector CDouble <- createVector B.arkSMax- btB2s :: V.Vector CDouble <- createVector B.arkSMax+ 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
+ src/Numeric/Sundials/Arkode.hsc view
@@ -0,0 +1,209 @@+{-# 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 <stdio.h>+#include <sundials/sundials_nvector.h>+#include <sundials/sundials_matrix.h>+#include <nvector/nvector_serial.h>+#include <sunmatrix/sunmatrix_dense.h>+#include <arkode/arkode.h>+#include <cvode/cvode.h>+++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_SUCCESS :: Int+cV_SUCCESS = #const CV_SUCCESS+cV_ROOT_RETURN :: Int+cV_ROOT_RETURN = #const CV_ROOT_RETURN++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
+ src/Numeric/Sundials/CVode/ODE.hs view
@@ -0,0 +1,793 @@+{-# 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.+--+-- <https://computation.llnl.gov/projects/sundials/sundials-software>+--+-- A simple example:+--+-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>+--+-- @+-- 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'+ , odeSolveRootVWith'+ , ODEMethod(..)+ , StepControl(..)+ , SolverResult(..)+ ) 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,+ vectorToC, cV_SUCCESS, cV_ROOT_RETURN)+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 "<stdlib.h>"+C.include "<stdio.h>"+C.include "<math.h>"+C.include "<cvode/cvode.h>" -- prototypes for CVODE fcts., consts.+C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros+C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix+C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver+C.include "<cvode/cvode_direct.h>" -- access to CVDls interface+C.include "<sundials/sundials_types.h>" -- definition of type realtype+C.include "<sundials/sundials_math.h>"+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)++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+ -> CInt -- ^ FIXME+ -> (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ FIXME+ -> V.Vector CDouble -- ^ Desired solution times+ -> SolverResult V.Vector V.Vector CInt CDouble+solveOdeC' maxErrTestFails maxNumSteps_ minStep_ method initStepSize+ jacH (aTols, rTol) fun f0 nr g 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.+ -- FIXME: The Haskell type is currently empty!+ let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt+ funIO t y f _ptr = do+ -- Convert the pointer we get from C (y) to a vector, and then+ -- apply the user-supplied function.+ fImm <- fun t <$> 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 nrPre = fromIntegral nr+ gResults :: V.Vector CInt <- createVector nrPre+ gResMut <- V.thaw gResults+ tRoot :: V.Vector CDouble <- createVector 1+ tRootMut <- V.thaw tRoot++ let gIO :: CDouble -> Ptr T.SunVector -> Ptr CDouble -> Ptr () -> IO CInt+ gIO x y f _ptr = do+ -- Convert the pointer we get from C (y) to a vector, and then+ -- apply the user-supplied function.+ gImm <- g x <$> getDataFromContents dim y+ -- Fill in the provided pointer with the resulting vector.+ vectorToC gImm nrPre 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 flagr; /* root finding 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;++ realtype tout;++ /* 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);++ /* Call CVodeRootInit to specify the root function g with nr components */+ flag = CVodeRootInit(cvode_mem, $(int nr), $fun:(int (* gIO) (double t, SunVector y[], double gout[], void * params)));++ if (check_flag(&flag, "CVodeRootInit", 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);+ }++ if (flag == CV_ROOT_RETURN) {+ flagr = CVodeGetRootInfo(cvode_mem, ($vec-ptr:(int *gResMut)));+ if (check_flag(&flagr, "CVodeGetRootInfo", 1)) return(1);+ ($vec-ptr:(double *tRootMut))[0] = t;+ flagr = flag;+ break;+ }+ }++ /* 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 */++ if (flag == CV_SUCCESS && flagr == CV_ROOT_RETURN) {+ return CV_ROOT_RETURN;+ }+ else {+ 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+ t <- V.freeze tRootMut+ rs <- V.freeze gResMut+ putStrLn $ show rs+ let f r | r == cV_SUCCESS = SolverSuccess m d+ | r == cV_ROOT_RETURN = SolverRoot (t V.!0) rs m d+ | otherwise = SolverError m res+ return $ f $ fromIntegral res++data SolverResult f g a b =+ SolverError (f b) a -- ^ Partial results and error code+ | SolverSuccess (f b) SundialsDiagnostics -- ^ Results and diagnostics+ | SolverRoot b (g a) (f b) SundialsDiagnostics -- ^ Time at which the root was found, the root itself and the+ -- results and diagnostics. NB the final result will be at the time+ -- at which the root was found not as specified by the times given+ -- to the solver.+ deriving Show++odeSolveRootVWith' ::+ 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+ -> Int -- ^ Dimension of the range of the roots function+ -> (Double -> V.Vector Double -> V.Vector Double) -- ^ Roots function+ -> V.Vector Double -- ^ Desired solution times+ -> SolverResult Matrix Vector Int Double+odeSolveRootVWith' opts method control initStepSize f y0 is gg tt =+ case solveOdeC' (fromIntegral $ maxFail opts)+ (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)+ (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)+ (coerce f) (coerce y0) (fromIntegral is) (coerce gg) (coerce tt) of+ SolverError v c -> SolverError (reshape l (coerce v)) (fromIntegral c)+ SolverSuccess v d -> SolverSuccess (reshape l (coerce v)) d+ SolverRoot t rs v d -> SolverRoot (coerce t) (V.map fromIntegral rs) (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++-- | 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}\)
+ src/Numeric/Sundials/ODEOpts.hs view
@@ -0,0 +1,32 @@+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+
− src/Types.hs
@@ -1,40 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--{-# LANGUAGE QuasiQuotes #-}-{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE MultiWayIf #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE EmptyDataDecls #-}--module Types where--import Foreign.C.Types--import qualified Language.Haskell.TH as TH-import qualified Language.C.Types as CT-import qualified Data.Map as Map-import Language.C.Inline.Context--import qualified Data.Vector.Storable as V---data SunVector-data SunMatrix = SunMatrix { rows :: CInt- , cols :: CInt- , vals :: V.Vector CDouble- }---- FIXME: Is this true?-type SunIndexType = CLong--sunTypesTable :: Map.Map CT.TypeSpecifier TH.TypeQ-sunTypesTable = Map.fromList- [- (CT.TypeName "sunindextype", [t| SunIndexType |] )- , (CT.TypeName "SunVector", [t| SunVector |] )- , (CT.TypeName "SunMatrix", [t| SunMatrix |] )- ]--sunCtx :: Context-sunCtx = mempty {ctxTypesTable = sunTypesTable}-