summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-03-12 09:50:59 +0000
committerDominic Steinitz <dominic@steinitz.org>2018-03-12 09:50:59 +0000
commitf2b1eae3d71c546abc71e099b4bd86010627f0fb (patch)
tree286ae3211ef401db83d02e03bf4044cf01d015c2 /packages
parent7c2337e093ecd7d367d30d567bf5172ee639666b (diff)
Now builds with stack and cabal
Diffstat (limited to 'packages')
-rw-r--r--packages/sundials/src/Main.hs (renamed from packages/sundials/src/Test.hs)38
-rw-r--r--packages/sundials/src/helpers.c36
-rw-r--r--packages/sundials/src/helpers.h26
-rw-r--r--packages/sundials/sundials.cabal3
4 files changed, 98 insertions, 5 deletions
diff --git a/packages/sundials/src/Test.hs b/packages/sundials/src/Main.hs
index a99582a..9e8bc63 100644
--- a/packages/sundials/src/Test.hs
+++ b/packages/sundials/src/Main.hs
@@ -150,15 +150,45 @@ C.include "helpers.h"
150-- } 150-- }
151 151
152main = do 152main = do
153 res <- [C.block| int { sunindextype NEQ = 1; /* number of dependent vars. */ 153 res <- [C.block| int { /* general problem variables */
154 N_Vector y = NULL; /* empty vector for storing solution */ 154 int flag; /* reusable error-checking flag */
155 void *arkode_mem = NULL; /* empty ARKode memory structure */ 155 N_Vector y = NULL; /* empty vector for storing solution */
156 SUNMatrix A = NULL; /* empty matrix for linear solver */
157 SUNLinearSolver LS = NULL; /* empty linear solver object */
158 void *arkode_mem = NULL; /* empty ARKode memory structure */
159 FILE *UFID;
160 realtype t, tout;
161 long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
162
163 /* general problem parameters */
164 realtype T0 = RCONST(0.0); /* initial time */
165 realtype Tf = RCONST(10.0); /* final time */
166 realtype dTout = RCONST(1.0); /* time between outputs */
167 sunindextype NEQ = 1; /* number of dependent vars. */
168 realtype reltol = 1.0e-6; /* tolerances */
169 realtype abstol = 1.0e-10;
170 realtype lamda = -100.0; /* stiffness parameter */
171
172 /* Initial diagnostics output */
173 printf("\nAnalytical ODE test problem:\n");
174 printf(" lamda = %"GSYM"\n", lamda);
175 printf(" reltol = %.1"ESYM"\n", reltol);
176 printf(" abstol = %.1"ESYM"\n\n",abstol);
177
178 /* Initialize data structures */
156 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */ 179 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
157 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1; 180 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
158
159 N_VConst(0.0, y); /* Specify initial condition */ 181 N_VConst(0.0, y); /* Specify initial condition */
160 arkode_mem = ARKodeCreate(); /* Create the solver memory */ 182 arkode_mem = ARKodeCreate(); /* Create the solver memory */
161 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1; 183 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
184
185 /* Call ARKodeInit to initialize the integrator memory and specify the */
186 /* right-hand side function in y'=f(t,y), the inital time T0, and */
187 /* the initial dependent variable vector y. Note: since this */
188 /* problem is fully implicit, we set f_E to NULL and f_I to f. */
189 flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
190 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
191
162 return 0; 192 return 0;
163 } |] 193 } |]
164 putStrLn $ show res 194 putStrLn $ show res
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
index 3498119..3e1a9a8 100644
--- a/packages/sundials/src/helpers.c
+++ b/packages/sundials/src/helpers.c
@@ -1,5 +1,12 @@
1#include <stdio.h> 1#include <stdio.h>
2#include "helpers.h" 2#include <math.h>
3#include <arkode/arkode.h> /* prototypes for ARKODE fcts., consts. */
4#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
5#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
6#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
7#include <arkode/arkode_direct.h> /* access to ARKDls interface */
8#include <sundials/sundials_types.h> /* definition of type realtype */
9#include <sundials/sundials_math.h>
3 10
4/* Check function return value... 11/* Check function return value...
5 opt == 0 means SUNDIALS function allocates memory so check if 12 opt == 0 means SUNDIALS function allocates memory so check if
@@ -35,3 +42,30 @@ int check_flag(void *flagvalue, const char *funcname, int opt)
35 42
36 return 0; 43 return 0;
37} 44}
45
46/* f routine to compute the ODE RHS function f(t,y). */
47int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
48{
49 realtype *rdata = (realtype *) user_data; /* cast user_data to realtype */
50 realtype lamda = rdata[0]; /* set shortcut for stiffness parameter */
51 realtype u = NV_Ith_S(y,0); /* access current solution value */
52
53 /* fill in the RHS function: "NV_Ith_S" accesses the 0th entry of ydot */
54 NV_Ith_S(ydot,0) = lamda*u + 1.0/(1.0+t*t) - lamda*atan(t);
55
56 return 0; /* return with success */
57}
58
59/* Jacobian routine to compute J(t,y) = df/dy. */
60int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
61 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
62{
63 realtype *rdata = (realtype *) user_data; /* cast user_data to realtype */
64 realtype lamda = rdata[0]; /* set shortcut for stiffness parameter */
65 realtype *Jdata = SUNDenseMatrix_Data(J);
66
67 /* Fill in Jacobian of f: set the first entry of the data array to set the (0,0) entry */
68 Jdata[0] = lamda;
69
70 return 0; /* return with success */
71}
diff --git a/packages/sundials/src/helpers.h b/packages/sundials/src/helpers.h
index d574ff5..ec11355 100644
--- a/packages/sundials/src/helpers.h
+++ b/packages/sundials/src/helpers.h
@@ -1 +1,27 @@
1#if defined(SUNDIALS_EXTENDED_PRECISION)
2#define GSYM "Lg"
3#define ESYM "Le"
4#define FSYM "Lf"
5#else
6#define GSYM "g"
7#define ESYM "e"
8#define FSYM "f"
9#endif
10
11/* Check function return value...
12 opt == 0 means SUNDIALS function allocates memory so check if
13 returned NULL pointer
14 opt == 1 means SUNDIALS function returns a flag so check if
15 flag >= 0
16 opt == 2 means function allocates memory so check if returned
17 NULL pointer
18*/
1int check_flag(void *flagvalue, const char *funcname, int opt); 19int check_flag(void *flagvalue, const char *funcname, int opt);
20
21/* f routine to compute the ODE RHS function f(t,y). */
22int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
23
24/* Jacobian routine to compute J(t,y) = df/dy. */
25int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
26 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
27
diff --git a/packages/sundials/sundials.cabal b/packages/sundials/sundials.cabal
index 85a96fa..a68f33a 100644
--- a/packages/sundials/sundials.cabal
+++ b/packages/sundials/sundials.cabal
@@ -16,6 +16,8 @@ build-type: Simple
16extra-source-files: ChangeLog.md 16extra-source-files: ChangeLog.md
17cabal-version: >=1.10 17cabal-version: >=1.10
18 18
19extra-source-files: src/helpers.c, src/helpers.h
20
19executable sundials 21executable sundials
20 main-is: Main.hs 22 main-is: Main.hs
21 -- other-modules: 23 -- other-modules:
@@ -28,3 +30,4 @@ executable sundials
28 hs-source-dirs: src 30 hs-source-dirs: src
29 default-language: Haskell2010 31 default-language: Haskell2010
30 extra-libraries: sundials_arkode 32 extra-libraries: sundials_arkode
33 C-sources: src/helpers.c src/helpers.h