summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-03-15 08:33:26 +0000
committerDominic Steinitz <dominic@steinitz.org>2018-03-15 08:33:26 +0000
commit0c864a1d758bdf39a2492e68130d02c05da06dc5 (patch)
tree8ffb9a3d3e1e163c6e58c817572c4123c3de469b /packages
parentd23f3abc8038e9669ef1aa6b7ab9fe5346f95410 (diff)
Do not handle Jacobians yet and tidy
Diffstat (limited to 'packages')
-rw-r--r--packages/sundials/src/Main.hs142
1 files changed, 1 insertions, 141 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index b6855cb..5972be7 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -111,14 +111,6 @@ solve lambda = unsafePerformIO $ do
111 111
112 /* Linear solver interface */ 112 /* Linear solver interface */
113 flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */ 113 flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */
114 if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
115 flag = ARKDlsSetJacFn(arkode_mem, Jac); /* Set Jacobian routine */
116 if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
117
118 /* Specify linearly implicit RHS, with non-time-dependent Jacobian */
119 flag = ARKodeSetLinear(arkode_mem, 0);
120 if (check_flag(&flag, "ARKodeSetLinear", 1)) return 1;
121
122 /* Open output stream for results, output comment line */ 114 /* Open output stream for results, output comment line */
123 UFID = fopen("solution.txt","w"); 115 UFID = fopen("solution.txt","w");
124 fprintf(UFID,"# t u\n"); 116 fprintf(UFID,"# t u\n");
@@ -193,137 +185,5 @@ solve lambda = unsafePerformIO $ do
193 return res 185 return res
194 186
195main = do 187main = do
196 let res = solve (coerce (100.0 :: Double)) 188 let res = solve undefined undefined (coerce (100.0 :: Double))
197 -- res <- [C.block| int { /* general problem variables */
198 -- int flag; /* reusable error-checking flag */
199 -- N_Vector y = NULL; /* empty vector for storing solution */
200 -- SUNMatrix A = NULL; /* empty matrix for linear solver */
201 -- SUNLinearSolver LS = NULL; /* empty linear solver object */
202 -- void *arkode_mem = NULL; /* empty ARKode memory structure */
203 -- FILE *UFID;
204 -- realtype t, tout;
205 -- long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
206
207 -- /* general problem parameters */
208 -- realtype T0 = RCONST(0.0); /* initial time */
209 -- realtype Tf = RCONST(10.0); /* final time */
210 -- realtype dTout = RCONST(1.0); /* time between outputs */
211 -- sunindextype NEQ = 1; /* number of dependent vars. */
212 -- realtype reltol = 1.0e-6; /* tolerances */
213 -- realtype abstol = 1.0e-10;
214 -- realtype lamda = -100.0; /* stiffness parameter */
215
216 -- /* Initial diagnostics output */
217 -- printf("\nAnalytical ODE test problem:\n");
218 -- printf(" lamda = %"GSYM"\n", lamda);
219 -- printf(" reltol = %.1"ESYM"\n", reltol);
220 -- printf(" abstol = %.1"ESYM"\n\n",abstol);
221
222 -- /* Initialize data structures */
223 -- y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
224 -- if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
225 -- N_VConst(0.0, y); /* Specify initial condition */
226 -- arkode_mem = ARKodeCreate(); /* Create the solver memory */
227 -- if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
228
229 -- /* Call ARKodeInit to initialize the integrator memory and specify the */
230 -- /* right-hand side function in y'=f(t,y), the inital time T0, and */
231 -- /* the initial dependent variable vector y. Note: since this */
232 -- /* problem is fully implicit, we set f_E to NULL and f_I to f. */
233 -- flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
234 -- if (check_flag(&flag, "ARKodeInit", 1)) return 1;
235
236 -- /* Set routines */
237 -- flag = ARKodeSetUserData(arkode_mem, (void *) &lamda); /* Pass lamda to user functions */
238 -- if (check_flag(&flag, "ARKodeSetUserData", 1)) return 1;
239 -- flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
240 -- if (check_flag(&flag, "ARKodeSStolerances", 1)) return 1;
241
242 -- /* Initialize dense matrix data structure and solver */
243 -- A = SUNDenseMatrix(NEQ, NEQ);
244 -- if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
245 -- LS = SUNDenseLinearSolver(y, A);
246 -- if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
247
248 -- /* Linear solver interface */
249 -- flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */
250 -- if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
251 -- flag = ARKDlsSetJacFn(arkode_mem, Jac); /* Set Jacobian routine */
252 -- if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
253
254 -- /* Specify linearly implicit RHS, with non-time-dependent Jacobian */
255 -- flag = ARKodeSetLinear(arkode_mem, 0);
256 -- if (check_flag(&flag, "ARKodeSetLinear", 1)) return 1;
257
258 -- /* Open output stream for results, output comment line */
259 -- UFID = fopen("solution.txt","w");
260 -- fprintf(UFID,"# t u\n");
261
262 -- /* output initial condition to disk */
263 -- fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", T0, NV_Ith_S(y,0));
264
265 -- /* Main time-stepping loop: calls ARKode to perform the integration, then
266 -- prints results. Stops when the final time has been reached */
267 -- t = T0;
268 -- tout = T0+dTout;
269 -- printf(" t u\n");
270 -- printf(" ---------------------\n");
271 -- while (Tf - t > 1.0e-15) {
272
273 -- flag = ARKode(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */
274 -- if (check_flag(&flag, "ARKode", 1)) break;
275 -- printf(" %10.6"FSYM" %10.6"FSYM"\n", t, NV_Ith_S(y,0)); /* access/print solution */
276 -- fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", t, NV_Ith_S(y,0));
277 -- if (flag >= 0) { /* successful solve: update time */
278 -- tout += dTout;
279 -- tout = (tout > Tf) ? Tf : tout;
280 -- } else { /* unsuccessful solve: break */
281 -- fprintf(stderr,"Solver failure, stopping integration\n");
282 -- break;
283 -- }
284 -- }
285 -- printf(" ---------------------\n");
286 -- fclose(UFID);
287
288 -- /* Get/print some final statistics on how the solve progressed */
289 -- flag = ARKodeGetNumSteps(arkode_mem, &nst);
290 -- check_flag(&flag, "ARKodeGetNumSteps", 1);
291 -- flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
292 -- check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
293 -- flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
294 -- check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
295 -- flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
296 -- check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
297 -- flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
298 -- check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
299 -- flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
300 -- check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
301 -- flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
302 -- check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
303 -- flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
304 -- check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
305 -- flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
306 -- check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
307
308 -- printf("\nFinal Solver Statistics:\n");
309 -- printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a);
310 -- printf(" Total RHS evals: Fe = %li, Fi = %li\n", nfe, nfi);
311 -- printf(" Total linear solver setups = %li\n", nsetups);
312 -- printf(" Total RHS evals for setting up the linear system = %li\n", nfeLS);
313 -- printf(" Total number of Jacobian evaluations = %li\n", nje);
314 -- printf(" Total number of Newton iterations = %li\n", nni);
315 -- printf(" Total number of linear solver convergence failures = %li\n", ncfn);
316 -- printf(" Total number of error test failures = %li\n\n", netf);
317
318 -- /* check the solution error */
319 -- flag = check_ans(y, t, reltol, abstol);
320
321 -- /* Clean up and return */
322 -- N_VDestroy(y); /* Free y vector */
323 -- ARKodeFree(&arkode_mem); /* Free integrator memory */
324 -- SUNLinSolFree(LS); /* Free linear solver */
325 -- SUNMatDestroy(A); /* Free A matrix */
326
327 -- return flag;
328 -- } |]
329 putStrLn $ show res 189 putStrLn $ show res
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-29 02:21:01 +0000