1 /* -----------------------------------------------------------------
2 * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, George D. Byrne,
3 * and Radu Serban @ LLNL
4 * -----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2021, Lawrence Livermore National Security
7 * and Southern Methodist University.
8 * All rights reserved.
9 *
10 * See the top-level LICENSE and NOTICE files for details.
11 *
12 * SPDX-License-Identifier: BSD-3-Clause
13 * SUNDIALS Copyright End
14 * -----------------------------------------------------------------
15 * Example problem:
16 *
17 * The following is a simple example problem, with the program for
18 * its solution by CVODES. The problem is the semi-discrete form of
19 * the advection-diffusion equation in 1-D:
20 * du/dt = q1 * d^2 u / dx^2 + q2 * du/dx
21 * on the interval 0 <= x <= 2, and the time interval 0 <= t <= 5.
22 * Homogeneous Dirichlet boundary conditions are posed, and the
23 * initial condition is:
24 * u(x,y,t=0) = x(2-x)exp(2x).
25 * The PDE is discretized on a uniform grid of size MX+2 with
26 * central differencing, and with boundary values eliminated,
27 * leaving an ODE system of size NEQ = MX.
28 * This program solves the problem with the option for nonstiff
29 * systems: ADAMS method and functional iteration.
30 * It uses scalar relative and absolute tolerances.
31 * Output is printed at t = .5, 1.0, ..., 5.
32 * Run statistics (optional outputs) are printed at the end.
33 *
34 * Optionally, CVODES can compute sensitivities with respect to the
35 * problem parameters q1 and q2.
36 * Any of three sensitivity methods (SIMULTANEOUS, STAGGERED, and
37 * STAGGERED1) can be used and sensitivities may be included in the
38 * error test or not (error control set on FULL or PARTIAL,
39 * respectively).
40 *
41 * Execution:
42 *
43 * If no sensitivities are desired:
44 * % cvsAdvDiff_FSA_non -nosensi
45 * If sensitivities are to be computed:
46 * % cvsAdvDiff_FSA_non -sensi sensi_meth err_con
47 * where sensi_meth is one of {sim, stg, stg1} and err_con is one of
48 * {t, f}.
49 * -----------------------------------------------------------------*/
50
51 #include <stdio.h>
52 #include <stdlib.h>
53 #include <string.h>
54 #include <math.h>
55
56 #include <cvodes/cvodes.h>
57 #include <nvector/nvector_serial.h>
58 #include <sundials/sundials_types.h>
59 #include <sundials/sundials_math.h>
60 #include "sunnonlinsol/sunnonlinsol_fixedpoint.h" /* access to the fixed point SUNNonlinearSolver */
61
62 /* Problem Constants */
63 #define XMAX RCONST(2.0) /* domain boundary */
64 #define MX 10 /* mesh dimension */
65 #define NEQ MX /* number of equations */
66 #define ATOL RCONST(1.e-5) /* scalar absolute tolerance */
67 #define T0 RCONST(0.0) /* initial time */
68 #define T1 RCONST(0.5) /* first output time */
69 #define DTOUT RCONST(0.5) /* output time increment */
70 #define NOUT 10 /* number of output times */
71
72 #define NP 2
73 #define NS 2
74
75 #define ZERO RCONST(0.0)
76
77 /* Type : UserData
78 contains problem parameters, grid constants, work array. */
79
80 typedef struct {
81 realtype *p;
82 realtype dx;
83 } *UserData;
84
85 /* Functions Called by the CVODES Solver */
86
87 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data);
88
89 /* Private Helper Functions */
90
91 static void ProcessArgs(int argc, char *argv[],
92 booleantype *sensi, int *sensi_meth,
93 booleantype *err_con);
94 static void WrongArgs(char *name);
95 static void SetIC(N_Vector u, realtype dx);
96 static void PrintOutput(void *cvode_mem, realtype t, N_Vector u);
97 static void PrintOutputS(N_Vector *uS);
98 static void PrintFinalStats(void *cvode_mem, booleantype sensi,
99 booleantype err_con, int sensi_meth);
100
101 static int check_retval(void *returnvalue, const char *funcname, int opt);
102
103 /*
104 *--------------------------------------------------------------------
105 * MAIN PROGRAM
106 *--------------------------------------------------------------------
107 */
108
main(int argc,char * argv[])109 int main(int argc, char *argv[])
110 {
111 void *cvode_mem;
112 UserData data;
113 realtype dx, reltol, abstol, t, tout;
114 N_Vector u;
115 int iout, retval;
116
117 realtype *pbar;
118 int is, *plist;
119 N_Vector *uS;
120 booleantype sensi, err_con;
121 int sensi_meth;
122
123 SUNNonlinearSolver NLS, NLSsens;
124
125 cvode_mem = NULL;
126 data = NULL;
127 u = NULL;
128 pbar = NULL;
129 plist = NULL;
130 uS = NULL;
131 NLS = NULL;
132 NLSsens = NULL;
133
134 /* Process arguments */
135 ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);
136
137 /* Set user data */
138 data = (UserData) malloc(sizeof *data); /* Allocate data memory */
139 if(check_retval((void *)data, "malloc", 2)) return(1);
140 data->p = (realtype *) malloc(NP * sizeof(realtype));
141 dx = data->dx = XMAX/((realtype)(MX+1));
142 data->p[0] = RCONST(1.0);
143 data->p[1] = RCONST(0.5);
144
145 /* Allocate and set initial states */
146 u = N_VNew_Serial(NEQ);
147 if(check_retval((void *)u, "N_VNew_Serial", 0)) return(1);
148 SetIC(u, dx);
149
150 /* Set integration tolerances */
151 reltol = ZERO;
152 abstol = ATOL;
153
154 /* Create CVODES object */
155 cvode_mem = CVodeCreate(CV_ADAMS);
156 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
157
158 retval = CVodeSetUserData(cvode_mem, data);
159 if(check_retval(&retval, "CVodeSetUserData", 1)) return(1);
160
161 /* Allocate CVODES memory */
162 retval = CVodeInit(cvode_mem, f, T0, u);
163 if(check_retval(&retval, "CVodeInit", 1)) return(1);
164
165 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
166 if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);
167
168 /* create fixed point nonlinear solver object */
169 NLS = SUNNonlinSol_FixedPoint(u, 0);
170 if(check_retval((void *)NLS, "SUNNonlinSol_FixedPoint", 0)) return(1);
171
172 /* attach nonlinear solver object to CVode */
173 retval = CVodeSetNonlinearSolver(cvode_mem, NLS);
174 if(check_retval(&retval, "CVodeSetNonlinearSolver", 1)) return(1);
175
176 printf("\n1-D advection-diffusion equation, mesh size =%3d\n", MX);
177
178 /* Sensitivity-related settings */
179 if(sensi) {
180
181 plist = (int *) malloc(NS * sizeof(int));
182 if(check_retval((void *)plist, "malloc", 2)) return(1);
183 for(is=0; is<NS; is++) plist[is] = is;
184
185 pbar = (realtype *) malloc(NS * sizeof(realtype));
186 if(check_retval((void *)pbar, "malloc", 2)) return(1);
187 for(is=0; is<NS; is++) pbar[is] = data->p[plist[is]];
188
189 uS = N_VCloneVectorArray(NS, u);
190 if(check_retval((void *)uS, "N_VCloneVectorArray", 0)) return(1);
191 for(is=0;is<NS;is++)
192 N_VConst(ZERO, uS[is]);
193
194 retval = CVodeSensInit1(cvode_mem, NS, sensi_meth, NULL, uS);
195 if(check_retval(&retval, "CVodeSensInit1", 1)) return(1);
196
197 retval = CVodeSensEEtolerances(cvode_mem);
198 if(check_retval(&retval, "CVodeSensEEtolerances", 1)) return(1);
199
200 retval = CVodeSetSensErrCon(cvode_mem, err_con);
201 if(check_retval(&retval, "CVodeSetSensErrCon", 1)) return(1);
202
203 retval = CVodeSetSensDQMethod(cvode_mem, CV_CENTERED, ZERO);
204 if(check_retval(&retval, "CVodeSetSensDQMethod", 1)) return(1);
205
206 retval = CVodeSetSensParams(cvode_mem, data->p, pbar, plist);
207 if(check_retval(&retval, "CVodeSetSensParams", 1)) return(1);
208
209 /* create sensitivity fixed point nonlinear solver object */
210 if (sensi_meth == CV_SIMULTANEOUS)
211 NLSsens = SUNNonlinSol_FixedPointSens(NS+1, u, 0);
212 else if(sensi_meth == CV_STAGGERED)
213 NLSsens = SUNNonlinSol_FixedPointSens(NS, u, 0);
214 else
215 NLSsens = SUNNonlinSol_FixedPoint(u, 0);
216 if(check_retval((void *)NLS, "SUNNonlinSol_FixedPoint", 0)) return(1);
217
218 /* attach nonlinear solver object to CVode */
219 if (sensi_meth == CV_SIMULTANEOUS)
220 retval = CVodeSetNonlinearSolverSensSim(cvode_mem, NLSsens);
221 else if(sensi_meth == CV_STAGGERED)
222 retval = CVodeSetNonlinearSolverSensStg(cvode_mem, NLSsens);
223 else
224 retval = CVodeSetNonlinearSolverSensStg1(cvode_mem, NLSsens);
225 if(check_retval(&retval, "CVodeSetNonlinearSolver", 1)) return(1);
226
227 printf("Sensitivity: YES ");
228 if(sensi_meth == CV_SIMULTANEOUS)
229 printf("( SIMULTANEOUS +");
230 else
231 if(sensi_meth == CV_STAGGERED) printf("( STAGGERED +");
232 else printf("( STAGGERED1 +");
233 if(err_con) printf(" FULL ERROR CONTROL )");
234 else printf(" PARTIAL ERROR CONTROL )");
235
236 } else {
237
238 printf("Sensitivity: NO ");
239
240 }
241
242 /* In loop over output points, call CVode, print results, test for error */
243
244 printf("\n\n");
245 printf("============================================================\n");
246 printf(" T Q H NST Max norm \n");
247 printf("============================================================\n");
248
249 for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {
250 retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
251 if(check_retval(&retval, "CVode", 1)) break;
252 PrintOutput(cvode_mem, t, u);
253 if (sensi) {
254 retval = CVodeGetSens(cvode_mem, &t, uS);
255 if(check_retval(&retval, "CVodeGetSens", 1)) break;
256 PrintOutputS(uS);
257 }
258 printf("------------------------------------------------------------\n");
259 }
260
261 /* Print final statistics */
262 PrintFinalStats(cvode_mem, sensi, err_con, sensi_meth);
263
264 /* Free memory */
265 N_VDestroy(u);
266 if (sensi) {
267 N_VDestroyVectorArray(uS, NS);
268 free(plist);
269 free(pbar);
270 }
271 free(data->p);
272 free(data);
273 CVodeFree(&cvode_mem);
274 SUNNonlinSolFree(NLS);
275 if (sensi) SUNNonlinSolFree(NLSsens);
276
277 return(0);
278 }
279
280 /*
281 *--------------------------------------------------------------------
282 * FUNCTIONS CALLED BY CVODES
283 *--------------------------------------------------------------------
284 */
285
286 /*
287 * f routine. Compute f(t,u).
288 */
289
f(realtype t,N_Vector u,N_Vector udot,void * user_data)290 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data)
291 {
292 realtype ui, ult, urt, hordc, horac, hdiff, hadv;
293 realtype dx;
294 realtype *udata, *dudata;
295 int i;
296 UserData data;
297
298 udata = N_VGetArrayPointer(u);
299 dudata = N_VGetArrayPointer(udot);
300
301 /* Extract needed problem constants from data */
302 data = (UserData) user_data;
303 dx = data->dx;
304 hordc = data->p[0]/(dx*dx);
305 horac = data->p[1]/(RCONST(2.0)*dx);
306
307 /* Loop over all grid points. */
308 for (i=0; i<NEQ; i++) {
309
310 /* Extract u at x_i and two neighboring points */
311 ui = udata[i];
312 if(i!=0)
313 ult = udata[i-1];
314 else
315 ult = ZERO;
316 if(i!=NEQ-1)
317 urt = udata[i+1];
318 else
319 urt = ZERO;
320
321 /* Set diffusion and advection terms and load into udot */
322 hdiff = hordc*(ult - RCONST(2.0)*ui + urt);
323 hadv = horac*(urt - ult);
324 dudata[i] = hdiff + hadv;
325 }
326
327 return(0);
328 }
329
330 /*
331 *--------------------------------------------------------------------
332 * PRIVATE FUNCTIONS
333 *--------------------------------------------------------------------
334 */
335
336 /*
337 * Process and verify arguments.
338 */
339
ProcessArgs(int argc,char * argv[],booleantype * sensi,int * sensi_meth,booleantype * err_con)340 static void ProcessArgs(int argc, char *argv[],
341 booleantype *sensi, int *sensi_meth, booleantype *err_con)
342 {
343 *sensi = SUNFALSE;
344 *sensi_meth = -1;
345 *err_con = SUNFALSE;
346
347 if (argc < 2) WrongArgs(argv[0]);
348
349 if (strcmp(argv[1],"-nosensi") == 0)
350 *sensi = SUNFALSE;
351 else if (strcmp(argv[1],"-sensi") == 0)
352 *sensi = SUNTRUE;
353 else
354 WrongArgs(argv[0]);
355
356 if (*sensi) {
357
358 if (argc != 4)
359 WrongArgs(argv[0]);
360
361 if (strcmp(argv[2],"sim") == 0)
362 *sensi_meth = CV_SIMULTANEOUS;
363 else if (strcmp(argv[2],"stg") == 0)
364 *sensi_meth = CV_STAGGERED;
365 else if (strcmp(argv[2],"stg1") == 0)
366 *sensi_meth = CV_STAGGERED1;
367 else
368 WrongArgs(argv[0]);
369
370 if (strcmp(argv[3],"t") == 0)
371 *err_con = SUNTRUE;
372 else if (strcmp(argv[3],"f") == 0)
373 *err_con = SUNFALSE;
374 else
375 WrongArgs(argv[0]);
376 }
377
378 }
379
WrongArgs(char * name)380 static void WrongArgs(char *name)
381 {
382 printf("\nUsage: %s [-nosensi] [-sensi sensi_meth err_con]\n",name);
383 printf(" sensi_meth = sim, stg, or stg1\n");
384 printf(" err_con = t or f\n");
385
386 exit(0);
387 }
388
389 /*
390 * Set initial conditions in u vector.
391 */
392
SetIC(N_Vector u,realtype dx)393 static void SetIC(N_Vector u, realtype dx)
394 {
395 int i;
396 realtype x;
397 realtype *udata;
398
399 /* Set pointer to data array and get local length of u. */
400 udata = N_VGetArrayPointer(u);
401
402 /* Load initial profile into u vector */
403 for (i=0; i<NEQ; i++) {
404 x = (i+1)*dx;
405 udata[i] = x*(XMAX - x)*SUNRexp(RCONST(2.0)*x);
406 }
407 }
408
409 /*
410 * Print current t, step count, order, stepsize, and max norm of solution
411 */
412
PrintOutput(void * cvode_mem,realtype t,N_Vector u)413 static void PrintOutput(void *cvode_mem, realtype t, N_Vector u)
414 {
415 long int nst;
416 int qu, retval;
417 realtype hu;
418
419 retval = CVodeGetNumSteps(cvode_mem, &nst);
420 check_retval(&retval, "CVodeGetNumSteps", 1);
421 retval = CVodeGetLastOrder(cvode_mem, &qu);
422 check_retval(&retval, "CVodeGetLastOrder", 1);
423 retval = CVodeGetLastStep(cvode_mem, &hu);
424 check_retval(&retval, "CVodeGetLastStep", 1);
425
426 #if defined(SUNDIALS_EXTENDED_PRECISION)
427 printf("%8.3Le %2d %8.3Le %5ld\n", t, qu, hu ,nst);
428 #elif defined(SUNDIALS_DOUBLE_PRECISION)
429 printf("%8.3e %2d %8.3e %5ld\n", t, qu, hu ,nst);
430 #else
431 printf("%8.3e %2d %8.3e %5ld\n", t, qu, hu ,nst);
432 #endif
433
434 printf(" Solution ");
435
436 #if defined(SUNDIALS_EXTENDED_PRECISION)
437 printf("%12.4Le \n", N_VMaxNorm(u));
438 #elif defined(SUNDIALS_DOUBLE_PRECISION)
439 printf("%12.4e \n", N_VMaxNorm(u));
440 #else
441 printf("%12.4e \n", N_VMaxNorm(u));
442 #endif
443 }
444
445 /*
446 * Print max norm of sensitivities
447 */
448
PrintOutputS(N_Vector * uS)449 static void PrintOutputS(N_Vector *uS)
450 {
451 printf(" Sensitivity 1 ");
452 #if defined(SUNDIALS_EXTENDED_PRECISION)
453 printf("%12.4Le \n", N_VMaxNorm(uS[0]));
454 #elif defined(SUNDIALS_DOUBLE_PRECISION)
455 printf("%12.4e \n", N_VMaxNorm(uS[0]));
456 #else
457 printf("%12.4e \n", N_VMaxNorm(uS[0]));
458 #endif
459
460 printf(" Sensitivity 2 ");
461 #if defined(SUNDIALS_EXTENDED_PRECISION)
462 printf("%12.4Le \n", N_VMaxNorm(uS[1]));
463 #elif defined(SUNDIALS_DOUBLE_PRECISION)
464 printf("%12.4e \n", N_VMaxNorm(uS[1]));
465 #else
466 printf("%12.4e \n", N_VMaxNorm(uS[1]));
467 #endif
468 }
469
470
471 /*
472 * Print some final statistics located in the CVODES memory
473 */
474
PrintFinalStats(void * cvode_mem,booleantype sensi,booleantype err_con,int sensi_meth)475 static void PrintFinalStats(void *cvode_mem, booleantype sensi,
476 booleantype err_con, int sensi_meth)
477 {
478 long int nst;
479 long int nfe, nsetups, nni, ncfn, netf;
480 long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;
481 int retval;
482
483 retval = CVodeGetNumSteps(cvode_mem, &nst);
484 check_retval(&retval, "CVodeGetNumSteps", 1);
485 retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
486 check_retval(&retval, "CVodeGetNumRhsEvals", 1);
487 retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
488 check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
489 retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
490 check_retval(&retval, "CVodeGetNumErrTestFails", 1);
491 retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
492 check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
493 retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
494 check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
495
496 if (sensi) {
497 retval = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);
498 check_retval(&retval, "CVodeGetSensNumRhsEvals", 1);
499 retval = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);
500 check_retval(&retval, "CVodeGetNumRhsEvalsSens", 1);
501 retval = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);
502 check_retval(&retval, "CVodeGetSensNumLinSolvSetups", 1);
503 if (err_con) {
504 retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);
505 check_retval(&retval, "CVodeGetSensNumErrTestFails", 1);
506 } else {
507 netfS = 0;
508 }
509 if ((sensi_meth == CV_STAGGERED) || (sensi_meth == CV_STAGGERED1)) {
510 retval = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);
511 check_retval(&retval, "CVodeGetSensNumNonlinSolvIters", 1);
512 retval = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);
513 check_retval(&retval, "CVodeGetSensNumNonlinSolvConvFails", 1);
514 } else {
515 nniS = 0;
516 ncfnS = 0;
517 }
518 }
519
520 printf("\nFinal Statistics\n\n");
521 printf("nst = %5ld\n\n", nst);
522 printf("nfe = %5ld\n", nfe);
523 printf("netf = %5ld nsetups = %5ld\n", netf, nsetups);
524 printf("nni = %5ld ncfn = %5ld\n", nni, ncfn);
525
526 if(sensi) {
527 printf("\n");
528 printf("nfSe = %5ld nfeS = %5ld\n", nfSe, nfeS);
529 printf("netfs = %5ld nsetupsS = %5ld\n", netfS, nsetupsS);
530 printf("nniS = %5ld ncfnS = %5ld\n", nniS, ncfnS);
531 }
532
533 }
534
535 /*
536 * Check function return value...
537 * opt == 0 means SUNDIALS function allocates memory so check if
538 * returned NULL pointer
539 * opt == 1 means SUNDIALS function returns an integer value so check if
540 * retval < 0
541 * opt == 2 means function allocates memory so check if returned
542 * NULL pointer
543 */
544
check_retval(void * returnvalue,const char * funcname,int opt)545 static int check_retval(void *returnvalue, const char *funcname, int opt)
546 {
547 int *retval;
548
549 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
550 if (opt == 0 && returnvalue == NULL) {
551 fprintf(stderr,
552 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
553 funcname);
554 return(1); }
555
556 /* Check if retval < 0 */
557 else if (opt == 1) {
558 retval = (int *) returnvalue;
559 if (*retval < 0) {
560 fprintf(stderr,
561 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
562 funcname, *retval);
563 return(1); }}
564
565 /* Check if function returned NULL pointer - no memory allocated */
566 else if (opt == 2 && returnvalue == NULL) {
567 fprintf(stderr,
568 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
569 funcname);
570 return(1); }
571
572 return(0);
573 }
574