diff options
Diffstat (limited to 'packages/sundials')
-rw-r--r-- | packages/sundials/LICENSE | 2 | ||||
-rw-r--r-- | packages/sundials/hmatrix-sundials.cabal | 14 | ||||
-rw-r--r-- | packages/sundials/src/Main.hs | 34 | ||||
-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.yaml | 66 |
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 @@ | |||
1 | Copyright (c) 2018 Novadiscovery, Dominic Steinitz | 1 | Copyright (c) 2018, Dominic Steinitz |
2 | 2 | ||
3 | All rights reserved. | 3 | All 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 | |||
1 | name: hmatrix-sundials | 4 | name: hmatrix-sundials |
2 | version: 0.1.0.0 | 5 | version: 0.1.0.0 |
3 | homepage: https://github.com/albertoruiz/hmatrix/tree/sundials | 6 | -- synopsis: |
7 | -- description: | ||
8 | homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials | ||
4 | license: BSD3 | 9 | license: BSD3 |
5 | license-file: LICENSE | 10 | license-file: LICENSE |
6 | author: Dominic Steinitz | 11 | author: Dominic Steinitz |
7 | maintainer: dominic@steinitz.org | 12 | maintainer: dominic@steinitz.org |
13 | -- copyright: | ||
8 | category: Math | 14 | category: Math |
9 | build-type: Simple | 15 | build-type: Simple |
10 | extra-source-files: ChangeLog.md | 16 | extra-source-files: ChangeLog.md |
11 | cabal-version: >=1.10 | 17 | cabal-version: >=1.10 |
12 | 18 | ||
13 | extra-source-files: | 19 | extra-source-files: src/helpers.c, src/helpers.h |
14 | src/helpers.c | 20 | |
15 | src/helpers.h | ||
16 | 21 | ||
17 | library | 22 | library |
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 @@ | |||
1 | module Main where | 1 | {-# OPTIONS_GHC -Wall #-} |
2 | 2 | ||
3 | import qualified Data.Vector.Storable as V | ||
3 | import Numeric.Sundials.Arkode.ODE | 4 | import Numeric.Sundials.Arkode.ODE |
4 | import Numeric.LinearAlgebra | 5 | import Numeric.LinearAlgebra |
5 | 6 | ||
@@ -8,10 +9,11 @@ import qualified Diagrams.Prelude as D | |||
8 | import Diagrams.Backend.Rasterific | 9 | import Diagrams.Backend.Rasterific |
9 | 10 | ||
10 | import Control.Lens | 11 | import Control.Lens |
12 | import Data.List (zip4) | ||
11 | 13 | ||
12 | brusselator :: a -> [Double] -> [Double] | 14 | |
13 | brusselator _t x = [ a - (w + 1) * u + v * u^(2::Int) | 15 | brusselator _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 | ||
25 | stiffish :: Double -> [Double] -> [Double] | ||
26 | stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] | 27 | stiffish 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 | ||
45 | main :: IO () | 46 | main :: IO () |
46 | main = do | 47 | main = 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 | ||
7 | module Numeric.Sundials.Arkode.ODE | 9 | module Numeric.Sundials.Arkode.ODE ( solveOde |
8 | ( SundialsDiagnostics(..) | 10 | , odeSolve |
9 | , solveOde | 11 | ) where |
10 | , odeSolve | ||
11 | ) where | ||
12 | 12 | ||
13 | import qualified Language.C.Inline as C | 13 | import qualified Language.C.Inline as C |
14 | import qualified Language.C.Inline.Unsafe as CU | 14 | import qualified Language.C.Inline.Unsafe as CU |
@@ -40,10 +40,10 @@ C.include "<stdio.h>" | |||
40 | C.include "<math.h>" | 40 | C.include "<math.h>" |
41 | C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. | 41 | C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts. |
42 | C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros | 42 | C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros |
43 | C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix | 43 | C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix |
44 | C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver | 44 | C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver |
45 | C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface | 45 | C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface |
46 | C.include "<sundials/sundials_types.h>" -- definition of type realtype | 46 | C.include "<sundials/sundials_types.h>" -- definition of type realtype |
47 | C.include "<sundials/sundials_math.h>" | 47 | C.include "<sundials/sundials_math.h>" |
48 | C.include "../../../helpers.h" | 48 | C.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 | ||
111 | solveOde :: | 111 | solveOde :: |
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 :: | |||
116 | solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of | 116 | solveOde 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 | ||
120 | solveOdeC :: | 120 | solveOdeC :: |
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" | ||
18 | resolver: 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. | ||
38 | packages: | ||
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 | ||