1 /* -----------------------------------------------------------------
2 * Programmer(s): Allan Taylor, Alan Hindmarsh and
3 * Radu Serban @ LLNL
4 * -----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2020, 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 for IDA: 2D heat equation, serial, GMRES.
16 *
17 * This example solves a discretized 2D heat equation problem.
18 * This version uses the Krylov solver Spgmr.
19 *
20 * The DAE system solved is a spatial discretization of the PDE
21 * du/dt = d^2u/dx^2 + d^2u/dy^2
22 * on the unit square. The boundary condition is u = 0 on all edges.
23 * Initial conditions are given by u = 16 x (1 - x) y (1 - y). The
24 * PDE is treated with central differences on a uniform M x M grid.
25 * The values of u at the interior points satisfy ODEs, and
26 * equations u = 0 at the boundaries are appended, to form a DAE
27 * system of size N = M^2. Here M = 10.
28 *
29 * The system is solved with IDA using the Krylov linear solver
30 * SPGMR. The preconditioner uses the diagonal elements of the
31 * Jacobian only. Routines for preconditioning, required by
32 * SPGMR, are supplied here. The constraints u >= 0 are posed
33 * for all components. Output is taken at t = 0, .01, .02, .04,
34 * ..., 10.24. Two cases are run -- with the Gram-Schmidt type
35 * being Modified in the first case, and Classical in the second.
36 * The second run uses IDAReInit.
37 * -----------------------------------------------------------------*/
38
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <math.h>
42
43 #include <ida/ida.h> /* prototypes for IDA fcts., consts. */
44 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
45 #include <sunlinsol/sunlinsol_spgmr.h> /* access to spgmr SUNLinearSolver */
46 #include <sundials/sundials_types.h> /* definition of type realtype */
47
48 /* Problem Constants */
49
50 #define NOUT 11
51 #define MGRID 10
52 #define NEQ MGRID*MGRID
53 #define ZERO RCONST(0.0)
54 #define ONE RCONST(1.0)
55 #define TWO RCONST(2.0)
56 #define FOUR RCONST(4.0)
57
58 /* User data type */
59
60 typedef struct {
61 sunindextype mm; /* number of grid points */
62 realtype dx;
63 realtype coeff;
64 N_Vector pp; /* vector of prec. diag. elements */
65 } *UserData;
66
67 /* Prototypes for functions called by IDA */
68
69 int resHeat(realtype tres, N_Vector uu, N_Vector up,
70 N_Vector resval, void *user_data);
71
72 int PsetupHeat(realtype tt,
73 N_Vector uu, N_Vector up, N_Vector rr,
74 realtype c_j, void *prec_data);
75
76 int PsolveHeat(realtype tt,
77 N_Vector uu, N_Vector up, N_Vector rr,
78 N_Vector rvec, N_Vector zvec,
79 realtype c_j, realtype delta, void *prec_data);
80
81 /* Prototypes for private functions */
82
83 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
84 N_Vector res);
85 static void PrintHeader(realtype rtol, realtype atol);
86 static void PrintOutput(void *mem, realtype t, N_Vector uu);
87 static int check_retval(void *returnvalue, const char *funcname, int opt);
88
89 /*
90 *--------------------------------------------------------------------
91 * MAIN PROGRAM
92 *--------------------------------------------------------------------
93 */
94
main()95 int main()
96 {
97 void *mem;
98 UserData data;
99 N_Vector uu, up, constraints, res;
100 int retval, iout;
101 realtype rtol, atol, t0, t1, tout, tret;
102 long int netf, ncfn, ncfl;
103 SUNLinearSolver LS;
104
105 mem = NULL;
106 data = NULL;
107 uu = up = constraints = res = NULL;
108 LS = NULL;
109
110 /* Allocate N-vectors and the user data structure. */
111
112 uu = N_VNew_Serial(NEQ);
113 if(check_retval((void *)uu, "N_VNew_Serial", 0)) return(1);
114
115 up = N_VNew_Serial(NEQ);
116 if(check_retval((void *)up, "N_VNew_Serial", 0)) return(1);
117
118 res = N_VNew_Serial(NEQ);
119 if(check_retval((void *)res, "N_VNew_Serial", 0)) return(1);
120
121 constraints = N_VNew_Serial(NEQ);
122 if(check_retval((void *)constraints, "N_VNew_Serial", 0)) return(1);
123
124 data = (UserData) malloc(sizeof *data);
125 data->pp = NULL;
126 if(check_retval((void *)data, "malloc", 2)) return(1);
127
128 /* Assign parameters in the user data structure. */
129
130 data->mm = MGRID;
131 data->dx = ONE/(MGRID-ONE);
132 data->coeff = ONE/(data->dx * data->dx);
133 data->pp = N_VNew_Serial(NEQ);
134 if(check_retval((void *)data->pp, "N_VNew_Serial", 0)) return(1);
135
136 /* Initialize uu, up. */
137
138 SetInitialProfile(data, uu, up, res);
139
140 /* Set constraints to all 1's for nonnegative solution values. */
141
142 N_VConst(ONE, constraints);
143
144 /* Assign various parameters. */
145
146 t0 = ZERO;
147 t1 = RCONST(0.01);
148 rtol = ZERO;
149 atol = RCONST(1.0e-3);
150
151 /* Call IDACreate and IDAMalloc to initialize solution */
152
153 mem = IDACreate();
154 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
155
156 retval = IDASetUserData(mem, data);
157 if(check_retval(&retval, "IDASetUserData", 1)) return(1);
158
159 retval = IDASetConstraints(mem, constraints);
160 if(check_retval(&retval, "IDASetConstraints", 1)) return(1);
161 N_VDestroy(constraints);
162
163 retval = IDAInit(mem, resHeat, t0, uu, up);
164 if(check_retval(&retval, "IDAInit", 1)) return(1);
165
166 retval = IDASStolerances(mem, rtol, atol);
167 if(check_retval(&retval, "IDASStolerances", 1)) return(1);
168
169 /* Create the linear solver SUNLinSol_SPGMR with left preconditioning
170 and the default Krylov dimension */
171 LS = SUNLinSol_SPGMR(uu, PREC_LEFT, 0);
172 if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);
173
174 /* IDA recommends allowing up to 5 restarts (default is 0) */
175 retval = SUNLinSol_SPGMRSetMaxRestarts(LS, 5);
176 if(check_retval(&retval, "SUNLinSol_SPGMRSetMaxRestarts", 1)) return(1);
177
178 /* Attach the linear solver */
179 retval = IDASetLinearSolver(mem, LS, NULL);
180 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
181
182 /* Set the preconditioner solve and setup functions */
183 retval = IDASetPreconditioner(mem, PsetupHeat, PsolveHeat);
184 if(check_retval(&retval, "IDASetPreconditioner", 1)) return(1);
185
186 /* Print output heading. */
187 PrintHeader(rtol, atol);
188
189 /*
190 * -------------------------------------------------------------------------
191 * CASE I
192 * -------------------------------------------------------------------------
193 */
194
195 /* Print case number, output table heading, and initial line of table. */
196
197 printf("\n\nCase 1: gsytpe = MODIFIED_GS\n");
198 printf("\n Output Summary (umax = max-norm of solution) \n\n");
199 printf(" time umax k nst nni nje nre nreLS h npe nps\n" );
200 printf("----------------------------------------------------------------------\n");
201
202 /* Loop over output times, call IDASolve, and print results. */
203
204 for (tout = t1,iout = 1; iout <= NOUT ; iout++, tout *= TWO) {
205 retval = IDASolve(mem, tout, &tret, uu, up, IDA_NORMAL);
206 if(check_retval(&retval, "IDASolve", 1)) return(1);
207 PrintOutput(mem, tret, uu);
208 }
209
210 /* Print remaining counters. */
211
212 retval = IDAGetNumErrTestFails(mem, &netf);
213 check_retval(&retval, "IDAGetNumErrTestFails", 1);
214
215 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
216 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
217
218 retval = IDAGetNumLinConvFails(mem, &ncfl);
219 check_retval(&retval, "IDAGetNumLinConvFails", 1);
220
221 printf("\nError test failures = %ld\n", netf);
222 printf("Nonlinear convergence failures = %ld\n", ncfn);
223 printf("Linear convergence failures = %ld\n", ncfl);
224
225 /*
226 * -------------------------------------------------------------------------
227 * CASE II
228 * -------------------------------------------------------------------------
229 */
230
231 /* Re-initialize uu, up. */
232
233 SetInitialProfile(data, uu, up, res);
234
235 /* Re-initialize IDA and SPGMR */
236
237 retval = IDAReInit(mem, t0, uu, up);
238 if(check_retval(&retval, "IDAReInit", 1)) return(1);
239
240 retval = SUNLinSol_SPGMRSetGSType(LS, CLASSICAL_GS);
241 if(check_retval(&retval, "SUNLinSol_SPGMRSetGSType",1)) return(1);
242
243 /* Print case number, output table heading, and initial line of table. */
244
245 printf("\n\nCase 2: gstype = CLASSICAL_GS\n");
246 printf("\n Output Summary (umax = max-norm of solution) \n\n");
247 printf(" time umax k nst nni nje nre nreLS h npe nps\n" );
248 printf("----------------------------------------------------------------------\n");
249
250 /* Loop over output times, call IDASolve, and print results. */
251
252 for (tout = t1,iout = 1; iout <= NOUT ; iout++, tout *= TWO) {
253 retval = IDASolve(mem, tout, &tret, uu, up, IDA_NORMAL);
254 if(check_retval(&retval, "IDASolve", 1)) return(1);
255 PrintOutput(mem, tret, uu);
256 }
257
258 /* Print remaining counters. */
259
260 retval = IDAGetNumErrTestFails(mem, &netf);
261 check_retval(&retval, "IDAGetNumErrTestFails", 1);
262
263 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
264 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
265
266 retval = IDAGetNumLinConvFails(mem, &ncfl);
267 check_retval(&retval, "IDAGetNumLinConvFails", 1);
268
269 printf("\nError test failures = %ld\n", netf);
270 printf("Nonlinear convergence failures = %ld\n", ncfn);
271 printf("Linear convergence failures = %ld\n", ncfl);
272
273 /* Free Memory */
274
275 IDAFree(&mem);
276 SUNLinSolFree(LS);
277
278 N_VDestroy(uu);
279 N_VDestroy(up);
280 N_VDestroy(res);
281
282 N_VDestroy(data->pp);
283 free(data);
284
285 return(0);
286 }
287
288 /*
289 *--------------------------------------------------------------------
290 * FUNCTIONS CALLED BY IDA
291 *--------------------------------------------------------------------
292 */
293
294 /*
295 * resHeat: heat equation system residual function (user-supplied)
296 * This uses 5-point central differencing on the interior points, and
297 * includes algebraic equations for the boundary values.
298 * So for each interior point, the residual component has the form
299 * res_i = u'_i - (central difference)_i
300 * while for each boundary point, it is res_i = u_i.
301 */
302
resHeat(realtype tt,N_Vector uu,N_Vector up,N_Vector rr,void * user_data)303 int resHeat(realtype tt,
304 N_Vector uu, N_Vector up, N_Vector rr,
305 void *user_data)
306 {
307 sunindextype i, j, offset, loc, mm;
308 realtype *uu_data, *up_data, *rr_data, coeff, dif1, dif2;
309 UserData data;
310
311 uu_data = N_VGetArrayPointer(uu);
312 up_data = N_VGetArrayPointer(up);
313 rr_data = N_VGetArrayPointer(rr);
314
315 data = (UserData) user_data;
316
317 coeff = data->coeff;
318 mm = data->mm;
319
320 /* Initialize rr to uu, to take care of boundary equations. */
321 N_VScale(ONE, uu, rr);
322
323 /* Loop over interior points; set res = up - (central difference). */
324 for (j = 1; j < MGRID-1; j++) {
325 offset = mm*j;
326 for (i = 1; i < mm-1; i++) {
327 loc = offset + i;
328 dif1 = uu_data[loc-1] + uu_data[loc+1] - TWO * uu_data[loc];
329 dif2 = uu_data[loc-mm] + uu_data[loc+mm] - TWO * uu_data[loc];
330 rr_data[loc]= up_data[loc] - coeff * ( dif1 + dif2 );
331 }
332 }
333
334 return(0);
335 }
336
337 /*
338 * PsetupHeat: setup for diagonal preconditioner for idaHeat2D_kry.
339 *
340 * The optional user-supplied functions PsetupHeat and
341 * PsolveHeat together must define the left preconditoner
342 * matrix P approximating the system Jacobian matrix
343 * J = dF/du + cj*dF/du'
344 * (where the DAE system is F(t,u,u') = 0), and solve the linear
345 * systems P z = r. This is done in this case by keeping only
346 * the diagonal elements of the J matrix above, storing them as
347 * inverses in a vector pp, when computed in PsetupHeat, for
348 * subsequent use in PsolveHeat.
349 *
350 * In this instance, only cj and data (user data structure, with
351 * pp etc.) are used from the PsetupdHeat argument list.
352 */
353
PsetupHeat(realtype tt,N_Vector uu,N_Vector up,N_Vector rr,realtype c_j,void * prec_data)354 int PsetupHeat(realtype tt,
355 N_Vector uu, N_Vector up, N_Vector rr,
356 realtype c_j, void *prec_data)
357 {
358
359 sunindextype i, j, offset, loc, mm;
360 realtype *ppv, pelinv;
361 UserData data;
362
363 data = (UserData) prec_data;
364 ppv = N_VGetArrayPointer(data->pp);
365 mm = data->mm;
366
367 /* Initialize the entire vector to 1., then set the interior points to the
368 correct value for preconditioning. */
369 N_VConst(ONE,data->pp);
370
371 /* Compute the inverse of the preconditioner diagonal elements. */
372 pelinv = ONE/(c_j + FOUR*data->coeff);
373
374 for (j = 1; j < mm-1; j++) {
375 offset = mm * j;
376 for (i = 1; i < mm-1; i++) {
377 loc = offset + i;
378 ppv[loc] = pelinv;
379 }
380 }
381
382 return(0);
383 }
384
385 /*
386 * PsolveHeat: solve preconditioner linear system.
387 * This routine multiplies the input vector rvec by the vector pp
388 * containing the inverse diagonal Jacobian elements (previously
389 * computed in PrecondHeateq), returning the result in zvec.
390 */
391
PsolveHeat(realtype tt,N_Vector uu,N_Vector up,N_Vector rr,N_Vector rvec,N_Vector zvec,realtype c_j,realtype delta,void * prec_data)392 int PsolveHeat(realtype tt,
393 N_Vector uu, N_Vector up, N_Vector rr,
394 N_Vector rvec, N_Vector zvec,
395 realtype c_j, realtype delta, void *prec_data)
396 {
397 UserData data;
398 data = (UserData) prec_data;
399 N_VProd(data->pp, rvec, zvec);
400 return(0);
401 }
402
403 /*
404 *--------------------------------------------------------------------
405 * PRIVATE FUNCTIONS
406 *--------------------------------------------------------------------
407 */
408
409 /*
410 * SetInitialProfile: routine to initialize u and up vectors.
411 */
412
SetInitialProfile(UserData data,N_Vector uu,N_Vector up,N_Vector res)413 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
414 N_Vector res)
415 {
416 sunindextype mm, mm1, i, j, offset, loc;
417 realtype xfact, yfact, *udata, *updata;
418
419 mm = data->mm;
420
421 udata = N_VGetArrayPointer(uu);
422 updata = N_VGetArrayPointer(up);
423
424 /* Initialize uu on all grid points. */
425 mm1 = mm - 1;
426 for (j = 0; j < mm; j++) {
427 yfact = data->dx * j;
428 offset = mm*j;
429 for (i = 0;i < mm; i++) {
430 xfact = data->dx * i;
431 loc = offset + i;
432 udata[loc] = RCONST(16.0) * xfact * (ONE - xfact) * yfact * (ONE - yfact);
433 }
434 }
435
436 /* Initialize up vector to 0. */
437 N_VConst(ZERO, up);
438
439 /* resHeat sets res to negative of ODE RHS values at interior points. */
440 resHeat(ZERO, uu, up, res, data);
441
442 /* Copy -res into up to get correct interior initial up values. */
443 N_VScale(-ONE, res, up);
444
445 /* Set up at boundary points to zero. */
446 for (j = 0; j < mm; j++) {
447 offset = mm*j;
448 for (i = 0; i < mm; i++) {
449 loc = offset + i;
450 if (j == 0 || j == mm1 || i == 0 || i == mm1 ) updata[loc] = ZERO;
451 }
452 }
453
454 return(0);
455 }
456
457 /*
458 * Print first lines of output (problem description)
459 */
460
PrintHeader(realtype rtol,realtype atol)461 static void PrintHeader(realtype rtol, realtype atol)
462 {
463 printf("\nidaHeat2D_kry: Heat equation, serial example problem for IDA \n");
464 printf(" Discretized heat equation on 2D unit square. \n");
465 printf(" Zero boundary conditions,");
466 printf(" polynomial initial conditions.\n");
467 printf(" Mesh dimensions: %d x %d", MGRID, MGRID);
468 printf(" Total system size: %d\n\n", NEQ);
469 #if defined(SUNDIALS_EXTENDED_PRECISION)
470 printf("Tolerance parameters: rtol = %Lg atol = %Lg\n", rtol, atol);
471 #elif defined(SUNDIALS_DOUBLE_PRECISION)
472 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
473 #else
474 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
475 #endif
476 printf("Constraints set to force all solution components >= 0. \n");
477 printf("Linear solver: SPGMR, preconditioner using diagonal elements. \n");
478 }
479
480 /*
481 * PrintOutput: print max norm of solution and current solver statistics
482 */
483
PrintOutput(void * mem,realtype t,N_Vector uu)484 static void PrintOutput(void *mem, realtype t, N_Vector uu)
485 {
486 realtype hused, umax;
487 long int nst, nni, nje, nre, nreLS, nli, npe, nps;
488 int kused, retval;
489
490 umax = N_VMaxNorm(uu);
491
492 retval = IDAGetLastOrder(mem, &kused);
493 check_retval(&retval, "IDAGetLastOrder", 1);
494 retval = IDAGetNumSteps(mem, &nst);
495 check_retval(&retval, "IDAGetNumSteps", 1);
496 retval = IDAGetNumNonlinSolvIters(mem, &nni);
497 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
498 retval = IDAGetNumResEvals(mem, &nre);
499 check_retval(&retval, "IDAGetNumResEvals", 1);
500 retval = IDAGetLastStep(mem, &hused);
501 check_retval(&retval, "IDAGetLastStep", 1);
502 retval = IDAGetNumJtimesEvals(mem, &nje);
503 check_retval(&retval, "IDAGetNumJtimesEvals", 1);
504 retval = IDAGetNumLinIters(mem, &nli);
505 check_retval(&retval, "IDAGetNumLinIters", 1);
506 retval = IDAGetNumLinResEvals(mem, &nreLS);
507 check_retval(&retval, "IDAGetNumLinResEvals", 1);
508 retval = IDAGetNumPrecEvals(mem, &npe);
509 check_retval(&retval, "IDAGetNumPrecEvals", 1);
510 retval = IDAGetNumPrecSolves(mem, &nps);
511 check_retval(&retval, "IDAGetNumPrecSolves", 1);
512
513 #if defined(SUNDIALS_EXTENDED_PRECISION)
514 printf(" %5.2Lf %13.5Le %d %3ld %3ld %3ld %4ld %4ld %9.2Le %3ld %3ld\n",
515 t, umax, kused, nst, nni, nje, nre, nreLS, hused, npe, nps);
516 #elif defined(SUNDIALS_DOUBLE_PRECISION)
517 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %4ld %9.2e %3ld %3ld\n",
518 t, umax, kused, nst, nni, nje, nre, nreLS, hused, npe, nps);
519 #else
520 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %4ld %9.2e %3ld %3ld\n",
521 t, umax, kused, nst, nni, nje, nre, nreLS, hused, npe, nps);
522 #endif
523 }
524
525 /*
526 * Check function return value...
527 * opt == 0 means SUNDIALS function allocates memory so check if
528 * returned NULL pointer
529 * opt == 1 means SUNDIALS function returns an integer value so check if
530 * retval < 0
531 * opt == 2 means function allocates memory so check if returned
532 * NULL pointer
533 */
534
check_retval(void * returnvalue,const char * funcname,int opt)535 static int check_retval(void *returnvalue, const char *funcname, int opt)
536 {
537 int *retval;
538
539 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
540 if (opt == 0 && returnvalue == NULL) {
541 fprintf(stderr,
542 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
543 funcname);
544 return(1);
545 } else if (opt == 1) {
546 /* Check if retval < 0 */
547 retval = (int *) returnvalue;
548 if (*retval < 0) {
549 fprintf(stderr,
550 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
551 funcname, *retval);
552 return(1);
553 }
554 } else if (opt == 2 && returnvalue == NULL) {
555 /* Check if function returned NULL pointer - no memory allocated */
556 fprintf(stderr,
557 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
558 funcname);
559 return(1);
560 }
561
562 return(0);
563 }
564