1 /*---------------------------------------------------------------
2 * Programmer(s): Shelby Lockhart @ LLNL
3 *---------------------------------------------------------------
4 * Based on the serial code ark_heat1D_adapt.c developed
5 * by Daniel R. Reynolds and parallelized with OpenMP 4.5
6 *---------------------------------------------------------------
7 * SUNDIALS Copyright Start
8 * Copyright (c) 2002-2021, Lawrence Livermore National Security
9 * and Southern Methodist University.
10 * All rights reserved.
11 *
12 * See the top-level LICENSE and NOTICE files for details.
13 *
14 * SPDX-License-Identifier: BSD-3-Clause
15 * SUNDIALS Copyright End
16 *---------------------------------------------------------------
17 * Example problem:
18 *
19 * The following test simulates a simple 1D heat equation,
20 * u_t = k*u_xx + f
21 * for t in [0, 10], x in [0, 1], with initial conditions
22 * u(0,x) = 0
23 * Dirichlet boundary conditions, i.e.
24 * u_t(t,0) = u_t(t,1) = 0,
25 * and a heating term of the form
26 * f = 2*exp(-200*(x-0.25)*(x-0.25))
27 * - exp(-400*(x-0.7)*(x-0.7))
28 * + exp(-500*(x-0.4)*(x-0.4))
29 * - 2*exp(-600*(x-0.55)*(x-0.55));
30 *
31 * The spatial derivatives are computed using a three-point
32 * centered stencil (second order for a uniform mesh). The data
33 * is initially uniformly distributed over N points in the interval
34 * [0, 1], but as the simulation proceeds the mesh is adapted.
35 *
36 * This program solves the problem with a DIRK method, solved with
37 * a Newton iteration and SUNLinSol_PCG linear solver, with a
38 * user-supplied Jacobian-vector product routine.
39 *---------------------------------------------------------------*/
40
41 /* Header files */
42 #include <stdio.h>
43 #include <stdlib.h>
44 #include <math.h>
45 #include <arkode/arkode_arkstep.h> /* prototypes for ARKStep fcts., consts */
46 #include <nvector/nvector_openmpdev.h> /* OpenMPDEV N_Vector types, fcts., macros */
47 #include <sunlinsol/sunlinsol_pcg.h> /* access to PCG SUNLinearSolver */
48 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype, etc */
49 #include <sundials/sundials_math.h> /* def. of SUNRsqrt, etc. */
50
51 #ifdef _OPENMP
52 #include <omp.h> /* OpenMP functions */
53 #endif
54
55 #if defined(SUNDIALS_EXTENDED_PRECISION)
56 #define GSYM "Lg"
57 #define ESYM "Le"
58 #define FSYM "Lf"
59 #else
60 #define GSYM "g"
61 #define ESYM "e"
62 #define FSYM "f"
63 #endif
64
65 /* constants */
66 #define ZERO RCONST(0.0)
67 #define PT25 RCONST(0.25)
68 #define PT4 RCONST(0.4)
69 #define PT5 RCONST(0.5)
70 #define PT55 RCONST(0.55)
71 #define PT7 RCONST(0.7)
72 #define ONE RCONST(1.0)
73 #define TWO RCONST(2.0)
74 #define TWOHUNDRED RCONST(200.0)
75 #define FOURHUNDRED RCONST(400.0)
76 #define FIVEHUNDRED RCONST(500.0)
77 #define SIXHUNDRED RCONST(600.0)
78
79 /* user data structure */
80 typedef struct {
81 sunindextype N; /* current number of intervals */
82 realtype *x_host; /* current mesh on host */
83 realtype *x_dev; /* current mesh on device */
84 realtype k; /* diffusion coefficient */
85 realtype refine_tol; /* adaptivity tolerance */
86 } *UserData;
87
88 /* User-supplied Functions Called by the Solver */
89 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
90 static int Jac(N_Vector v, N_Vector Jv, realtype t, N_Vector y,
91 N_Vector fy, void *user_data, N_Vector tmp);
92
93 /* Private function to check function return values */
94 realtype * adapt_mesh(N_Vector y, sunindextype *Nnew, UserData udata);
95 static int project(sunindextype Nold, realtype *xold, N_Vector yold,
96 sunindextype Nnew, realtype *xnew, N_Vector ynew);
97 static int check_flag(void *flagvalue, const char *funcname, int opt);
98
99 /* Main Program */
main()100 int main() {
101
102 /* general problem parameters */
103 realtype T0 = RCONST(0.0); /* initial time */
104 realtype Tf = RCONST(1.0); /* final time */
105 realtype rtol = RCONST(1.e-3); /* relative tolerance */
106 realtype atol = RCONST(1.e-10); /* absolute tolerance */
107 realtype hscale = RCONST(1.0); /* time step change factor on resizes */
108 UserData udata = NULL;
109 realtype *data;
110 sunindextype N = 21; /* initial spatial mesh size */
111 realtype refine = RCONST(3.0e-3); /* adaptivity refinement tolerance */
112 realtype k = RCONST(0.5); /* heat conductivity */
113 sunindextype i;
114 long int nni, nni_tot=0, nli, nli_tot=0;
115 int iout=0;
116
117 /* general problem variables */
118 int flag; /* reusable error-checking flag */
119 N_Vector y = NULL; /* empty vector for storing solution */
120 N_Vector y2 = NULL; /* empty vector for storing solution */
121 N_Vector yt = NULL; /* empty vector for swapping */
122 SUNLinearSolver LS = NULL; /* empty linear solver object */
123 void *arkode_mem = NULL; /* empty ARKode memory structure */
124 FILE *XFID, *UFID;
125 realtype t, olddt, newdt;
126 realtype *xnew_host = NULL;
127 realtype *xnew_dev = NULL;
128 sunindextype Nnew;
129 int dev, host;
130
131 /* get host and offloading device */
132 dev = omp_get_default_device();
133 host = omp_get_initial_device();
134
135 /* allocate and fill initial udata structure */
136 udata = (UserData) malloc(sizeof(*udata));
137 udata->N = N;
138 udata->k = k;
139 udata->refine_tol = refine;
140 udata->x_host = malloc(N * sizeof(realtype));
141 for (i=0; i<N; i++) udata->x_host[i] = ONE*i/(N-1);
142 udata->x_dev = omp_target_alloc(N * sizeof(realtype), dev);
143 omp_target_memcpy(udata->x_dev, udata->x_host, N * sizeof(realtype), 0, 0, dev, host);
144
145 /* Initial problem output */
146 printf("\n1D adaptive Heat PDE test problem:\n");
147 printf(" diffusion coefficient: k = %"GSYM"\n", udata->k);
148 printf(" initial N = %li\n", (long int) udata->N);
149
150 /* Initialize data structures */
151 y = N_VNew_OpenMPDEV(N); /* Create initial OpenMPDEV vector for solution */
152 if (check_flag((void *) y, "N_VNew_OpenMPDEV", 0)) return 1;
153 N_VConst(ZERO, y); /* Set initial conditions */
154
155 /* output mesh to disk */
156 XFID=fopen("heat_mesh.txt","w");
157
158 /* output initial mesh to disk */
159 for (i=0; i<udata->N; i++) fprintf(XFID," %.16"ESYM, udata->x_host[i]);
160 fprintf(XFID,"\n");
161
162 /* Open output stream for results, access data array */
163 UFID=fopen("heat1D.txt","w");
164
165 /* output initial condition to disk */
166 N_VCopyFromDevice_OpenMPDEV(y);
167 data = N_VGetHostArrayPointer_OpenMPDEV(y);
168 for (i=0; i<udata->N; i++) fprintf(UFID," %.16"ESYM, data[i]);
169 fprintf(UFID,"\n");
170
171 /* Initialize the ARK timestepper */
172 arkode_mem = ARKStepCreate(NULL, f, T0, y);
173 if (check_flag((void *) arkode_mem, "ARKStepCreate", 0)) return 1;
174
175 /* Set routines */
176 flag = ARKStepSetUserData(arkode_mem, (void *) udata); /* Pass udata to user functions */
177 if (check_flag(&flag, "ARKStepSetUserData", 1)) return 1;
178 flag = ARKStepSetMaxNumSteps(arkode_mem, 10000); /* Increase max num steps */
179 if (check_flag(&flag, "ARKStepSetMaxNumSteps", 1)) return 1;
180 flag = ARKStepSStolerances(arkode_mem, rtol, atol); /* Specify tolerances */
181 if (check_flag(&flag, "ARKStepSStolerances", 1)) return 1;
182 flag = ARKStepSetAdaptivityMethod(arkode_mem, 2, 1, 0, NULL); /* Set adaptivity method */
183 if (check_flag(&flag, "ARKStepSetAdaptivityMethod", 1)) return 1;
184 flag = ARKStepSetPredictorMethod(arkode_mem, 0); /* Set predictor method */
185 if (check_flag(&flag, "ARKStepSetPredictorMethod", 1)) return 1;
186
187 /* Specify linearly implicit RHS, with time-dependent Jacobian */
188 flag = ARKStepSetLinear(arkode_mem, 1);
189 if (check_flag(&flag, "ARKStepSetLinear", 1)) return 1;
190
191 /* Initialize PCG solver -- no preconditioning, with up to N iterations */
192 LS = SUNLinSol_PCG(y, 0, (int) N);
193 if (check_flag((void *)LS, "SUNLinSol_PCG", 0)) return 1;
194
195 /* Linear solver interface -- set user-supplied J*v routine (no 'jtsetup' required) */
196 flag = ARKStepSetLinearSolver(arkode_mem, LS, NULL); /* Attach linear solver to ARKStep */
197 if (check_flag(&flag, "ARKStepSetLinearSolver", 1)) return 1;
198 flag = ARKStepSetJacTimes(arkode_mem, NULL, Jac); /* Set the Jacobian routine */
199 if (check_flag(&flag, "ARKStepSetJacTimes", 1)) return 1;
200
201 /* Main time-stepping loop: calls ARKStepEvolve to perform the integration, then
202 prints results. Stops when the final time has been reached */
203 t = T0;
204 olddt = ZERO;
205 newdt = ZERO;
206 printf(" iout dt_old dt_new ||u||_rms N NNI NLI\n");
207 printf(" ----------------------------------------------------------------------------------------\n");
208 printf(" %4i %19.15"ESYM" %19.15"ESYM" %19.15"ESYM" %li %2i %3i\n",
209 iout, olddt, newdt, SUNRsqrt(N_VDotProd(y,y)/udata->N),
210 (long int) udata->N, 0, 0);
211 while (t < Tf) {
212
213 /* "set" routines */
214 flag = ARKStepSetStopTime(arkode_mem, Tf);
215 if (check_flag(&flag, "ARKStepSetStopTime", 1)) return 1;
216 flag = ARKStepSetInitStep(arkode_mem, newdt);
217 if (check_flag(&flag, "ARKStepSetInitStep", 1)) return 1;
218
219 /* call integrator */
220 flag = ARKStepEvolve(arkode_mem, Tf, y, &t, ARK_ONE_STEP);
221 if (check_flag(&flag, "ARKStepEvolve", 1)) return 1;
222
223 /* "get" routines */
224 flag = ARKStepGetLastStep(arkode_mem, &olddt);
225 if (check_flag(&flag, "ARKStepGetLastStep", 1)) return 1;
226 flag = ARKStepGetCurrentStep(arkode_mem, &newdt);
227 if (check_flag(&flag, "ARKStepGetCurrentStep", 1)) return 1;
228 flag = ARKStepGetNumNonlinSolvIters(arkode_mem, &nni);
229 if (check_flag(&flag, "ARKStepGetNumNonlinSolvIters", 1)) return 1;
230 flag = ARKStepGetNumLinIters(arkode_mem, &nli);
231 if (check_flag(&flag, "ARKStepGetNumLinIters", 1)) return 1;
232
233 /* print current solution stats */
234 iout++;
235 printf(" %4i %19.15"ESYM" %19.15"ESYM" %19.15"ESYM" %li %2li %3li\n",
236 iout, olddt, newdt, SUNRsqrt(N_VDotProd(y,y)/udata->N),
237 (long int) udata->N, nni, nli);
238 nni_tot += nni;
239 nli_tot += nli;
240
241 /* output results and current mesh to disk */
242 N_VCopyFromDevice_OpenMPDEV(y);
243 data = N_VGetHostArrayPointer_OpenMPDEV(y);
244 for (i=0; i<udata->N; i++) fprintf(UFID," %.16"ESYM, data[i]);
245 fprintf(UFID,"\n");
246 for (i=0; i<udata->N; i++) fprintf(XFID," %.16"ESYM, udata->x_host[i]);
247 fprintf(XFID,"\n");
248
249 /* adapt the spatial mesh */
250 xnew_host = adapt_mesh(y, &Nnew, udata);
251 if (check_flag(xnew_host, "ark_adapt", 0)) return 1;
252
253 /* create N_Vector of new length */
254 y2 = N_VNew_OpenMPDEV(Nnew);
255 if (check_flag((void *) y2, "N_VNew_OpenMPDEV", 0)) return 1;
256
257 /* copy new mesh from host array to device array */
258 xnew_dev = omp_target_alloc(Nnew * sizeof(realtype), dev);
259 omp_target_memcpy(xnew_dev, xnew_host, Nnew*sizeof(realtype), 0, 0, dev, host);
260
261 /* project solution onto new mesh */
262 flag = project(udata->N, udata->x_dev, y, Nnew, xnew_dev, y2);
263 if (check_flag(&flag, "project", 1)) return 1;
264
265 /* delete old vector, old mesh */
266 N_VDestroy(y);
267 free(udata->x_host);
268 omp_target_free(udata->x_dev, dev);
269
270 /* swap x and xnew so that new mesh is stored in udata structure */
271 udata->x_host = xnew_host;
272 xnew_host = NULL;
273 udata->N = Nnew; /* store size of new mesh */
274 udata->x_dev = xnew_dev;
275 xnew_dev = NULL;
276
277 /* swap y and y2 so that y holds new solution */
278 yt = y;
279 y = y2;
280 y2 = yt;
281
282 /* call ARKStepResize to notify integrator of change in mesh */
283 flag = ARKStepResize(arkode_mem, y, hscale, t, NULL, NULL);
284 if (check_flag(&flag, "ARKStepResize", 1)) return 1;
285
286 /* destroy and re-allocate linear solver memory; reattach to ARKStep interface */
287 SUNLinSolFree(LS);
288 LS = SUNLinSol_PCG(y, 0, (int) N);
289 if (check_flag((void *)LS, "SUNLinSol_PCG", 0)) return 1;
290 flag = ARKStepSetLinearSolver(arkode_mem, LS, NULL);
291 if (check_flag(&flag, "ARKStepSetLinearSolver", 1)) return 1;
292 flag = ARKStepSetJacTimes(arkode_mem, NULL, Jac);
293 if (check_flag(&flag, "ARKStepSetJacTimes", 1)) return 1;
294
295 }
296 printf(" ----------------------------------------------------------------------------------------\n");
297
298 /* print some final statistics */
299 printf(" Final solver statistics:\n");
300 printf(" Total number of time steps = %i\n", iout);
301 printf(" Total nonlinear iterations = %li\n", nni_tot);
302 printf(" Total linear iterations = %li\n\n", nli_tot);
303
304 /* Clean up and return with successful completion */
305 fclose(UFID);
306 fclose(XFID);
307 N_VDestroy(y); /* Free vectors */
308 free(udata->x_host); /* Free user data */
309 omp_target_free(udata->x_dev, dev);
310 free(udata);
311 ARKStepFree(&arkode_mem); /* Free integrator memory */
312 SUNLinSolFree(LS); /* Free linear solver */
313
314 return 0;
315 }
316
317 /*--------------------------------
318 * Functions called by the solver
319 *--------------------------------*/
320
321 /* f routine to compute the ODE RHS function f(t,y). */
f(realtype t,N_Vector y,N_Vector ydot,void * user_data)322 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
323 {
324 UserData udata = (UserData) user_data; /* access problem data */
325 sunindextype N = udata->N; /* set variable shortcuts */
326 realtype k = udata->k;
327 realtype *x = udata->x_dev;
328 realtype *Y=NULL, *Ydot=NULL;
329 realtype dxL, dxR;
330 sunindextype i;
331 int dev;
332
333 dev = omp_get_default_device();
334
335 /* access data arrays */
336 Y = N_VGetDeviceArrayPointer_OpenMPDEV(y);
337 if (check_flag((void *) Y, "N_VGetDeviceArrayPointer", 0)) return 1;
338 Ydot = N_VGetDeviceArrayPointer_OpenMPDEV(ydot);
339 if (check_flag((void *) Ydot, "N_VGetDeviceArrayPointer", 0)) return 1;
340
341 /* Initialize ydot to zero - also handles boundary conditions */
342 N_VConst(ZERO, ydot);
343
344 /* iterate over domain interior, computing all equations */
345 #pragma omp target map(to:N) is_device_ptr(x, Ydot, Y) device(dev)
346 #pragma omp teams distribute parallel for schedule(static, 1)
347 for (i=1; i<N-1; i++) { /* interior */
348 dxL = x[i]-x[i-1];
349 dxR = x[i+1]-x[i];
350 Ydot[i] = Y[i-1]*k*TWO/(dxL*(dxL+dxR))
351 - Y[i]*k*TWO/(dxL*dxR)
352 + Y[i+1]*k*TWO/(dxR*(dxL+dxR))
353 + TWO*SUNRexp(-TWOHUNDRED*(x[i]-PT25)*(x[i]-PT25)) /* source term */
354 - SUNRexp(-FOURHUNDRED*(x[i]-PT7)*(x[i]-PT7))
355 + SUNRexp(-FIVEHUNDRED*(x[i]-PT4)*(x[i]-PT4))
356 - TWO*SUNRexp(-SIXHUNDRED*(x[i]-PT55)*(x[i]-PT55));
357 }
358
359 return 0; /* Return with success */
360 }
361
362 /* Jacobian routine to compute J(t,y) = df/dy. */
Jac(N_Vector v,N_Vector Jv,realtype t,N_Vector y,N_Vector fy,void * user_data,N_Vector tmp)363 static int Jac(N_Vector v, N_Vector Jv, realtype t, N_Vector y,
364 N_Vector fy, void *user_data, N_Vector tmp)
365 {
366 UserData udata = (UserData) user_data; /* variable shortcuts */
367 sunindextype N = udata->N;
368 realtype k = udata->k;
369 realtype *x = udata->x_dev;
370 realtype *V=NULL, *JV=NULL;
371 realtype dxL, dxR;
372 sunindextype i;
373 int dev;
374
375 dev = omp_get_default_device();
376
377 /* access data arrays */
378 V = N_VGetDeviceArrayPointer_OpenMPDEV(v);
379 if (check_flag((void *) V, "N_VGetDeviceArrayPointer", 0)) return 1;
380 JV = N_VGetDeviceArrayPointer_OpenMPDEV(Jv);
381 if (check_flag((void *) JV, "N_VGetDeviceArrayPointer", 0)) return 1;
382
383 /* initialize Jv product to zero - also handles boundary conditions */
384 N_VConst(ZERO, Jv);
385
386 /* iterate over domain, computing all Jacobian-vector products */
387 #pragma omp target map(to:N) is_device_ptr(x, JV, V) device(dev)
388 #pragma omp teams distribute parallel for schedule(static, 1)
389 for (i=1; i<N-1; i++) {
390 dxL = x[i]-x[i-1];
391 dxR = x[i+1]-x[i];
392 JV[i] = V[i-1]*k*TWO/(dxL*(dxL+dxR))
393 - V[i]*k*TWO/(dxL*dxR)
394 + V[i+1]*k*TWO/(dxR*(dxL+dxR));
395 }
396
397 return 0; /* Return with success */
398 }
399
400 /*-------------------------------
401 * Private helper functions
402 *-------------------------------*/
403
404 /* Adapts the current mesh, using a simple adaptivity strategy of
405 refining when an approximation of the scaled second-derivative is
406 too large. We only do this in one sweep, so no attempt is made to
407 ensure the resulting mesh meets these same criteria after adaptivity:
408 y [input] -- the current solution vector
409 Nnew [output] -- the size of the new mesh
410 udata [input] -- the current system information
411 The return for this function is a pointer to the new mesh. */
adapt_mesh(N_Vector y,sunindextype * Nnew,UserData udata)412 realtype* adapt_mesh(N_Vector y, sunindextype *Nnew, UserData udata)
413 {
414 sunindextype i, j;
415 int *marks=NULL;
416 realtype ydd, *xold=NULL, *Y=NULL, *xnew=NULL;
417 sunindextype num_refine, N_new;
418
419 /* Access current solution and mesh arrays */
420 xold = udata->x_host;
421 Y = N_VGetHostArrayPointer_OpenMPDEV(y); /* assumes copy to host already done */
422 if (check_flag((void *) Y, "N_VGetHostArrayPointer_OpenMPDEV", 0)) return NULL;
423
424 /* create marking array */
425 marks = calloc(udata->N-1, sizeof(int));
426
427 /* perform marking:
428 0 -> leave alone
429 1 -> refine */
430 for (i=1; i<udata->N-1; i++) {
431
432 /* approximate scaled second-derivative */
433 ydd = Y[i-1] - TWO*Y[i] + Y[i+1];
434
435 /* check for refinement */
436 if (fabs(ydd) > udata->refine_tol) {
437 marks[i-1] = 1;
438 marks[i] = 1;
439 }
440
441 }
442
443 /* allocate new mesh */
444 num_refine = 0;
445 for (i=0; i<udata->N-1; i++)
446 if (marks[i] == 1) num_refine++;
447 N_new = udata->N + num_refine;
448 *Nnew = N_new; /* Store new array length */
449 xnew = malloc((N_new) * sizeof(realtype));
450
451
452 /* fill new mesh */
453 xnew[0] = xold[0]; /* store endpoints */
454 xnew[N_new-1] = xold[udata->N-1];
455 j=1;
456 /* iterate over old intervals */
457 for (i=0; i<udata->N-1; i++) {
458 /* if mark is 0, reuse old interval */
459 if (marks[i] == 0) {
460 xnew[j++] = xold[i+1];
461 continue;
462 }
463
464 /* if mark is 1, refine old interval */
465 if (marks[i] == 1) {
466 xnew[j++] = PT5*(xold[i]+xold[i+1]);
467 xnew[j++] = xold[i+1];
468 continue;
469 }
470 }
471
472 /* verify that new mesh is legal */
473 for (i=0; i<N_new-1; i++) {
474 if (xnew[i+1] <= xnew[i]) {
475 fprintf(stderr,"adapt_mesh error: illegal mesh created\n");
476 free(xnew);
477 return NULL;
478 }
479 }
480
481 free(marks); /* Delete marking array */
482 return xnew; /* Return with success */
483 }
484
485
486 /* Projects one vector onto another:
487 Nold [input] -- the size of the old mesh
488 xold [input] -- the old mesh
489 yold [input] -- the vector defined over the old mesh
490 Nnew [input] -- the size of the new mesh
491 xnew [input] -- the new mesh
492 ynew [output] -- the vector defined over the new mesh
493 (allocated prior to calling project) */
project(sunindextype Nold,realtype * xold,N_Vector yold,sunindextype Nnew,realtype * xnew,N_Vector ynew)494 static int project(sunindextype Nold, realtype *xold, N_Vector yold,
495 sunindextype Nnew, realtype *xnew, N_Vector ynew)
496 {
497 sunindextype iv, i, j;
498 realtype *Yold=NULL, *Ynew=NULL;
499
500 int dev = omp_get_default_device();
501
502 /* Access data arrays */
503 Yold = N_VGetDeviceArrayPointer_OpenMPDEV(yold); /* access data arrays */
504 if (check_flag((void *) Yold, "N_VGetDeviceArrayPointer_OpenMPDEV", 0)) return 1;
505 Ynew = N_VGetDeviceArrayPointer_OpenMPDEV(ynew);
506 if (check_flag((void *) Ynew, "N_VGetDeviceArrayPointer_OpenMPDEV", 0)) return 1;
507
508 /* loop over new mesh, finding corresponding interval within old mesh,
509 and perform piecewise linear interpolation from yold to ynew */
510 iv=0;
511 #pragma omp target map(to:iv) is_device_ptr(Yold,Ynew,xnew,xold) device(dev)
512 #pragma omp teams distribute parallel for schedule(static, 1)
513 {
514 for (i=0; i<Nnew; i++) {
515
516 /* find old interval, start with previous value since sorted */
517 for (j=iv; j<Nold-1; j++) {
518 if (xnew[i] >= xold[j] && xnew[i] <= xold[j+1]) {
519 iv = j;
520 break;
521 }
522 iv = Nold-1; /* just in case it wasn't found above */
523 }
524
525 /* perform interpolation */
526 Ynew[i] = Yold[iv]*(xnew[i]-xold[iv+1])/(xold[iv]-xold[iv+1])
527 + Yold[iv+1]*(xnew[i]-xold[iv])/(xold[iv+1]-xold[iv]);
528 }
529 }
530
531 return 0; /* Return with success */
532 }
533
534
535 /* Check function return value...
536 opt == 0 means SUNDIALS function allocates memory so check if
537 returned NULL pointer
538 opt == 1 means SUNDIALS function returns a flag so check if
539 flag >= 0
540 opt == 2 means function allocates memory so check if returned
541 NULL pointer
542 */
check_flag(void * flagvalue,const char * funcname,int opt)543 static int check_flag(void *flagvalue, const char *funcname, int opt)
544 {
545 int *errflag;
546
547 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
548 if (opt == 0 && flagvalue == NULL) {
549 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
550 funcname);
551 return 1; }
552
553 /* Check if flag < 0 */
554 else if (opt == 1) {
555 errflag = (int *) flagvalue;
556 if (*errflag < 0) {
557 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
558 funcname, *errflag);
559 return 1; }}
560
561 /* Check if function returned NULL pointer - no memory allocated */
562 else if (opt == 2 && flagvalue == NULL) {
563 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
564 funcname);
565 return 1; }
566
567 return 0;
568 }
569
570
571 /*---- end of file ----*/
572