summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-03-13 15:37:13 +0000
committerDominic Steinitz <dominic@steinitz.org>2018-03-13 15:37:13 +0000
commitff9f12bc128010b555ccc45280a576b560e33571 (patch)
treec0646bba5db8c40af7a4763ee9ea902f2914e93d
parentf2b1eae3d71c546abc71e099b4bd86010627f0fb (diff)
Almost all the C example
-rw-r--r--packages/sundials/src/Main.hs52
1 files changed, 52 insertions, 0 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 9e8bc63..7388b2b 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -189,6 +189,58 @@ main = do
189 flag = ARKodeInit(arkode_mem, NULL, f, T0, y); 189 flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
190 if (check_flag(&flag, "ARKodeInit", 1)) return 1; 190 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
191 191
192 /* Set routines */
193 flag = ARKodeSetUserData(arkode_mem, (void *) &lamda); /* Pass lamda to user functions */
194 if (check_flag(&flag, "ARKodeSetUserData", 1)) return 1;
195 flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
196 if (check_flag(&flag, "ARKodeSStolerances", 1)) return 1;
197
198 /* Initialize dense matrix data structure and solver */
199 A = SUNDenseMatrix(NEQ, NEQ);
200 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
201 LS = SUNDenseLinearSolver(y, A);
202 if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
203
204 /* Linear solver interface */
205 flag = ARKDlsSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */
206 if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
207 flag = ARKDlsSetJacFn(arkode_mem, Jac); /* Set Jacobian routine */
208 if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
209
210 /* Specify linearly implicit RHS, with non-time-dependent Jacobian */
211 flag = ARKodeSetLinear(arkode_mem, 0);
212 if (check_flag(&flag, "ARKodeSetLinear", 1)) return 1;
213
214 /* Open output stream for results, output comment line */
215 UFID = fopen("solution.txt","w");
216 fprintf(UFID,"# t u\n");
217
218 /* output initial condition to disk */
219 fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", T0, NV_Ith_S(y,0));
220
221 /* Main time-stepping loop: calls ARKode to perform the integration, then
222 prints results. Stops when the final time has been reached */
223 t = T0;
224 tout = T0+dTout;
225 printf(" t u\n");
226 printf(" ---------------------\n");
227 while (Tf - t > 1.0e-15) {
228
229 flag = ARKode(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */
230 if (check_flag(&flag, "ARKode", 1)) break;
231 printf(" %10.6"FSYM" %10.6"FSYM"\n", t, NV_Ith_S(y,0)); /* access/print solution */
232 fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", t, NV_Ith_S(y,0));
233 if (flag >= 0) { /* successful solve: update time */
234 tout += dTout;
235 tout = (tout > Tf) ? Tf : tout;
236 } else { /* unsuccessful solve: break */
237 fprintf(stderr,"Solver failure, stopping integration\n");
238 break;
239 }
240 }
241 printf(" ---------------------\n");
242 fclose(UFID);
243
192 return 0; 244 return 0;
193 } |] 245 } |]
194 putStrLn $ show res 246 putStrLn $ show res