1 /* -----------------------------------------------------------------
2 * Programmer(s): Allan Taylor, Alan Hindmarsh and
3 * 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 for IDA: 2D heat equation, serial, banded.
16 *
17 * This example solves a discretized 2D heat equation problem.
18 * This version uses the band solver and IDACalcIC.
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).
24 * The PDE is treated with central differences on a uniform M x M
25 * grid. 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 banded linear system
30 * solver, half-bandwidths equal to M, and default
31 * difference-quotient Jacobian. For purposes of illustration,
32 * IDACalcIC is called to compute correct values at the boundary,
33 * given incorrect values as input initial guesses. The constraints
34 * u >= 0 are posed for all components. Output is taken at
35 * t = 0, .01, .02, .04, ..., 10.24. (Output at t = 0 is for
36 * IDACalcIC cost statistics only.)
37 * -----------------------------------------------------------------*/
38
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <math.h>
42
43 #include <idas/idas.h> /* prototypes for IDA fcts., consts. */
44 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
45 #include <sunmatrix/sunmatrix_band.h> /* access to band SUNMatrix */
46 #include <sunlinsol/sunlinsol_band.h> /* access to band SUNLinearSolver */
47 #include <sundials/sundials_types.h> /* definition of type realtype */
48
49 /* Problem Constants */
50
51 #define NOUT 11
52 #define MGRID 10
53 #define NEQ MGRID*MGRID
54 #define ZERO RCONST(0.0)
55 #define ONE RCONST(1.0)
56 #define TWO RCONST(2.0)
57 #define BVAL RCONST(0.1)
58
59 /* Type: UserData */
60
61 typedef struct {
62 sunindextype mm;
63 realtype dx;
64 realtype coeff;
65 } *UserData;
66
67 /* Prototypes of functions called by IDA */
68
69 int heatres(realtype tres, N_Vector uu, N_Vector up, N_Vector resval, void *user_data);
70
71 /* Prototypes of private functions */
72
73 static void PrintHeader(realtype rtol, realtype atol);
74 static void PrintOutput(void *mem, realtype t, N_Vector u);
75 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
76 N_Vector id, N_Vector res);
77
78 static int check_retval(void *returnvalue, const char *funcname, int opt);
79
80 /*
81 *--------------------------------------------------------------------
82 * MAIN PROGRAM
83 *--------------------------------------------------------------------
84 */
85
main(void)86 int main(void)
87 {
88 void *mem;
89 UserData data;
90 N_Vector uu, up, constraints, id, res;
91 int retval, iout;
92 long int netf, ncfn;
93 sunindextype mu, ml;
94 realtype rtol, atol, t0, t1, tout, tret;
95 SUNMatrix A;
96 SUNLinearSolver LS;
97
98 mem = NULL;
99 data = NULL;
100 uu = up = constraints = id = res = NULL;
101 A = NULL;
102 LS = NULL;
103
104 /* Create vectors uu, up, res, constraints, id. */
105 uu = N_VNew_Serial(NEQ);
106 if(check_retval((void *)uu, "N_VNew_Serial", 0)) return(1);
107 up = N_VNew_Serial(NEQ);
108 if(check_retval((void *)up, "N_VNew_Serial", 0)) return(1);
109 res = N_VNew_Serial(NEQ);
110 if(check_retval((void *)res, "N_VNew_Serial", 0)) return(1);
111 constraints = N_VNew_Serial(NEQ);
112 if(check_retval((void *)constraints, "N_VNew_Serial", 0)) return(1);
113 id = N_VNew_Serial(NEQ);
114 if(check_retval((void *)id, "N_VNew_Serial", 0)) return(1);
115
116 /* Create and load problem data block. */
117 data = (UserData) malloc(sizeof *data);
118 if(check_retval((void *)data, "malloc", 2)) return(1);
119 data->mm = MGRID;
120 data->dx = ONE/(MGRID - ONE);
121 data->coeff = ONE/( (data->dx) * (data->dx) );
122
123 /* Initialize uu, up, id. */
124 SetInitialProfile(data, uu, up, id, res);
125
126 /* Set constraints to all 1's for nonnegative solution values. */
127 N_VConst(ONE, constraints);
128
129 /* Set remaining input parameters. */
130 t0 = ZERO;
131 t1 = RCONST(0.01);
132 rtol = ZERO;
133 atol = RCONST(1.0e-3);
134
135 /* Call IDACreate and IDAMalloc to initialize solution */
136 mem = IDACreate();
137 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
138
139 retval = IDASetUserData(mem, data);
140 if(check_retval(&retval, "IDASetUserData", 1)) return(1);
141
142 /* Set which components are algebraic or differential */
143 retval = IDASetId(mem, id);
144 if(check_retval(&retval, "IDASetId", 1)) return(1);
145
146 retval = IDASetConstraints(mem, constraints);
147 if(check_retval(&retval, "IDASetConstraints", 1)) return(1);
148 N_VDestroy(constraints);
149
150 retval = IDAInit(mem, heatres, t0, uu, up);
151 if(check_retval(&retval, "IDAInit", 1)) return(1);
152
153 retval = IDASStolerances(mem, rtol, atol);
154 if(check_retval(&retval, "IDASStolerances", 1)) return(1);
155
156 /* Create banded SUNMatrix for use in linear solves */
157 mu = MGRID; ml = MGRID;
158 A = SUNBandMatrix(NEQ, mu, ml);
159 if(check_retval((void *)A, "SUNBandMatrix", 0)) return(1);
160
161 /* Create banded SUNLinearSolver object */
162 LS = SUNLinSol_Band(uu, A);
163 if(check_retval((void *)LS, "SUNLinSol_Band", 0)) return(1);
164
165 /* Attach the matrix and linear solver */
166 retval = IDASetLinearSolver(mem, LS, A);
167 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
168
169 /* Call IDACalcIC to correct the initial values. */
170
171 retval = IDACalcIC(mem, IDA_YA_YDP_INIT, t1);
172 if(check_retval(&retval, "IDACalcIC", 1)) return(1);
173
174 /* Print output heading. */
175 PrintHeader(rtol, atol);
176
177 PrintOutput(mem, t0, uu);
178
179
180 /* Loop over output times, call IDASolve, and print results. */
181
182 for (tout = t1, iout = 1; iout <= NOUT; iout++, tout *= TWO) {
183
184 retval = IDASolve(mem, tout, &tret, uu, up, IDA_NORMAL);
185 if(check_retval(&retval, "IDASolve", 1)) return(1);
186
187 PrintOutput(mem, tret, uu);
188
189 }
190
191 /* Print remaining counters and free memory. */
192 retval = IDAGetNumErrTestFails(mem, &netf);
193 check_retval(&retval, "IDAGetNumErrTestFails", 1);
194 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
195 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
196 printf("\n netf = %ld, ncfn = %ld \n", netf, ncfn);
197
198 IDAFree(&mem);
199 SUNLinSolFree(LS);
200 SUNMatDestroy(A);
201 N_VDestroy(uu);
202 N_VDestroy(up);
203 N_VDestroy(id);
204 N_VDestroy(res);
205 free(data);
206
207 return(0);
208 }
209
210 /*
211 *--------------------------------------------------------------------
212 * FUNCTIONS CALLED BY IDA
213 *--------------------------------------------------------------------
214 */
215
216 /*
217 * heatres: heat equation system residual function
218 * This uses 5-point central differencing on the interior points, and
219 * includes algebraic equations for the boundary values.
220 * So for each interior point, the residual component has the form
221 * res_i = u'_i - (central difference)_i
222 * while for each boundary point, it is res_i = u_i.
223 */
224
heatres(realtype tres,N_Vector uu,N_Vector up,N_Vector resval,void * user_data)225 int heatres(realtype tres, N_Vector uu, N_Vector up, N_Vector resval,
226 void *user_data)
227 {
228 sunindextype mm, i, j, offset, loc;
229 realtype *uv, *upv, *resv, coeff;
230 UserData data;
231
232 uv = N_VGetArrayPointer(uu); upv = N_VGetArrayPointer(up); resv = N_VGetArrayPointer(resval);
233
234 data = (UserData)user_data;
235 mm = data->mm;
236 coeff = data->coeff;
237
238 /* Initialize resval to uu, to take care of boundary equations. */
239 N_VScale(ONE, uu, resval);
240
241 /* Loop over interior points; set res = up - (central difference). */
242 for (j = 1; j < mm-1; j++) {
243 offset = mm*j;
244 for (i = 1; i < mm-1; i++) {
245 loc = offset + i;
246 resv[loc] = upv[loc] - coeff *
247 (uv[loc-1] + uv[loc+1] + uv[loc-mm] + uv[loc+mm] - RCONST(4.0)*uv[loc]);
248 }
249 }
250
251 return(0);
252
253 }
254
255 /*
256 *--------------------------------------------------------------------
257 * PRIVATE FUNCTIONS
258 *--------------------------------------------------------------------
259 */
260
261 /*
262 * SetInitialProfile: routine to initialize u, up, and id vectors.
263 */
264
SetInitialProfile(UserData data,N_Vector uu,N_Vector up,N_Vector id,N_Vector res)265 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
266 N_Vector id, N_Vector res)
267 {
268 realtype xfact, yfact, *udata, *updata, *iddata;
269 sunindextype mm, mm1, i, j, offset, loc;
270
271 mm = data->mm;
272 mm1 = mm - 1;
273
274 udata = N_VGetArrayPointer(uu);
275 updata = N_VGetArrayPointer(up);
276 iddata = N_VGetArrayPointer(id);
277
278 /* Initialize id to 1's. */
279 N_VConst(ONE, id);
280
281 /* Initialize uu on all grid points. */
282 for (j = 0; j < mm; j++) {
283 yfact = data->dx * j;
284 offset = mm*j;
285 for (i = 0;i < mm; i++) {
286 xfact = data->dx * i;
287 loc = offset + i;
288 udata[loc] = RCONST(16.0) * xfact * (ONE - xfact) * yfact * (ONE - yfact);
289 }
290 }
291
292 /* Initialize up vector to 0. */
293 N_VConst(ZERO, up);
294
295 /* heatres sets res to negative of ODE RHS values at interior points. */
296 heatres(ZERO, uu, up, res, data);
297
298 /* Copy -res into up to get correct interior initial up values. */
299 N_VScale(-ONE, res, up);
300
301 /* Finally, set values of u, up, and id at boundary points. */
302 for (j = 0; j < mm; j++) {
303 offset = mm*j;
304 for (i = 0;i < mm; i++) {
305 loc = offset + i;
306 if (j == 0 || j == mm1 || i == 0 || i == mm1 ) {
307 udata[loc] = BVAL; updata[loc] = ZERO; iddata[loc] = ZERO; }
308 }
309 }
310
311 return(0);
312
313 }
314
315 /*
316 * Print first lines of output (problem description)
317 */
318
PrintHeader(realtype rtol,realtype atol)319 static void PrintHeader(realtype rtol, realtype atol)
320 {
321 printf("\nidasHeat2D_bnd: Heat equation, serial example problem for IDA\n");
322 printf(" Discretized heat equation on 2D unit square.\n");
323 printf(" Zero boundary conditions,");
324 printf(" polynomial initial conditions.\n");
325 printf(" Mesh dimensions: %d x %d", MGRID, MGRID);
326 printf(" Total system size: %d\n\n", NEQ);
327 #if defined(SUNDIALS_EXTENDED_PRECISION)
328 printf("Tolerance parameters: rtol = %Lg atol = %Lg\n", rtol, atol);
329 #elif defined(SUNDIALS_DOUBLE_PRECISION)
330 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
331 #else
332 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
333 #endif
334 printf("Constraints set to force all solution components >= 0. \n");
335 printf("Linear solver: BAND, banded direct solver \n");
336 printf(" difference quotient Jacobian, half-bandwidths = %d \n",MGRID);
337 #if defined(SUNDIALS_EXTENDED_PRECISION)
338 printf("IDACalcIC called with input boundary values = %Lg \n",BVAL);
339 #elif defined(SUNDIALS_DOUBLE_PRECISION)
340 printf("IDACalcIC called with input boundary values = %g \n",BVAL);
341 #else
342 printf("IDACalcIC called with input boundary values = %g \n",BVAL);
343 #endif
344 /* Print output table heading and initial line of table. */
345 printf("\n Output Summary (umax = max-norm of solution) \n\n");
346 printf(" time umax k nst nni nje nre nreLS h \n" );
347 printf(" . . . . . . . . . . . . . . . . . . . . . \n");
348 }
349
350 /*
351 * Print Output
352 */
353
PrintOutput(void * mem,realtype t,N_Vector uu)354 static void PrintOutput(void *mem, realtype t, N_Vector uu)
355 {
356 int retval;
357 realtype umax, hused;
358 long int nst, nni, nje, nre, nreLS;
359 int kused;
360
361 umax = N_VMaxNorm(uu);
362
363 retval = IDAGetLastOrder(mem, &kused);
364 check_retval(&retval, "IDAGetLastOrder", 1);
365 retval = IDAGetNumSteps(mem, &nst);
366 check_retval(&retval, "IDAGetNumSteps", 1);
367 retval = IDAGetNumNonlinSolvIters(mem, &nni);
368 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
369 retval = IDAGetNumResEvals(mem, &nre);
370 check_retval(&retval, "IDAGetNumResEvals", 1);
371 retval = IDAGetLastStep(mem, &hused);
372 check_retval(&retval, "IDAGetLastStep", 1);
373 retval = IDAGetNumJacEvals(mem, &nje);
374 check_retval(&retval, "IDAGetNumJacEvals", 1);
375 retval = IDAGetNumLinResEvals(mem, &nreLS);
376 check_retval(&retval, "IDAGetNumLinResEvals", 1);
377
378 #if defined(SUNDIALS_EXTENDED_PRECISION)
379 printf(" %5.2Lf %13.5Le %d %3ld %3ld %3ld %4ld %4ld %9.2Le \n",
380 t, umax, kused, nst, nni, nje, nre, nreLS, hused);
381 #elif defined(SUNDIALS_DOUBLE_PRECISION)
382 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %4ld %9.2e \n",
383 t, umax, kused, nst, nni, nje, nre, nreLS, hused);
384 #else
385 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %4ld %9.2e \n",
386 t, umax, kused, nst, nni, nje, nre, nreLS, hused);
387 #endif
388
389 }
390
391 /*
392 * Check function return value...
393 * opt == 0 means SUNDIALS function allocates memory so check if
394 * returned NULL pointer
395 * opt == 1 means SUNDIALS function returns an integer value so check if
396 * retval < 0
397 * opt == 2 means function allocates memory so check if returned
398 * NULL pointer
399 */
400
check_retval(void * returnvalue,const char * funcname,int opt)401 static int check_retval(void *returnvalue, const char *funcname, int opt)
402 {
403 int *retval;
404
405 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
406 if (opt == 0 && returnvalue == NULL) {
407 fprintf(stderr,
408 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
409 funcname);
410 return(1);
411 } else if (opt == 1) {
412 /* Check if retval < 0 */
413 retval = (int *) returnvalue;
414 if (*retval < 0) {
415 fprintf(stderr,
416 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
417 funcname, *retval);
418 return(1);
419 }
420 } else if (opt == 2 && returnvalue == NULL) {
421 /* Check if function returned NULL pointer - no memory allocated */
422 fprintf(stderr,
423 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
424 funcname);
425 return(1);
426 }
427
428 return(0);
429 }
430