1 /* -----------------------------------------------------------------
2 * Programmer(s): Chris Nguyen @ LLNL
3 * based on idaHeat2D_bnd.c and idaRoberts_klu.c
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, sparse.
16 *
17 * This example solves a discretized 2D heat equation problem.
18 * This version uses the KLU 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 MGRID x MGRID
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 = MGRID^2. Here MGRID = 10.
28 *
29 * The system is solved with IDA using the sparse linear system
30 * solver and a user supplied Jacobian.
31 * 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 <ida/ida.h> /* prototypes for IDA fcts., consts. */
44 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
45 #include <sunmatrix/sunmatrix_sparse.h> /* access to sparse SUNMatrix */
46 #include <sunlinsol/sunlinsol_klu.h> /* access to KLU linear solver */
47 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
48 #include <sundials/sundials_math.h> /* defs. of SUNRabs, SUNRexp, etc. */
49
50 /* Problem Constants */
51
52 #define NOUT 11
53 #define MGRID 10
54 #define NEQ MGRID*MGRID
55 #define ZERO RCONST(0.0)
56 #define ONE RCONST(1.0)
57 #define TWO RCONST(2.0)
58 #define BVAL RCONST(0.0)
59 #define TOTAL 4*MGRID+8*(MGRID-2)+(MGRID-4)*(MGRID+4*(MGRID-2)) /* total num of nonzero elements */
60
61 /* Type: UserData */
62
63 typedef struct {
64 sunindextype mm;
65 realtype dx;
66 realtype coeff;
67 } *UserData;
68
69 /* Prototypes of functions called by IDA */
70
71 int heatres(realtype tres, N_Vector uu, N_Vector up, N_Vector resval, void *user_data);
72
73 int jacHeat(realtype tt, realtype cj,
74 N_Vector yy, N_Vector yp, N_Vector resvec,
75 SUNMatrix JJ, void *user_data,
76 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
77
78 /* Exact same setup as jacHeat. Function needed for special case MGRID=3 */
79 int jacHeat3(realtype tt, realtype cj,
80 N_Vector yy, N_Vector yp, N_Vector resvec,
81 SUNMatrix JJ, void *user_data,
82 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
83
84 /* Prototypes of private functions */
85
86 static void PrintHeader(realtype rtol, realtype atol);
87 static void PrintOutput(void *mem, realtype t, N_Vector u);
88 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
89 N_Vector id, N_Vector res);
90
91 static int check_retval(void *returnvalue, const char *funcname, int opt);
92
93 /*
94 *--------------------------------------------------------------------
95 * MAIN PROGRAM
96 *--------------------------------------------------------------------
97 */
98
main(void)99 int main(void)
100 {
101 void *mem;
102 UserData data;
103 N_Vector uu, up, constraints, id, res;
104 int retval, iout;
105 long int netf, ncfn;
106 realtype rtol, atol, t0, t1, tout, tret;
107 SUNMatrix A;
108 SUNLinearSolver LS;
109 sunindextype nnz;
110
111 mem = NULL;
112 data = NULL;
113 uu = up = constraints = id = res = NULL;
114 A = NULL;
115 LS = NULL;
116
117 /* Create vectors uu, up, res, constraints, id. */
118 uu = N_VNew_Serial(NEQ);
119 if(check_retval((void *)uu, "N_VNew_Serial", 0)) return(1);
120 up = N_VNew_Serial(NEQ);
121 if(check_retval((void *)up, "N_VNew_Serial", 0)) return(1);
122 res = N_VNew_Serial(NEQ);
123 if(check_retval((void *)res, "N_VNew_Serial", 0)) return(1);
124 constraints = N_VNew_Serial(NEQ);
125 if(check_retval((void *)constraints, "N_VNew_Serial", 0)) return(1);
126 id = N_VNew_Serial(NEQ);
127 if(check_retval((void *)id, "N_VNew_Serial", 0)) return(1);
128
129 /* Create and load problem data block. */
130 data = (UserData) malloc(sizeof *data);
131 if(check_retval((void *)data, "malloc", 2)) return(1);
132 data->mm = MGRID;
133 data->dx = ONE/(MGRID - ONE);
134 data->coeff = ONE/( (data->dx) * (data->dx) );
135
136 /* Initialize uu, up, id. */
137 SetInitialProfile(data, uu, up, id, res);
138
139 /* Set constraints to all 1's for nonnegative solution values. */
140 N_VConst(ONE, constraints);
141
142 /* Set remaining input parameters. */
143 t0 = ZERO;
144 t1 = RCONST(0.01);
145 rtol = ZERO;
146 atol = RCONST(1.0e-8);
147
148 /* Call IDACreate and IDAMalloc to initialize solution */
149 mem = IDACreate();
150 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
151
152 retval = IDASetUserData(mem, data);
153 if(check_retval(&retval, "IDASetUserData", 1)) return(1);
154
155 /* Set which components are algebraic or differential */
156 retval = IDASetId(mem, id);
157 if(check_retval(&retval, "IDASetId", 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, heatres, 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 sparse SUNMatrix for use in linear solves */
170 nnz = NEQ*NEQ;
171 A = SUNSparseMatrix(NEQ, NEQ, nnz, CSC_MAT);
172 if(check_retval((void*)A, "SUNSparseMtarix", 0)) return(1);
173
174 /* Create KLU SUNLinearSolver object */
175 LS = SUNLinSol_KLU(uu, A);
176 if(check_retval((void *)LS, "SUNLinSol_KLU", 0)) return(1);
177
178 /* Attach the matrix and linear solver */
179 retval = IDASetLinearSolver(mem, LS, A);
180 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
181
182 /* Set the user-supplied Jacobian routine */
183 if(MGRID >= 4){
184 retval = IDASetJacFn(mem, jacHeat);
185 } else if(MGRID == 3) {
186 retval = IDASetJacFn(mem, jacHeat3);
187 } else {
188 /* MGRID<=2 is pure boundary points, nothing to solve */
189 printf("MGRID size is too small to run.\n");
190 return(1);
191 }
192 if(check_retval(&retval, "IDASetJacFn", 1)) return(1);
193
194 /* Call IDACalcIC to correct the initial values. */
195
196 retval = IDACalcIC(mem, IDA_YA_YDP_INIT, t1);
197 if(check_retval(&retval, "IDACalcIC", 1)) return(1);
198
199 /* Print output heading. */
200 PrintHeader(rtol, atol);
201
202 PrintOutput(mem, t0, uu);
203
204
205 /* Loop over output times, call IDASolve, and print results. */
206
207 for (tout = t1, iout = 1; iout <= NOUT; iout++, tout *= TWO) {
208
209 retval = IDASolve(mem, tout, &tret, uu, up, IDA_NORMAL);
210 if(check_retval(&retval, "IDASolve", 1)) return(1);
211
212 PrintOutput(mem, tret, uu);
213
214 }
215
216 /* Print remaining counters and free memory. */
217 retval = IDAGetNumErrTestFails(mem, &netf);
218 check_retval(&retval, "IDAGetNumErrTestFails", 1);
219 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
220 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
221 printf("\n netf = %ld, ncfn = %ld \n", netf, ncfn);
222
223 IDAFree(&mem);
224 SUNLinSolFree(LS);
225 SUNMatDestroy(A);
226 N_VDestroy(uu);
227 N_VDestroy(up);
228 N_VDestroy(id);
229 N_VDestroy(res);
230 free(data);
231
232 return(0);
233 }
234
235 /*
236 *--------------------------------------------------------------------
237 * FUNCTIONS CALLED BY IDA
238 *--------------------------------------------------------------------
239 */
240
241 /*
242 * heatres: heat equation system residual function
243 * This uses 5-point central differencing on the interior points, and
244 * includes algebraic equations for the boundary values.
245 * So for each interior point, the residual component has the form
246 * res_i = u'_i - (central difference)_i
247 * while for each boundary point, it is res_i = u_i.
248 */
249
heatres(realtype tres,N_Vector uu,N_Vector up,N_Vector resval,void * user_data)250 int heatres(realtype tres, N_Vector uu, N_Vector up, N_Vector resval,
251 void *user_data)
252 {
253 sunindextype mm, i, j, offset, loc;
254 realtype *uv, *upv, *resv, coeff;
255 UserData data;
256
257 uv = N_VGetArrayPointer(uu); upv = N_VGetArrayPointer(up); resv = N_VGetArrayPointer(resval);
258
259 data = (UserData)user_data;
260 mm = data->mm;
261 coeff = data->coeff;
262
263 /* Initialize resval to uu, to take care of boundary equations. */
264 N_VScale(ZERO, uu, resval);
265
266 /* Loop over interior points; set res = up - (central difference). */
267 for (j = 1; j < mm-1; j++) {
268 offset = mm*j;
269 for (i = 1; i < mm-1; i++) {
270 loc = offset + i;
271 resv[loc] = upv[loc] - coeff *
272 (uv[loc-1] + uv[loc+1] + uv[loc-mm] + uv[loc+mm] - RCONST(4.0)*uv[loc]);
273 }
274 }
275
276 return(0);
277
278 }
279
280 /* Jacobian matrix setup for MGRID=3 */
jacHeat3(realtype tt,realtype cj,N_Vector yy,N_Vector yp,N_Vector resvec,SUNMatrix JJ,void * user_data,N_Vector tempv1,N_Vector tempv2,N_Vector tempv3)281 int jacHeat3(realtype tt, realtype cj,
282 N_Vector yy, N_Vector yp, N_Vector resvec,
283 SUNMatrix JJ, void *user_data,
284 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3)
285 {
286 realtype dx = ONE/(MGRID - ONE);
287 realtype beta = RCONST(4.0)/(dx*dx) + cj;
288
289 sunindextype *colptrs = SUNSparseMatrix_IndexPointers(JJ);
290 sunindextype *rowvals = SUNSparseMatrix_IndexValues(JJ);
291 realtype *data = SUNSparseMatrix_Data(JJ);
292
293 SUNMatZero(JJ);
294
295 /*
296 * set up number of elements in each column
297 */
298 colptrs[0] = 0;
299 colptrs[1] = 1;
300 colptrs[2] = 3;
301 colptrs[3] = 4;
302 colptrs[4] = 6;
303 colptrs[5] = 7;
304 colptrs[6] = 9;
305 colptrs[7] = 10;
306 colptrs[8] = 12;
307 colptrs[9] = 13;
308
309 /*
310 * set up data and row values stored
311 */
312
313 data[0] = ONE;
314 rowvals[0] = 0;
315 data[1] = ONE;
316 rowvals[1] = 1;
317 data[2] = -ONE/(dx*dx);
318 rowvals[2] = 4;
319 data[3] = ONE;
320 rowvals[3] = 2;
321 data[4] = ONE;
322 rowvals[4] = 3;
323 data[5] = -ONE/(dx*dx);
324 rowvals[5] = 4;
325 data[6] = beta;
326 rowvals[6] = 4;
327 data[7] = -ONE/(dx*dx);
328 rowvals[7] = 4;
329 data[8] = ONE;
330 rowvals[8] = 5;
331 data[9] = ONE;
332 rowvals[9] = 6;
333 data[10] = -ONE/(dx*dx);
334 rowvals[10] = 4;
335 data[11] = ONE;
336 rowvals[11] = 7;
337 data[12] = ONE;
338 rowvals[12] = 8;
339
340 return(0);
341 }
342
343 /* Jacobian matrix setup for MGRID>=4 */
jacHeat(realtype tt,realtype cj,N_Vector yy,N_Vector yp,N_Vector resvec,SUNMatrix JJ,void * user_data,N_Vector tempv1,N_Vector tempv2,N_Vector tempv3)344 int jacHeat(realtype tt, realtype cj,
345 N_Vector yy, N_Vector yp, N_Vector resvec,
346 SUNMatrix JJ, void *user_data,
347 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3)
348 {
349 realtype dx = ONE/(MGRID - ONE);
350 realtype beta = RCONST(4.0)/(dx*dx) + cj;
351 int i,j, repeat=0;
352
353 sunindextype *colptrs = SUNSparseMatrix_IndexPointers(JJ);
354 sunindextype *rowvals = SUNSparseMatrix_IndexValues(JJ);
355 realtype *data = SUNSparseMatrix_Data(JJ);
356
357 SUNMatZero(JJ);
358
359 /*
360 *-----------------------------------------------
361 * set up number of elements in each column
362 *-----------------------------------------------
363 */
364
365 /**** first column block ****/
366 colptrs[0] = 0;
367 colptrs[1] = 1;
368 /* count by twos in the middle */
369 for(i=2;i<MGRID;i++) colptrs[i] = (colptrs[i-1])+2;
370 colptrs[MGRID] = 2*MGRID-2;
371
372 /**** second column block ****/
373 colptrs[MGRID+1] = 2*MGRID;
374 colptrs[MGRID+2] = 2*MGRID+3;
375 /* count by fours in the middle */
376 for(i=0;i<MGRID-4;i++) colptrs[MGRID+3+i] = (colptrs[MGRID+3+i-1])+4;
377 colptrs[2*MGRID-1] = 2*MGRID+4*(MGRID-2)-2;
378 colptrs[2*MGRID] = 2*MGRID+4*(MGRID-2);
379
380 /**** repeated (MGRID-4 times) middle column blocks ****/
381 for(i=0;i<MGRID-4;i++){
382 colptrs[2*MGRID+1+repeat] = (colptrs[2*MGRID+1+repeat-1])+2;
383 colptrs[2*MGRID+1+repeat+1] = (colptrs[2*MGRID+1+repeat])+4;
384
385 /* count by fives in the middle */
386 for(j=0;j<MGRID-4;j++) colptrs[2*MGRID+1+repeat+2+j] =
387 (colptrs[2*MGRID+1+repeat+1+j])+5;
388
389 colptrs[2*MGRID+1+repeat+(MGRID-4)+2] = (colptrs[2*MGRID+1+repeat+(MGRID-4)+1])+4;
390 colptrs[2*MGRID+1+repeat+(MGRID-4)+3] = (colptrs[2*MGRID+1+repeat+(MGRID-4)+2])+2;
391
392 repeat+=MGRID; /* shift that accounts for accumulated number of columns */
393 }
394
395 /**** last-1 column block ****/
396 colptrs[MGRID*MGRID-2*MGRID+1] = TOTAL-2*MGRID-4*(MGRID-2)+2;
397 colptrs[MGRID*MGRID-2*MGRID+2] = TOTAL-2*MGRID-4*(MGRID-2)+5;
398 /* count by fours in the middle */
399 for(i=0;i<MGRID-4;i++) colptrs[MGRID*MGRID-2*MGRID+3+i] =
400 (colptrs[MGRID*MGRID-2*MGRID+3+i-1])+4;
401 colptrs[MGRID*MGRID-MGRID-1] = TOTAL-2*MGRID;
402 colptrs[MGRID*MGRID-MGRID] = TOTAL-2*MGRID+2;
403
404 /**** last column block ****/
405 colptrs[MGRID*MGRID-MGRID+1] = TOTAL-MGRID-(MGRID-2)+1;
406 /* count by twos in the middle */
407 for(i=0;i<MGRID-2;i++) colptrs[MGRID*MGRID-MGRID+2+i] =
408 (colptrs[MGRID*MGRID-MGRID+2+i-1])+2;
409 colptrs[MGRID*MGRID-1] = TOTAL-1;
410 colptrs[MGRID*MGRID] = TOTAL;
411
412
413 /*
414 *-----------------------------------------------
415 * set up data stored
416 *-----------------------------------------------
417 */
418
419 /**** first column block ****/
420 data[0] = ONE;
421 /* alternating pattern in data, separate loop for each pattern */
422 for(i=1;i<MGRID+(MGRID-2) ;i+=2) data[i] = ONE;
423 for(i=2;i<MGRID+(MGRID-2)-1;i+=2) data[i] = -ONE/(dx*dx);
424
425 /**** second column block ****/
426 data[MGRID+MGRID-2] = ONE;
427 data[MGRID+MGRID-1] = -ONE/(dx*dx);
428 data[MGRID+MGRID] = beta;
429 data[MGRID+MGRID+1] = -ONE/(dx*dx);
430 data[MGRID+MGRID+2] = -ONE/(dx*dx);
431 /* middle data elements */
432 for(i=0;i<(MGRID-4);i++) data[MGRID+MGRID+3+4*i] = -ONE/(dx*dx);
433 for(i=0;i<(MGRID-4);i++) data[MGRID+MGRID+4+4*i] = beta;
434 for(i=0;i<(MGRID-4);i++) data[MGRID+MGRID+5+4*i] = -ONE/(dx*dx);
435 for(i=0;i<(MGRID-4);i++) data[MGRID+MGRID+6+4*i] = -ONE/(dx*dx);
436 data[2*MGRID+4*(MGRID-2)-5] = -ONE/(dx*dx);
437 data[2*MGRID+4*(MGRID-2)-4] = beta;
438 data[2*MGRID+4*(MGRID-2)-3] = -ONE/(dx*dx);
439 data[2*MGRID+4*(MGRID-2)-2] = -ONE/(dx*dx);
440 data[2*MGRID+4*(MGRID-2)-1] = ONE;
441
442 /**** repeated (MGRID-4 times) middle column blocks ****/
443 repeat=0;
444 for(i=0;i<MGRID-4;i++){
445 data[2*MGRID+4*(MGRID-2)+repeat] = ONE;
446 data[2*MGRID+4*(MGRID-2)+repeat+1] = -ONE/(dx*dx);
447
448 data[2*MGRID+4*(MGRID-2)+repeat+2] = -ONE/(dx*dx);
449 data[2*MGRID+4*(MGRID-2)+repeat+3] = beta;
450 data[2*MGRID+4*(MGRID-2)+repeat+4] = -ONE/(dx*dx);
451 data[2*MGRID+4*(MGRID-2)+repeat+5] = -ONE/(dx*dx);
452
453 /* 5 in 5*j chosen since there are 5 elements in each column */
454 /* this column loops MGRID-4 times within the outer loop */
455 for(j=0;j<MGRID-4;j++){
456 data[2*MGRID+4*(MGRID-2)+repeat+6+5*j] = -ONE/(dx*dx);
457 data[2*MGRID+4*(MGRID-2)+repeat+7+5*j] = -ONE/(dx*dx);
458 data[2*MGRID+4*(MGRID-2)+repeat+8+5*j] = beta;
459 data[2*MGRID+4*(MGRID-2)+repeat+9+5*j] = -ONE/(dx*dx);
460 data[2*MGRID+4*(MGRID-2)+repeat+10+5*j] = -ONE/(dx*dx);
461 }
462
463 data[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+6] = -ONE/(dx*dx);
464 data[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+7] = -ONE/(dx*dx);
465 data[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+8] = beta;
466 data[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+9] = -ONE/(dx*dx);
467
468 data[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+10] = -ONE/(dx*dx);
469 data[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+11] = ONE;
470
471 repeat+=MGRID+4*(MGRID-2); /* shift that accounts for accumulated columns and elements */
472 }
473
474 /**** last-1 column block ****/
475 data[TOTAL-6*(MGRID-2)-4] = ONE;
476 data[TOTAL-6*(MGRID-2)-3] = -ONE/(dx*dx);
477 data[TOTAL-6*(MGRID-2)-2] = -ONE/(dx*dx);
478 data[TOTAL-6*(MGRID-2)-1] = beta;
479 data[TOTAL-6*(MGRID-2) ] = -ONE/(dx*dx);
480 /* middle data elements */
481 for(i=0;i<(MGRID-4);i++) data[TOTAL-6*(MGRID-2)+1+4*i] = -ONE/(dx*dx);
482 for(i=0;i<(MGRID-4);i++) data[TOTAL-6*(MGRID-2)+2+4*i] = -ONE/(dx*dx);
483 for(i=0;i<(MGRID-4);i++) data[TOTAL-6*(MGRID-2)+3+4*i] = beta;
484 for(i=0;i<(MGRID-4);i++) data[TOTAL-6*(MGRID-2)+4+4*i] = -ONE/(dx*dx);
485 data[TOTAL-2*(MGRID-2)-7] = -ONE/(dx*dx);
486 data[TOTAL-2*(MGRID-2)-6] = -ONE/(dx*dx);
487 data[TOTAL-2*(MGRID-2)-5] = beta;
488 data[TOTAL-2*(MGRID-2)-4] = -ONE/(dx*dx);
489 data[TOTAL-2*(MGRID-2)-3] = ONE;
490
491 /**** last column block ****/
492 data[TOTAL-2*(MGRID-2)-2] = ONE;
493 /* alternating pattern in data, separate loop for each pattern */
494 for(i=TOTAL-2*(MGRID-2)-1;i<TOTAL-2;i+=2) data[i] = -ONE/(dx*dx);
495 for(i=TOTAL-2*(MGRID-2) ;i<TOTAL-1;i+=2) data[i] = ONE;
496 data[TOTAL-1] = ONE;
497
498 /*
499 *-----------------------------------------------
500 * row values
501 *-----------------------------------------------
502 */
503
504 /**** first block ****/
505 rowvals[0] = 0;
506 /* alternating pattern in data, separate loop for each pattern */
507 for(i=1;i<MGRID+(MGRID-2) ;i+=2) rowvals[i] = (i+1)/2;
508 for(i=2;i<MGRID+(MGRID-2)-1;i+=2) rowvals[i] = i/2+MGRID; /* i+1 unnecessary here */
509
510 /**** second column block ****/
511 rowvals[MGRID+MGRID-2] = MGRID;
512 rowvals[MGRID+MGRID-1] = MGRID+1;
513 rowvals[MGRID+MGRID] = MGRID+1;
514 rowvals[MGRID+MGRID+1] = MGRID+2;
515 rowvals[MGRID+MGRID+2] = 2*MGRID+1;
516 /* middle row values */
517 for(i=0;i<(MGRID-4);i++) rowvals[MGRID+MGRID+3+4*i] = MGRID+1+i;
518 for(i=0;i<(MGRID-4);i++) rowvals[MGRID+MGRID+4+4*i] = MGRID+2+i;
519 for(i=0;i<(MGRID-4);i++) rowvals[MGRID+MGRID+5+4*i] = MGRID+3+i;
520 for(i=0;i<(MGRID-4);i++) rowvals[MGRID+MGRID+6+4*i] = 2*MGRID+2+i;
521 rowvals[2*MGRID+4*(MGRID-2)-5] = MGRID+(MGRID-2)-1;
522 rowvals[2*MGRID+4*(MGRID-2)-4] = MGRID+(MGRID-2); /* starting from here, add two diag patterns */
523 rowvals[2*MGRID+4*(MGRID-2)-3] = 2*MGRID+(MGRID-2);
524 rowvals[2*MGRID+4*(MGRID-2)-2] = MGRID+(MGRID-2);
525 rowvals[2*MGRID+4*(MGRID-2)-1] = MGRID+(MGRID-2)+1;
526
527 /**** repeated (MGRID-4 times) middle column blocks ****/
528 repeat=0;
529 for(i=0;i<MGRID-4;i++){
530 rowvals[2*MGRID+4*(MGRID-2)+repeat] = MGRID+(MGRID-2)+2+MGRID*i;
531 rowvals[2*MGRID+4*(MGRID-2)+repeat+1] = MGRID+(MGRID-2)+2+MGRID*i+1;
532
533 rowvals[2*MGRID+4*(MGRID-2)+repeat+2] = MGRID+(MGRID-2)+2+MGRID*i+1-MGRID;
534 rowvals[2*MGRID+4*(MGRID-2)+repeat+3] = MGRID+(MGRID-2)+2+MGRID*i+1;
535 rowvals[2*MGRID+4*(MGRID-2)+repeat+4] = MGRID+(MGRID-2)+2+MGRID*i+2; /* *this */
536 rowvals[2*MGRID+4*(MGRID-2)+repeat+5] = MGRID+(MGRID-2)+2+MGRID*i+1+MGRID;
537
538 /* 5 in 5*j chosen since there are 5 elements in each column */
539 /* column repeats MGRID-4 times within the outer loop */
540 for(j=0;j<MGRID-4;j++){
541 rowvals[2*MGRID+4*(MGRID-2)+repeat+6+5*j] = MGRID+(MGRID-2)+2+MGRID*i+1-MGRID+1+j;
542 rowvals[2*MGRID+4*(MGRID-2)+repeat+7+5*j] = MGRID+(MGRID-2)+2+MGRID*i+1+j;
543 rowvals[2*MGRID+4*(MGRID-2)+repeat+8+5*j] = MGRID+(MGRID-2)+2+MGRID*i+2+j;
544 rowvals[2*MGRID+4*(MGRID-2)+repeat+9+5*j] = MGRID+(MGRID-2)+2+MGRID*i+2+1+j;
545 rowvals[2*MGRID+4*(MGRID-2)+repeat+10+5*j] = MGRID+(MGRID-2)+2+MGRID*i+1+MGRID+1+j;
546 }
547
548 rowvals[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+6] = MGRID+(MGRID-2)+2+MGRID*i-2;
549 rowvals[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+7] = MGRID+(MGRID-2)+2+MGRID*i-2+MGRID-1;
550 rowvals[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+8] = MGRID+(MGRID-2)+2+MGRID*i-2+MGRID; /* *this+MGRID */
551 rowvals[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+9] = MGRID+(MGRID-2)+2+MGRID*i-2+2*MGRID;
552
553 rowvals[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+10] = MGRID+(MGRID-2)+2+MGRID*i-2+MGRID;
554 rowvals[2*MGRID+4*(MGRID-2)+repeat+(MGRID-4)*5+11] = MGRID+(MGRID-2)+2+MGRID*i-2+MGRID+1;
555
556 repeat+=MGRID+4*(MGRID-2); /* shift that accounts for accumulated columns and elements */
557 }
558
559
560 /**** last-1 column block ****/
561 rowvals[TOTAL-6*(MGRID-2)-4] = MGRID*MGRID-1-2*(MGRID-1)-1;
562 rowvals[TOTAL-6*(MGRID-2)-3] = MGRID*MGRID-1-2*(MGRID-1); /* starting with this as base */
563 rowvals[TOTAL-6*(MGRID-2)-2] = MGRID*MGRID-1-2*(MGRID-1)-MGRID;
564 rowvals[TOTAL-6*(MGRID-2)-1] = MGRID*MGRID-1-2*(MGRID-1);
565 rowvals[TOTAL-6*(MGRID-2) ] = MGRID*MGRID-1-2*(MGRID-1)+1;
566 /* middle row values */
567 for(i=0;i<(MGRID-4);i++) rowvals[TOTAL-6*(MGRID-2)+1+4*i] = MGRID*MGRID-1-2*(MGRID-1)-MGRID+1+i;
568 for(i=0;i<(MGRID-4);i++) rowvals[TOTAL-6*(MGRID-2)+2+4*i] = MGRID*MGRID-1-2*(MGRID-1)+i;
569 for(i=0;i<(MGRID-4);i++) rowvals[TOTAL-6*(MGRID-2)+3+4*i] = MGRID*MGRID-1-2*(MGRID-1)+1+i;/*copied above*/
570 for(i=0;i<(MGRID-4);i++) rowvals[TOTAL-6*(MGRID-2)+4+4*i] = MGRID*MGRID-1-2*(MGRID-1)+2+i;
571 rowvals[TOTAL-2*(MGRID-2)-7] = MGRID*MGRID-2*MGRID-2;
572 rowvals[TOTAL-2*(MGRID-2)-6] = MGRID*MGRID-MGRID-3;
573 rowvals[TOTAL-2*(MGRID-2)-5] = MGRID*MGRID-MGRID-2;
574 rowvals[TOTAL-2*(MGRID-2)-4] = MGRID*MGRID-MGRID-2;
575 rowvals[TOTAL-2*(MGRID-2)-3] = MGRID*MGRID-MGRID-1;
576
577 /* last column block */
578 rowvals[TOTAL-2*(MGRID-2)-2] = MGRID*MGRID-MGRID;
579 /* alternating pattern in data, separate loop for each pattern */
580 for(i=0;i<(MGRID-2);i++) rowvals[TOTAL-2*(MGRID-2)-1+2*i] = MGRID*MGRID-2*MGRID+1+i;
581 for(i=0;i<(MGRID-2);i++) rowvals[TOTAL-2*(MGRID-2) +2*i] = MGRID*MGRID-MGRID+1+i;
582 rowvals[TOTAL-1] = MGRID*MGRID-1;
583
584 return(0);
585 }
586
587
588 /*
589 *--------------------------------------------------------------------
590 * PRIVATE FUNCTIONS
591 *--------------------------------------------------------------------
592 */
593
594 /*
595 * SetInitialProfile: routine to initialize u, up, and id vectors.
596 */
597
SetInitialProfile(UserData data,N_Vector uu,N_Vector up,N_Vector id,N_Vector res)598 static int SetInitialProfile(UserData data, N_Vector uu, N_Vector up,
599 N_Vector id, N_Vector res)
600 {
601 realtype xfact, yfact, *udata, *updata, *iddata;
602 sunindextype mm, mm1, i, j, offset, loc;
603
604 mm = data->mm;
605 mm1 = mm - 1;
606
607 udata = N_VGetArrayPointer(uu);
608 updata = N_VGetArrayPointer(up);
609 iddata = N_VGetArrayPointer(id);
610
611 /* Initialize id to 1's. */
612 N_VConst(ONE, id);
613
614 /* Initialize uu on all grid points. */
615 for (j = 0; j < mm; j++) {
616 yfact = data->dx * j;
617 offset = mm*j;
618 for (i = 0;i < mm; i++) {
619 xfact = data->dx * i;
620 loc = offset + i;
621 udata[loc] = RCONST(16.0) * xfact * (ONE - xfact) * yfact * (ONE - yfact);
622 }
623 }
624
625 /* Initialize up vector to 0. */
626 N_VConst(ZERO, up);
627
628 /* heatres sets res to negative of ODE RHS values at interior points. */
629 heatres(ZERO, uu, up, res, data);
630
631 /* Copy -res into up to get correct interior initial up values. */
632 N_VScale(-ONE, res, up);
633
634 /* Finally, set values of u, up, and id at boundary points. */
635 for (j = 0; j < mm; j++) {
636 offset = mm*j;
637 for (i = 0;i < mm; i++) {
638 loc = offset + i;
639 if (j == 0 || j == mm1 || i == 0 || i == mm1 ) {
640 udata[loc] = BVAL; updata[loc] = ZERO; iddata[loc] = ZERO; }
641 }
642 }
643
644 return(0);
645
646 }
647
648 /*
649 * Print first lines of output (problem description)
650 */
651
PrintHeader(realtype rtol,realtype atol)652 static void PrintHeader(realtype rtol, realtype atol)
653 {
654 printf("\nidaHeat2D_klu: Heat equation, serial example problem for IDA\n");
655 printf(" Discretized heat equation on 2D unit square.\n");
656 printf(" Zero boundary conditions,");
657 printf(" polynomial initial conditions.\n");
658 printf(" Mesh dimensions: %d x %d", MGRID, MGRID);
659 printf(" Total system size: %d\n\n", NEQ);
660 #if defined(SUNDIALS_EXTENDED_PRECISION)
661 printf("Tolerance parameters: rtol = %Lg atol = %Lg\n", rtol, atol);
662 #elif defined(SUNDIALS_DOUBLE_PRECISION)
663 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
664 #else
665 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
666 #endif
667 printf("Constraints set to force all solution components >= 0. \n");
668 printf("Linear solver: KLU, sparse direct solver \n");
669 printf(" difference quotient Jacobian\n");
670 #if defined(SUNDIALS_EXTENDED_PRECISION)
671 printf("IDACalcIC called with input boundary values = %Lg \n",BVAL);
672 #elif defined(SUNDIALS_DOUBLE_PRECISION)
673 printf("IDACalcIC called with input boundary values = %g \n",BVAL);
674 #else
675 printf("IDACalcIC called with input boundary values = %g \n",BVAL);
676 #endif
677 /* Print output table heading and initial line of table. */
678 printf("\n Output Summary (umax = max-norm of solution) \n\n");
679 printf(" time umax k nst nni nje nre h \n" );
680 printf(" . . . . . . . . . . . . . . . . . . . \n");
681 }
682
683 /*
684 * Print Output
685 */
686
PrintOutput(void * mem,realtype t,N_Vector uu)687 static void PrintOutput(void *mem, realtype t, N_Vector uu)
688 {
689 int retval;
690 realtype umax, hused;
691 long int nst, nni, nje, nre;
692 int kused;
693
694 umax = N_VMaxNorm(uu);
695
696 retval = IDAGetLastOrder(mem, &kused);
697 check_retval(&retval, "IDAGetLastOrder", 1);
698 retval = IDAGetNumSteps(mem, &nst);
699 check_retval(&retval, "IDAGetNumSteps", 1);
700 retval = IDAGetNumNonlinSolvIters(mem, &nni);
701 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
702 retval = IDAGetNumResEvals(mem, &nre);
703 check_retval(&retval, "IDAGetNumResEvals", 1);
704 retval = IDAGetLastStep(mem, &hused);
705 check_retval(&retval, "IDAGetLastStep", 1);
706 retval = IDAGetNumJacEvals(mem, &nje);
707 check_retval(&retval, "IDAGetNumJacEvals", 1);
708
709 #if defined(SUNDIALS_EXTENDED_PRECISION)
710 printf(" %5.2Lf %13.5Le %d %3ld %3ld %3ld %4ld %9.2Le \n",
711 t, umax, kused, nst, nni, nje, nre, hused);
712 #elif defined(SUNDIALS_DOUBLE_PRECISION)
713 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %9.2e \n",
714 t, umax, kused, nst, nni, nje, nre, hused);
715 #else
716 printf(" %5.2f %13.5e %d %3ld %3ld %3ld %4ld %9.2e \n",
717 t, umax, kused, nst, nni, nje, nre, hused);
718 #endif
719
720 }
721
722 /*
723 * Check function return value...
724 * opt == 0 means SUNDIALS function allocates memory so check if
725 * returned NULL pointer
726 * opt == 1 means SUNDIALS function returns an integer value so check if
727 * retval < 0
728 * opt == 2 means function allocates memory so check if returned
729 * NULL pointer
730 */
731
check_retval(void * returnvalue,const char * funcname,int opt)732 static int check_retval(void *returnvalue, const char *funcname, int opt)
733 {
734 int *retval;
735
736 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
737 if (opt == 0 && returnvalue == NULL) {
738 fprintf(stderr,
739 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
740 funcname);
741 return(1);
742 } else if (opt == 1) {
743 /* Check if retval < 0 */
744 retval = (int *) returnvalue;
745 if (*retval < 0) {
746 fprintf(stderr,
747 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
748 funcname, *retval);
749 return(1);
750 }
751 } else if (opt == 2 && returnvalue == NULL) {
752 /* Check if function returned NULL pointer - no memory allocated */
753 fprintf(stderr,
754 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
755 funcname);
756 return(1);
757 }
758
759 return(0);
760 }
761