1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
|
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE OverloadedStrings #-}
import qualified Language.C.Inline as C
import qualified Language.C.Inline.Unsafe as CU
import Data.Monoid ((<>))
import Foreign.C.Types
import Foreign.Ptr (Ptr)
import Foreign.Marshal.Array
import qualified Data.Vector.Storable as V
import Data.Coerce (coerce)
import Data.Monoid ((<>))
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Mutable as VM
import Foreign.C.Types
import Foreign.ForeignPtr (newForeignPtr_)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable)
import qualified Language.C.Inline as C
import qualified Language.C.Inline.Unsafe as CU
import System.IO.Unsafe (unsafePerformIO)
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
C.context (C.baseCtx <> C.vecCtx <> C.funCtx)
-- C includes
C.include "<stdio.h>"
C.include "<math.h>"
C.include "<arkode/arkode.h>" -- prototypes for ARKODE 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 "<arkode/arkode_direct.h>" -- access to ARKDls interface
C.include "<sundials/sundials_types.h>" -- definition of type realtype
C.include "<sundials/sundials_math.h>"
C.include "helpers.h"
-- 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
solve :: CDouble -> CInt
solve lambda = unsafePerformIO $ do
res <- [C.block| int { /* general problem variables */
int flag; /* reusable error-checking flag */
N_Vector y = NULL; /* empty vector for storing solution */
SUNMatrix A = NULL; /* empty matrix for linear solver */
SUNLinearSolver LS = NULL; /* empty linear solver object */
void *arkode_mem = NULL; /* empty ARKode memory structure */
FILE *UFID;
realtype t, tout;
long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
/* general problem parameters */
realtype T0 = RCONST(0.0); /* initial time */
realtype Tf = RCONST(10.0); /* final time */
realtype dTout = RCONST(1.0); /* time between outputs */
sunindextype NEQ = 1; /* number of dependent vars. */
realtype reltol = 1.0e-6; /* tolerances */
realtype abstol = 1.0e-10;
realtype lamda = -100.0; /* stiffness parameter */
/* Initial diagnostics output */
printf("\nAnalytical ODE test problem:\n");
printf(" lamda = %"GSYM"\n", lamda);
printf(" lambda = %"GSYM"\n", $(double lambda));
printf(" reltol = %.1"ESYM"\n", reltol);
printf(" abstol = %.1"ESYM"\n\n",abstol);
/* Initialize data structures */
y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
N_VConst(0.0, y); /* Specify initial condition */
arkode_mem = ARKodeCreate(); /* Create the solver memory */
if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
/* Call ARKodeInit to initialize the integrator memory and specify the */
/* right-hand side function in y'=f(t,y), the inital time T0, and */
/* the initial dependent variable vector y. Note: since this */
/* problem is fully implicit, we set f_E to NULL and f_I to f. */
flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
if (check_flag(&flag, "ARKodeInit", 1)) return 1;
/* Set routines */
flag = ARKodeSetUserData(arkode_mem, (void *) &lamda); /* Pass lamda to user functions */
if (check_flag(&flag, "ARKodeSetUserData", 1)) return 1;
flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
if (check_flag(&flag, "ARKodeSStolerances", 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;
/* Linear solver interface */
flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */
if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
flag = ARKDlsSetJacFn(arkode_mem, Jac); /* Set Jacobian routine */
if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
/* Specify linearly implicit RHS, with non-time-dependent Jacobian */
flag = ARKodeSetLinear(arkode_mem, 0);
if (check_flag(&flag, "ARKodeSetLinear", 1)) return 1;
/* Open output stream for results, output comment line */
UFID = fopen("solution.txt","w");
fprintf(UFID,"# t u\n");
/* output initial condition to disk */
fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", T0, NV_Ith_S(y,0));
/* Main time-stepping loop: calls ARKode to perform the integration, then
prints results. Stops when the final time has been reached */
t = T0;
tout = T0+dTout;
printf(" t u\n");
printf(" ---------------------\n");
while (Tf - t > 1.0e-15) {
flag = ARKode(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */
if (check_flag(&flag, "ARKode", 1)) break;
printf(" %10.6"FSYM" %10.6"FSYM"\n", t, NV_Ith_S(y,0)); /* access/print solution */
fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", t, NV_Ith_S(y,0));
if (flag >= 0) { /* successful solve: update time */
tout += dTout;
tout = (tout > Tf) ? Tf : tout;
} else { /* unsuccessful solve: break */
fprintf(stderr,"Solver failure, stopping integration\n");
break;
}
}
printf(" ---------------------\n");
fclose(UFID);
/* Get/print some final statistics on how the solve progressed */
flag = ARKodeGetNumSteps(arkode_mem, &nst);
check_flag(&flag, "ARKodeGetNumSteps", 1);
flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
printf("\nFinal Solver Statistics:\n");
printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a);
printf(" Total RHS evals: Fe = %li, Fi = %li\n", nfe, nfi);
printf(" Total linear solver setups = %li\n", nsetups);
printf(" Total RHS evals for setting up the linear system = %li\n", nfeLS);
printf(" Total number of Jacobian evaluations = %li\n", nje);
printf(" Total number of Newton iterations = %li\n", nni);
printf(" Total number of linear solver convergence failures = %li\n", ncfn);
printf(" Total number of error test failures = %li\n\n", netf);
/* check the solution error */
flag = check_ans(y, t, reltol, abstol);
/* Clean up and return */
N_VDestroy(y); /* Free y vector */
ARKodeFree(&arkode_mem); /* Free integrator memory */
SUNLinSolFree(LS); /* Free linear solver */
SUNMatDestroy(A); /* Free A matrix */
return flag;
} |]
return res
main = do
let res = solve (coerce (100.0 :: Double))
-- res <- [C.block| int { /* general problem variables */
-- int flag; /* reusable error-checking flag */
-- N_Vector y = NULL; /* empty vector for storing solution */
-- SUNMatrix A = NULL; /* empty matrix for linear solver */
-- SUNLinearSolver LS = NULL; /* empty linear solver object */
-- void *arkode_mem = NULL; /* empty ARKode memory structure */
-- FILE *UFID;
-- realtype t, tout;
-- long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
-- /* general problem parameters */
-- realtype T0 = RCONST(0.0); /* initial time */
-- realtype Tf = RCONST(10.0); /* final time */
-- realtype dTout = RCONST(1.0); /* time between outputs */
-- sunindextype NEQ = 1; /* number of dependent vars. */
-- realtype reltol = 1.0e-6; /* tolerances */
-- realtype abstol = 1.0e-10;
-- realtype lamda = -100.0; /* stiffness parameter */
-- /* Initial diagnostics output */
-- printf("\nAnalytical ODE test problem:\n");
-- printf(" lamda = %"GSYM"\n", lamda);
-- printf(" reltol = %.1"ESYM"\n", reltol);
-- printf(" abstol = %.1"ESYM"\n\n",abstol);
-- /* Initialize data structures */
-- y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
-- if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
-- N_VConst(0.0, y); /* Specify initial condition */
-- arkode_mem = ARKodeCreate(); /* Create the solver memory */
-- if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
-- /* Call ARKodeInit to initialize the integrator memory and specify the */
-- /* right-hand side function in y'=f(t,y), the inital time T0, and */
-- /* the initial dependent variable vector y. Note: since this */
-- /* problem is fully implicit, we set f_E to NULL and f_I to f. */
-- flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
-- if (check_flag(&flag, "ARKodeInit", 1)) return 1;
-- /* Set routines */
-- flag = ARKodeSetUserData(arkode_mem, (void *) &lamda); /* Pass lamda to user functions */
-- if (check_flag(&flag, "ARKodeSetUserData", 1)) return 1;
-- flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
-- if (check_flag(&flag, "ARKodeSStolerances", 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;
-- /* Linear solver interface */
-- flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */
-- if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
-- flag = ARKDlsSetJacFn(arkode_mem, Jac); /* Set Jacobian routine */
-- if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
-- /* Specify linearly implicit RHS, with non-time-dependent Jacobian */
-- flag = ARKodeSetLinear(arkode_mem, 0);
-- if (check_flag(&flag, "ARKodeSetLinear", 1)) return 1;
-- /* Open output stream for results, output comment line */
-- UFID = fopen("solution.txt","w");
-- fprintf(UFID,"# t u\n");
-- /* output initial condition to disk */
-- fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", T0, NV_Ith_S(y,0));
-- /* Main time-stepping loop: calls ARKode to perform the integration, then
-- prints results. Stops when the final time has been reached */
-- t = T0;
-- tout = T0+dTout;
-- printf(" t u\n");
-- printf(" ---------------------\n");
-- while (Tf - t > 1.0e-15) {
-- flag = ARKode(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */
-- if (check_flag(&flag, "ARKode", 1)) break;
-- printf(" %10.6"FSYM" %10.6"FSYM"\n", t, NV_Ith_S(y,0)); /* access/print solution */
-- fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", t, NV_Ith_S(y,0));
-- if (flag >= 0) { /* successful solve: update time */
-- tout += dTout;
-- tout = (tout > Tf) ? Tf : tout;
-- } else { /* unsuccessful solve: break */
-- fprintf(stderr,"Solver failure, stopping integration\n");
-- break;
-- }
-- }
-- printf(" ---------------------\n");
-- fclose(UFID);
-- /* Get/print some final statistics on how the solve progressed */
-- flag = ARKodeGetNumSteps(arkode_mem, &nst);
-- check_flag(&flag, "ARKodeGetNumSteps", 1);
-- flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
-- check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
-- flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
-- check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
-- flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
-- check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
-- flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
-- check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
-- flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
-- check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
-- flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
-- check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
-- flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
-- check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
-- flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
-- check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
-- printf("\nFinal Solver Statistics:\n");
-- printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a);
-- printf(" Total RHS evals: Fe = %li, Fi = %li\n", nfe, nfi);
-- printf(" Total linear solver setups = %li\n", nsetups);
-- printf(" Total RHS evals for setting up the linear system = %li\n", nfeLS);
-- printf(" Total number of Jacobian evaluations = %li\n", nje);
-- printf(" Total number of Newton iterations = %li\n", nni);
-- printf(" Total number of linear solver convergence failures = %li\n", ncfn);
-- printf(" Total number of error test failures = %li\n\n", netf);
-- /* check the solution error */
-- flag = check_ans(y, t, reltol, abstol);
-- /* Clean up and return */
-- N_VDestroy(y); /* Free y vector */
-- ARKodeFree(&arkode_mem); /* Free integrator memory */
-- SUNLinSolFree(LS); /* Free linear solver */
-- SUNMatDestroy(A); /* Free A matrix */
-- return flag;
-- } |]
putStrLn $ show res
|
