summaryrefslogtreecommitdiff
path: root/packages/sundials
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials')
-rw-r--r--packages/sundials/LICENSE2
-rw-r--r--packages/sundials/hmatrix-sundials.cabal14
-rw-r--r--packages/sundials/src/Main.hs34
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs (renamed from packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs)31
-rw-r--r--packages/sundials/stack.yaml66
5 files changed, 110 insertions, 37 deletions
diff --git a/packages/sundials/LICENSE b/packages/sundials/LICENSE
index e0cc65d..18e5956 100644
--- a/packages/sundials/LICENSE
+++ b/packages/sundials/LICENSE
@@ -1,4 +1,4 @@
1Copyright (c) 2018 Novadiscovery, Dominic Steinitz 1Copyright (c) 2018, Dominic Steinitz
2 2
3All rights reserved. 3All rights reserved.
4 4
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
index 1ddede3..762537e 100644
--- a/packages/sundials/hmatrix-sundials.cabal
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -1,18 +1,23 @@
1-- Initial sundials.cabal generated by cabal init. For further
2-- documentation, see http://haskell.org/cabal/users-guide/
3
1name: hmatrix-sundials 4name: hmatrix-sundials
2version: 0.1.0.0 5version: 0.1.0.0
3homepage: https://github.com/albertoruiz/hmatrix/tree/sundials 6-- synopsis:
7-- description:
8homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials
4license: BSD3 9license: BSD3
5license-file: LICENSE 10license-file: LICENSE
6author: Dominic Steinitz 11author: Dominic Steinitz
7maintainer: dominic@steinitz.org 12maintainer: dominic@steinitz.org
13-- copyright:
8category: Math 14category: Math
9build-type: Simple 15build-type: Simple
10extra-source-files: ChangeLog.md 16extra-source-files: ChangeLog.md
11cabal-version: >=1.10 17cabal-version: >=1.10
12 18
13extra-source-files: 19extra-source-files: src/helpers.c, src/helpers.h
14 src/helpers.c 20
15 src/helpers.h
16 21
17library 22library
18 build-depends: base >=4.10 && <4.11, 23 build-depends: base >=4.10 && <4.11,
@@ -36,7 +41,6 @@ executable sundials
36 template-haskell >=2.12 && <2.13, 41 template-haskell >=2.12 && <2.13,
37 containers >=0.5 && <0.6, 42 containers >=0.5 && <0.6,
38 hmatrix>=0.18, 43 hmatrix>=0.18,
39 lens,
40 plots, 44 plots,
41 diagrams-lib, 45 diagrams-lib,
42 diagrams-rasterific 46 diagrams-rasterific
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index d1f6755..5e51372 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -1,5 +1,6 @@
1module Main where 1{-# OPTIONS_GHC -Wall #-}
2 2
3import qualified Data.Vector.Storable as V
3import Numeric.Sundials.Arkode.ODE 4import Numeric.Sundials.Arkode.ODE
4import Numeric.LinearAlgebra 5import Numeric.LinearAlgebra
5 6
@@ -8,10 +9,11 @@ import qualified Diagrams.Prelude as D
8import Diagrams.Backend.Rasterific 9import Diagrams.Backend.Rasterific
9 10
10import Control.Lens 11import Control.Lens
12import Data.List (zip4)
11 13
12brusselator :: a -> [Double] -> [Double] 14
13brusselator _t x = [ a - (w + 1) * u + v * u^(2::Int) 15brusselator _t x = [ a - (w + 1) * u + v * u^2
14 , w * u - v * u^(2::Int) 16 , w * u - v * u^2
15 , (b - w) / eps - w * u 17 , (b - w) / eps - w * u
16 ] 18 ]
17 where 19 where
@@ -22,7 +24,6 @@ brusselator _t x = [ a - (w + 1) * u + v * u^(2::Int)
22 v = x !! 1 24 v = x !! 1
23 w = x !! 2 25 w = x !! 2
24 26
25stiffish :: Double -> [Double] -> [Double]
26stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] 27stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
27 where 28 where
28 lamda = -100.0 29 lamda = -100.0
@@ -44,14 +45,15 @@ kSaxis xs = P.r2Axis &~ do
44 45
45main :: IO () 46main :: IO ()
46main = do 47main = do
47 do let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) 48 let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
48 putStrLn $ show res 49 putStrLn $ show res
49 renderRasterific "diagrams/brusselator.png" 50 renderRasterific "diagrams/brusselator.png"
50 (D.dims2D 500.0 500.0) 51 (D.dims2D 500.0 500.0)
51 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res)) 52 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res))
52 53
53 do let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) 54 let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
54 putStrLn $ show res 55 putStrLn $ show res
55 renderRasterific "diagrams/stiffish.png" 56 renderRasterific "diagrams/stiffish.png"
56 (D.dims2D 500.0 500.0) 57 (D.dims2D 500.0 500.0)
57 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res)) 58 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res))
59
diff --git a/packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 6d9a1b2..f432951 100644
--- a/packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -1,14 +1,14 @@
1{-# OPTIONS_GHC -Wall #-}
2
3{-# LANGUAGE QuasiQuotes #-}
4{-# LANGUAGE TemplateHaskell #-}
1{-# LANGUAGE MultiWayIf #-} 5{-# LANGUAGE MultiWayIf #-}
2{-# LANGUAGE OverloadedStrings #-} 6{-# LANGUAGE OverloadedStrings #-}
3{-# LANGUAGE QuasiQuotes #-}
4{-# LANGUAGE ScopedTypeVariables #-} 7{-# LANGUAGE ScopedTypeVariables #-}
5{-# LANGUAGE TemplateHaskell #-}
6 8
7module Numeric.Sundials.Arkode.ODE 9module Numeric.Sundials.Arkode.ODE ( solveOde
8 ( SundialsDiagnostics(..) 10 , odeSolve
9 , solveOde 11 ) where
10 , odeSolve
11 ) where
12 12
13import qualified Language.C.Inline as C 13import qualified Language.C.Inline as C
14import qualified Language.C.Inline.Unsafe as CU 14import qualified Language.C.Inline.Unsafe as CU
@@ -40,10 +40,10 @@ C.include "<stdio.h>"
40C.include "<math.h>" 40C.include "<math.h>"
41C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. 41C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts.
42C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros 42C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros
43C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix 43C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
44C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver 44C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
45C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface 45C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface
46C.include "<sundials/sundials_types.h>" -- definition of type realtype 46C.include "<sundials/sundials_types.h>" -- definition of type realtype
47C.include "<sundials/sundials_math.h>" 47C.include "<sundials/sundials_math.h>"
48C.include "../../../helpers.h" 48C.include "../../../helpers.h"
49 49
@@ -108,7 +108,7 @@ odeSolve f y0 ts = case solveOde g (V.fromList y0) (V.fromList $ toList ts) of
108 nC = length y0 108 nC = length y0
109 g t x0 = V.fromList $ f t (V.toList x0) 109 g t x0 = V.fromList $ f t (V.toList x0)
110 110
111solveOde :: 111solveOde ::
112 (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 112 (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
113 -> V.Vector Double -- ^ Initial conditions 113 -> V.Vector Double -- ^ Initial conditions
114 -> V.Vector Double -- ^ Desired solution times 114 -> V.Vector Double -- ^ Desired solution times
@@ -116,7 +116,7 @@ solveOde ::
116solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of 116solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of
117 Left c -> Left $ fromIntegral c 117 Left c -> Left $ fromIntegral c
118 Right (v, d) -> Right (coerce v, d) 118 Right (v, d) -> Right (coerce v, d)
119 119
120solveOdeC :: 120solveOdeC ::
121 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 121 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
122 -> V.Vector CDouble -- ^ Initial conditions 122 -> V.Vector CDouble -- ^ Initial conditions
@@ -285,13 +285,13 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
285 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS); 285 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
286 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1); 286 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
287 ($vec-ptr:(long int *diagMut))[9] = ncfn; 287 ($vec-ptr:(long int *diagMut))[9] = ncfn;
288 288
289 /* Clean up and return */ 289 /* Clean up and return */
290 N_VDestroy(y); /* Free y vector */ 290 N_VDestroy(y); /* Free y vector */
291 ARKodeFree(&arkode_mem); /* Free integrator memory */ 291 ARKodeFree(&arkode_mem); /* Free integrator memory */
292 SUNLinSolFree(LS); /* Free linear solver */ 292 SUNLinSolFree(LS); /* Free linear solver */
293 SUNMatDestroy(A); /* Free A matrix */ 293 SUNMatDestroy(A); /* Free A matrix */
294 294
295 return flag; 295 return flag;
296 } |] 296 } |]
297 if res == 0 297 if res == 0
@@ -311,3 +311,4 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
311 return $ Right (m, d) 311 return $ Right (m, d)
312 else do 312 else do
313 return $ Left res 313 return $ Left res
314
diff --git a/packages/sundials/stack.yaml b/packages/sundials/stack.yaml
new file mode 100644
index 0000000..9c6b17c
--- /dev/null
+++ b/packages/sundials/stack.yaml
@@ -0,0 +1,66 @@
1# This file was automatically generated by 'stack init'
2#
3# Some commonly used options have been documented as comments in this file.
4# For advanced use and comprehensive documentation of the format, please see:
5# https://docs.haskellstack.org/en/stable/yaml_configuration/
6
7# Resolver to choose a 'specific' stackage snapshot or a compiler version.
8# A snapshot resolver dictates the compiler version and the set of packages
9# to be used for project dependencies. For example:
10#
11# resolver: lts-3.5
12# resolver: nightly-2015-09-21
13# resolver: ghc-7.10.2
14# resolver: ghcjs-0.1.0_ghc-7.10.2
15# resolver:
16# name: custom-snapshot
17# location: "./custom-snapshot.yaml"
18resolver: lts-10.9
19
20# User packages to be built.
21# Various formats can be used as shown in the example below.
22#
23# packages:
24# - some-directory
25# - https://example.com/foo/bar/baz-0.0.2.tar.gz
26# - location:
27# git: https://github.com/commercialhaskell/stack.git
28# commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a
29# - location: https://github.com/commercialhaskell/stack/commit/e7b331f14bcffb8367cd58fbfc8b40ec7642100a
30# extra-dep: true
31# subdirs:
32# - auto-update
33# - wai
34#
35# A package marked 'extra-dep: true' will only be built if demanded by a
36# non-dependency (i.e. a user package), and its test suites and benchmarks
37# will not be run. This is useful for tweaking upstream packages.
38packages:
39- .
40# Dependency packages to be pulled from upstream that are not in the resolver
41# (e.g., acme-missiles-0.3)
42# extra-deps: []
43
44# Override default flag values for local packages and extra-deps
45# flags: {}
46
47# Extra package databases containing global packages
48# extra-package-dbs: []
49
50# Control whether we use the GHC we find on the path
51# system-ghc: true
52#
53# Require a specific version of stack, using version ranges
54# require-stack-version: -any # Default
55# require-stack-version: ">=1.6"
56#
57# Override the architecture used by stack, especially useful on Windows
58# arch: i386
59# arch: x86_64
60#
61# Extra directories used by stack for building
62# extra-include-dirs: [/path/to/dir]
63# extra-lib-dirs: [/path/to/dir]
64#
65# Allow a newer minor version of GHC than the snapshot specifies
66# compiler-check: newer-minor \ No newline at end of file