Actual source code: ex31.c

  1: static char help[] = "Solves the ordinary differential equations (IVPs) using explicit and implicit time-integration methods.\n";

  3: /*

  5:   Concepts:   TS
  6:   Useful command line parameters:
  7:   -problem <hull1972a1>: choose which problem to solve (see references
  8:                       for complete listing of problems).
  9:   -ts_type <euler>: specify time-integrator
 10:   -ts_adapt_type <basic>: specify time-step adapting (none,basic,advanced)
 11:   -refinement_levels <1>: number of refinement levels for convergence analysis
 12:   -refinement_factor <2.0>: factor to refine time step size by for convergence analysis
 13:   -dt <0.01>: specify time step (initial time step for convergence analysis)

 15: */

 17: /*
 18: List of cases and their names in the code:-
 19:   From Hull, T.E., Enright, W.H., Fellen, B.M., and Sedgwick, A.E.,
 20:       "Comparing Numerical Methods for Ordinary Differential
 21:        Equations", SIAM J. Numer. Anal., 9(4), 1972, pp. 603 - 635
 22:     A1 -> "hull1972a1" (exact solution available)
 23:     A2 -> "hull1972a2" (exact solution available)
 24:     A3 -> "hull1972a3" (exact solution available)
 25:     A4 -> "hull1972a4" (exact solution available)
 26:     A5 -> "hull1972a5"
 27:     B1 -> "hull1972b1"
 28:     B2 -> "hull1972b2"
 29:     B3 -> "hull1972b3"
 30:     B4 -> "hull1972b4"
 31:     B5 -> "hull1972b5"
 32:     C1 -> "hull1972c1"
 33:     C2 -> "hull1972c2"
 34:     C3 -> "hull1972c3"
 35:     C4 -> "hull1972c4"

 37:  From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
 38:        https://arxiv.org/abs/1503.05166, 2016

 40:     Kulikov2013I -> "kulik2013i"

 42: */

 44: #include <petscts.h>

 46: /* Function declarations */
 47: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
 48: PetscErrorCode (*RHSJacobian) (TS,PetscReal,Vec,Mat,Mat,void*);
 49: PetscErrorCode (*IFunction)   (TS,PetscReal,Vec,Vec,Vec,void*);
 50: PetscErrorCode (*IJacobian)   (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

 52: /* Returns the size of the system of equations depending on problem specification */
 53: PetscInt GetSize(const char *p)
 54: {
 55:   if      ((!strcmp(p,"hull1972a1"))
 56:          ||(!strcmp(p,"hull1972a2"))
 57:          ||(!strcmp(p,"hull1972a3"))
 58:          ||(!strcmp(p,"hull1972a4"))
 59:          ||(!strcmp(p,"hull1972a5"))) return 1;
 60:   else if  (!strcmp(p,"hull1972b1")) return 2;
 61:   else if ((!strcmp(p,"hull1972b2"))
 62:          ||(!strcmp(p,"hull1972b3"))
 63:          ||(!strcmp(p,"hull1972b4"))
 64:          ||(!strcmp(p,"hull1972b5"))) return 3;
 65:   else if ((!strcmp(p,"kulik2013i"))) return 4;
 66:   else if ((!strcmp(p,"hull1972c1"))
 67:          ||(!strcmp(p,"hull1972c2"))
 68:          ||(!strcmp(p,"hull1972c3"))) return 10;
 69:   else if  (!strcmp(p,"hull1972c4")) return 51;
 70:   else PetscFunctionReturn(-1);
 71: }

 73: /****************************************************************/

 75: /* Problem specific functions */

 77: /* Hull, 1972, Problem A1 */

 79: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
 80: {
 81:   PetscScalar       *f;
 82:   const PetscScalar *y;

 84:   VecGetArrayRead(Y,&y);
 85:   VecGetArray(F,&f);
 86:   f[0] = -y[0];
 87:   VecRestoreArrayRead(Y,&y);
 88:   VecRestoreArray(F,&f);
 89:   return 0;
 90: }

 92: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
 93: {
 94:   const PetscScalar *y;
 95:   PetscInt          row = 0,col = 0;
 96:   PetscScalar       value = -1.0;

 98:   VecGetArrayRead(Y,&y);
 99:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
100:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
102:   VecRestoreArrayRead(Y,&y);
103:   return 0;
104: }

106: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
107: {
108:   const PetscScalar *y;
109:   PetscScalar       *f;

111:   VecGetArrayRead(Y,&y);
112:   VecGetArray(F,&f);
113:   f[0] = -y[0];
114:   VecRestoreArrayRead(Y,&y);
115:   VecRestoreArray(F,&f);
116:   /* Left hand side = ydot - f(y) */
117:   VecAYPX(F,-1.0,Ydot);
118:   return 0;
119: }

121: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
122: {
123:   const PetscScalar *y;
124:   PetscInt          row = 0,col = 0;
125:   PetscScalar       value = a - 1.0;

127:   VecGetArrayRead(Y,&y);
128:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
129:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
131:   VecRestoreArrayRead(Y,&y);
132:   return 0;
133: }

135: /* Hull, 1972, Problem A2 */

137: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
138: {
139:   const PetscScalar *y;
140:   PetscScalar       *f;

142:   VecGetArrayRead(Y,&y);
143:   VecGetArray(F,&f);
144:   f[0] = -0.5*y[0]*y[0]*y[0];
145:   VecRestoreArrayRead(Y,&y);
146:   VecRestoreArray(F,&f);
147:   return 0;
148: }

150: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
151: {
152:   const PetscScalar *y;
153:   PetscInt          row = 0,col = 0;
154:   PetscScalar       value;

156:   VecGetArrayRead(Y,&y);
157:   value = -0.5*3.0*y[0]*y[0];
158:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
159:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
160:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
161:   VecRestoreArrayRead(Y,&y);
162:   return 0;
163: }

165: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
166: {
167:   PetscScalar       *f;
168:   const PetscScalar *y;

170:   VecGetArrayRead(Y,&y);
171:   VecGetArray(F,&f);
172:   f[0] = -0.5*y[0]*y[0]*y[0];
173:   VecRestoreArrayRead(Y,&y);
174:   VecRestoreArray(F,&f);
175:   /* Left hand side = ydot - f(y) */
176:   VecAYPX(F,-1.0,Ydot);
177:   return 0;
178: }

180: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
181: {
182:   const PetscScalar *y;
183:   PetscInt          row = 0,col = 0;
184:   PetscScalar       value;

186:   VecGetArrayRead(Y,&y);
187:   value = a + 0.5*3.0*y[0]*y[0];
188:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
189:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
190:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
191:   VecRestoreArrayRead(Y,&y);
192:   return 0;
193: }

195: /* Hull, 1972, Problem A3 */

197: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
198: {
199:   const PetscScalar *y;
200:   PetscScalar       *f;

202:   VecGetArrayRead(Y,&y);
203:   VecGetArray(F,&f);
204:   f[0] = y[0]*PetscCosReal(t);
205:   VecRestoreArrayRead(Y,&y);
206:   VecRestoreArray(F,&f);
207:   return 0;
208: }

210: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
211: {
212:   const PetscScalar *y;
213:   PetscInt          row = 0,col = 0;
214:   PetscScalar       value = PetscCosReal(t);

216:   VecGetArrayRead(Y,&y);
217:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
218:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
219:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
220:   VecRestoreArrayRead(Y,&y);
221:   return 0;
222: }

224: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
225: {
226:   PetscScalar       *f;
227:   const PetscScalar *y;

229:   VecGetArrayRead(Y,&y);
230:   VecGetArray(F,&f);
231:   f[0] = y[0]*PetscCosReal(t);
232:   VecRestoreArrayRead(Y,&y);
233:   VecRestoreArray(F,&f);
234:   /* Left hand side = ydot - f(y) */
235:   VecAYPX(F,-1.0,Ydot);
236:   return 0;
237: }

239: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
240: {
241:   const PetscScalar *y;
242:   PetscInt          row = 0,col = 0;
243:   PetscScalar       value = a - PetscCosReal(t);

245:   VecGetArrayRead(Y,&y);
246:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
247:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
248:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
249:   VecRestoreArrayRead(Y,&y);
250:   return 0;
251: }

253: /* Hull, 1972, Problem A4 */

255: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
256: {
257:   PetscScalar       *f;
258:   const PetscScalar *y;

260:   VecGetArrayRead(Y,&y);
261:   VecGetArray(F,&f);
262:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
263:   VecRestoreArrayRead(Y,&y);
264:   VecRestoreArray(F,&f);
265:   return 0;
266: }

268: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
269: {
270:   const PetscScalar *y;
271:   PetscInt          row = 0,col = 0;
272:   PetscScalar       value;

274:   VecGetArrayRead(Y,&y);
275:   value = 0.25*(1.0-0.05*y[0]) - (0.25*y[0])*0.05;
276:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
277:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
278:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
279:   VecRestoreArrayRead(Y,&y);
280:   return 0;
281: }

283: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
284: {
285:   PetscScalar       *f;
286:   const PetscScalar *y;

288:   VecGetArrayRead(Y,&y);
289:   VecGetArray(F,&f);
290:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
291:   VecRestoreArrayRead(Y,&y);
292:   VecRestoreArray(F,&f);
293:   /* Left hand side = ydot - f(y) */
294:   VecAYPX(F,-1.0,Ydot);
295:   return 0;
296: }

298: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
299: {
300:   const PetscScalar *y;
301:   PetscInt          row = 0,col = 0;
302:   PetscScalar       value;

304:   VecGetArrayRead(Y,&y);
305:   value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
306:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
307:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
308:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
309:   VecRestoreArrayRead(Y,&y);
310:   return 0;
311: }

313: /* Hull, 1972, Problem A5 */

315: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
316: {
317:   PetscScalar       *f;
318:   const PetscScalar *y;

320:   VecGetArrayRead(Y,&y);
321:   VecGetArray(F,&f);
322:   f[0] = (y[0]-t)/(y[0]+t);
323:   VecRestoreArrayRead(Y,&y);
324:   VecRestoreArray(F,&f);
325:   return 0;
326: }

328: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
329: {
330:   const PetscScalar *y;
331:   PetscInt          row = 0,col = 0;
332:   PetscScalar       value;

334:   VecGetArrayRead(Y,&y);
335:   value = 2*t/((t+y[0])*(t+y[0]));
336:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
337:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
338:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
339:   VecRestoreArrayRead(Y,&y);
340:   return 0;
341: }

343: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
344: {
345:   PetscScalar       *f;
346:   const PetscScalar *y;

348:   VecGetArrayRead(Y,&y);
349:   VecGetArray(F,&f);
350:   f[0] = (y[0]-t)/(y[0]+t);
351:   VecRestoreArrayRead(Y,&y);
352:   VecRestoreArray(F,&f);
353:   /* Left hand side = ydot - f(y) */
354:   VecAYPX(F,-1.0,Ydot);
355:   return 0;
356: }

358: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
359: {
360:   const PetscScalar *y;
361:   PetscInt          row = 0,col = 0;
362:   PetscScalar       value;

364:   VecGetArrayRead(Y,&y);
365:   value = a - 2*t/((t+y[0])*(t+y[0]));
366:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
367:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
368:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
369:   VecRestoreArrayRead(Y,&y);
370:   return 0;
371: }

373: /* Hull, 1972, Problem B1 */

375: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
376: {
377:   PetscScalar       *f;
378:   const PetscScalar *y;

380:   VecGetArrayRead(Y,&y);
381:   VecGetArray(F,&f);
382:   f[0] = 2.0*(y[0] - y[0]*y[1]);
383:   f[1] = -(y[1]-y[0]*y[1]);
384:   VecRestoreArrayRead(Y,&y);
385:   VecRestoreArray(F,&f);
386:   return 0;
387: }

389: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
390: {
391:   PetscScalar       *f;
392:   const PetscScalar *y;

394:   VecGetArrayRead(Y,&y);
395:   VecGetArray(F,&f);
396:   f[0] = 2.0*(y[0] - y[0]*y[1]);
397:   f[1] = -(y[1]-y[0]*y[1]);
398:   VecRestoreArrayRead(Y,&y);
399:   VecRestoreArray(F,&f);
400:   /* Left hand side = ydot - f(y) */
401:   VecAYPX(F,-1.0,Ydot);
402:   return 0;
403: }

405: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
406: {
407:   const PetscScalar *y;
408:   PetscInt          row[2] = {0,1};
409:   PetscScalar       value[2][2];

411:   VecGetArrayRead(Y,&y);
412:   value[0][0] = a - 2.0*(1.0-y[1]);    value[0][1] = 2.0*y[0];
413:   value[1][0] = -y[1];                 value[1][1] = a + 1.0 - y[0];
414:   MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
415:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
416:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
417:   VecRestoreArrayRead(Y,&y);
418:   return 0;
419: }

421: /* Hull, 1972, Problem B2 */

423: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
424: {
425:   PetscScalar       *f;
426:   const PetscScalar *y;

428:   VecGetArrayRead(Y,&y);
429:   VecGetArray(F,&f);
430:   f[0] = -y[0] + y[1];
431:   f[1] = y[0] - 2.0*y[1] + y[2];
432:   f[2] = y[1] - y[2];
433:   VecRestoreArrayRead(Y,&y);
434:   VecRestoreArray(F,&f);
435:   return 0;
436: }

438: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
439: {
440:   PetscScalar       *f;
441:   const PetscScalar *y;

443:   VecGetArrayRead(Y,&y);
444:   VecGetArray(F,&f);
445:   f[0] = -y[0] + y[1];
446:   f[1] = y[0] - 2.0*y[1] + y[2];
447:   f[2] = y[1] - y[2];
448:   VecRestoreArrayRead(Y,&y);
449:   VecRestoreArray(F,&f);
450:   /* Left hand side = ydot - f(y) */
451:   VecAYPX(F,-1.0,Ydot);
452:   return 0;
453: }

455: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
456: {
457:   const PetscScalar *y;
458:   PetscInt          row[3] = {0,1,2};
459:   PetscScalar       value[3][3];

461:   VecGetArrayRead(Y,&y);
462:   value[0][0] = a + 1.0;  value[0][1] = -1.0;     value[0][2] = 0;
463:   value[1][0] = -1.0;     value[1][1] = a + 2.0;  value[1][2] = -1.0;
464:   value[2][0] = 0;        value[2][1] = -1.0;     value[2][2] = a + 1.0;
465:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
466:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
467:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
468:   VecRestoreArrayRead(Y,&y);
469:   return 0;
470: }

472: /* Hull, 1972, Problem B3 */

474: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
475: {
476:   PetscScalar       *f;
477:   const PetscScalar *y;

479:   VecGetArrayRead(Y,&y);
480:   VecGetArray(F,&f);
481:   f[0] = -y[0];
482:   f[1] = y[0] - y[1]*y[1];
483:   f[2] = y[1]*y[1];
484:   VecRestoreArrayRead(Y,&y);
485:   VecRestoreArray(F,&f);
486:   return 0;
487: }

489: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
490: {
491:   PetscScalar       *f;
492:   const PetscScalar *y;

494:   VecGetArrayRead(Y,&y);
495:   VecGetArray(F,&f);
496:   f[0] = -y[0];
497:   f[1] = y[0] - y[1]*y[1];
498:   f[2] = y[1]*y[1];
499:   VecRestoreArrayRead(Y,&y);
500:   VecRestoreArray(F,&f);
501:   /* Left hand side = ydot - f(y) */
502:   VecAYPX(F,-1.0,Ydot);
503:   return 0;
504: }

506: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
507: {
508:   const PetscScalar *y;
509:   PetscInt          row[3] = {0,1,2};
510:   PetscScalar       value[3][3];

512:   VecGetArrayRead(Y,&y);
513:   value[0][0] = a + 1.0; value[0][1] = 0;             value[0][2] = 0;
514:   value[1][0] = -1.0;    value[1][1] = a + 2.0*y[1];  value[1][2] = 0;
515:   value[2][0] = 0;       value[2][1] = -2.0*y[1];     value[2][2] = a;
516:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
517:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
518:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
519:   VecRestoreArrayRead(Y,&y);
520:   return 0;
521: }

523: /* Hull, 1972, Problem B4 */

525: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
526: {
527:   PetscScalar       *f;
528:   const PetscScalar *y;

530:   VecGetArrayRead(Y,&y);
531:   VecGetArray(F,&f);
532:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
533:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
534:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
535:   VecRestoreArrayRead(Y,&y);
536:   VecRestoreArray(F,&f);
537:   return 0;
538: }

540: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
541: {
542:   PetscScalar       *f;
543:   const PetscScalar *y;

545:   VecGetArrayRead(Y,&y);
546:   VecGetArray(F,&f);
547:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
548:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
549:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
550:   VecRestoreArrayRead(Y,&y);
551:   VecRestoreArray(F,&f);
552:   /* Left hand side = ydot - f(y) */
553:   VecAYPX(F,-1.0,Ydot);
554:   return 0;
555: }

557: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
558: {
559:   const PetscScalar *y;
560:   PetscInt          row[3] = {0,1,2};
561:   PetscScalar       value[3][3],fac,fac2;

563:   VecGetArrayRead(Y,&y);
564:   fac  = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
565:   fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
566:   value[0][0] = a + (y[1]*y[1]*y[2])*fac;
567:   value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
568:   value[0][2] = y[0]*fac2;
569:   value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
570:   value[1][1] = a + y[0]*y[0]*y[2]*fac;
571:   value[1][2] = y[1]*fac2;
572:   value[2][0] = -y[1]*y[1]*fac;
573:   value[2][1] = y[0]*y[1]*fac;
574:   value[2][2] = a;
575:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
576:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
577:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
578:   VecRestoreArrayRead(Y,&y);
579:   return 0;
580: }

582: /* Hull, 1972, Problem B5 */

584: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
585: {
586:   PetscScalar       *f;
587:   const PetscScalar *y;

589:   VecGetArrayRead(Y,&y);
590:   VecGetArray(F,&f);
591:   f[0] = y[1]*y[2];
592:   f[1] = -y[0]*y[2];
593:   f[2] = -0.51*y[0]*y[1];
594:   VecRestoreArrayRead(Y,&y);
595:   VecRestoreArray(F,&f);
596:   return 0;
597: }

599: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
600: {
601:   PetscScalar       *f;
602:   const PetscScalar *y;

604:   VecGetArrayRead(Y,&y);
605:   VecGetArray(F,&f);
606:   f[0] = y[1]*y[2];
607:   f[1] = -y[0]*y[2];
608:   f[2] = -0.51*y[0]*y[1];
609:   VecRestoreArrayRead(Y,&y);
610:   VecRestoreArray(F,&f);
611:   /* Left hand side = ydot - f(y) */
612:   VecAYPX(F,-1.0,Ydot);
613:   return 0;
614: }

616: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
617: {
618:   const PetscScalar *y;
619:   PetscInt          row[3] = {0,1,2};
620:   PetscScalar       value[3][3];

622:   VecGetArrayRead(Y,&y);
623:   value[0][0] = a;          value[0][1] = -y[2];      value[0][2] = -y[1];
624:   value[1][0] = y[2];       value[1][1] = a;          value[1][2] = y[0];
625:   value[2][0] = 0.51*y[1];  value[2][1] = 0.51*y[0];  value[2][2] = a;
626:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
627:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
628:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
629:   VecRestoreArrayRead(Y,&y);
630:   return 0;
631: }

633: /* Kulikov, 2013, Problem I */

635: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
636: {
637:   PetscScalar       *f;
638:   const PetscScalar *y;

640:   VecGetArrayRead(Y,&y);
641:   VecGetArray(F,&f);
642:   f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
643:   f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
644:   f[2] = 2.*t*y[3];
645:   f[3] = -2.*t*PetscLogScalar(y[0]);
646:   VecRestoreArrayRead(Y,&y);
647:   VecRestoreArray(F,&f);
648:   return 0;
649: }

651: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
652: {
653:   const PetscScalar *y;
654:   PetscInt          row[4] = {0,1,2,3};
655:   PetscScalar       value[4][4];
656:   PetscScalar       m1,m2;
657:   VecGetArrayRead(Y,&y);
658:   m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
659:   m2=2.*t*PetscPowScalar(y[1],1./5.);
660:   value[0][0] = 0. ;        value[0][1] = m1; value[0][2] = 0.;  value[0][3] = m2;
661:   m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
662:   m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
663:   value[1][0] = 0.;        value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
664:   value[2][0] = 0.;        value[2][1] = 0.;  value[2][2] = 0.; value[2][3] = 2*t;
665:   value[3][0] = -2.*t/y[0];value[3][1] = 0.;  value[3][2] = 0.; value[3][3] = 0.;
666:   MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
667:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
668:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
669:   VecRestoreArrayRead(Y,&y);
670:   return 0;
671: }

673: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
674: {
675:   PetscScalar       *f;
676:   const PetscScalar *y;

678:   VecGetArrayRead(Y,&y);
679:   VecGetArray(F,&f);
680:   f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
681:   f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
682:   f[2] = 2.*t*y[3];
683:   f[3] = -2.*t*PetscLogScalar(y[0]);
684:   VecRestoreArrayRead(Y,&y);
685:   VecRestoreArray(F,&f);
686:   /* Left hand side = ydot - f(y) */
687:   VecAYPX(F,-1.0,Ydot);
688:   return 0;
689: }

691: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
692: {
693:   const PetscScalar *y;
694:   PetscInt          row[4] = {0,1,2,3};
695:   PetscScalar       value[4][4];
696:   PetscScalar       m1,m2;

698:   VecGetArrayRead(Y,&y);
699:   m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
700:   m2=2.*t*PetscPowScalar(y[1],1./5.);
701:   value[0][0] = a ;        value[0][1] = m1;  value[0][2] = 0.; value[0][3] = m2;
702:   m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
703:   m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
704:   value[1][0] = 0.;        value[1][1] = a ;  value[1][2] = m1; value[1][3] = m2;
705:   value[2][0] = 0.;        value[2][1] = 0.;  value[2][2] = a;  value[2][3] = 2*t;
706:   value[3][0] = -2.*t/y[0];value[3][1] = 0.;  value[3][2] = 0.; value[3][3] = a;
707:   MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
708:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
709:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
710:   VecRestoreArrayRead(Y,&y);
711:   return 0;
712: }

714: /* Hull, 1972, Problem C1 */

716: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
717: {
718:   PetscScalar       *f;
719:   const PetscScalar *y;
720:   PetscInt          N,i;

722:   VecGetSize (Y,&N);
723:   VecGetArrayRead(Y,&y);
724:   VecGetArray(F,&f);
725:   f[0] = -y[0];
726:   for (i = 1; i < N-1; i++) {
727:     f[i] = y[i-1] - y[i];
728:   }
729:   f[N-1] = y[N-2];
730:   VecRestoreArrayRead(Y,&y);
731:   VecRestoreArray(F,&f);
732:   return 0;
733: }

735: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
736: {
737:   PetscScalar       *f;
738:   const PetscScalar *y;
739:   PetscInt          N,i;

741:   VecGetSize (Y,&N);
742:   VecGetArrayRead(Y,&y);
743:   VecGetArray(F,&f);
744:   f[0] = -y[0];
745:   for (i = 1; i < N-1; i++) {
746:     f[i] = y[i-1] - y[i];
747:   }
748:   f[N-1] = y[N-2];
749:   VecRestoreArrayRead(Y,&y);
750:   VecRestoreArray(F,&f);
751:   /* Left hand side = ydot - f(y) */
752:   VecAYPX(F,-1.0,Ydot);
753:   return 0;
754: }

756: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
757: {
758:   const PetscScalar *y;
759:   PetscInt          N,i,col[2];
760:   PetscScalar       value[2];

762:   VecGetSize (Y,&N);
763:   VecGetArrayRead(Y,&y);
764:   i = 0;
765:   value[0] = a+1; col[0] = 0;
766:   value[1] =  0;  col[1] = 1;
767:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
768:   for (i = 0; i < N; i++) {
769:     value[0] =  -1; col[0] = i-1;
770:     value[1] = a+1; col[1] = i;
771:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
772:   }
773:   i = N-1;
774:   value[0] = -1;  col[0] = N-2;
775:   value[1] = a;   col[1] = N-1;
776:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
777:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
778:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
779:   VecRestoreArrayRead(Y,&y);
780:   return 0;
781: }

783: /* Hull, 1972, Problem C2 */

785: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
786: {
787:   const PetscScalar *y;
788:   PetscScalar       *f;
789:   PetscInt          N,i;

791:   VecGetSize (Y,&N);
792:   VecGetArrayRead(Y,&y);
793:   VecGetArray(F,&f);
794:   f[0] = -y[0];
795:   for (i = 1; i < N-1; i++) {
796:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
797:   }
798:   f[N-1] = (PetscReal)(N-1)*y[N-2];
799:   VecRestoreArrayRead(Y,&y);
800:   VecRestoreArray(F,&f);
801:   return 0;
802: }

804: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
805: {
806:   PetscScalar       *f;
807:   const PetscScalar *y;
808:   PetscInt          N,i;

810:   VecGetSize (Y,&N);
811:   VecGetArrayRead(Y,&y);
812:   VecGetArray(F,&f);
813:   f[0] = -y[0];
814:   for (i = 1; i < N-1; i++) {
815:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
816:   }
817:   f[N-1] = (PetscReal)(N-1)*y[N-2];
818:   VecRestoreArrayRead(Y,&y);
819:   VecRestoreArray(F,&f);
820:   /* Left hand side = ydot - f(y) */
821:   VecAYPX(F,-1.0,Ydot);
822:   return 0;
823: }

825: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
826: {
827:   const PetscScalar *y;
828:   PetscInt          N,i,col[2];
829:   PetscScalar       value[2];

831:   VecGetSize (Y,&N);
832:   VecGetArrayRead(Y,&y);
833:   i = 0;
834:   value[0] = a+1;                 col[0] = 0;
835:   value[1] = 0;                   col[1] = 1;
836:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
837:   for (i = 0; i < N; i++) {
838:     value[0] = -(PetscReal) i;      col[0] = i-1;
839:     value[1] = a+(PetscReal)(i+1);  col[1] = i;
840:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
841:   }
842:   i = N-1;
843:   value[0] = -(PetscReal) (N-1);  col[0] = N-2;
844:   value[1] = a;                   col[1] = N-1;
845:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
846:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
847:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
848:   VecRestoreArrayRead(Y,&y);
849:   return 0;
850: }

852: /* Hull, 1972, Problem C3 and C4 */

854: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
855: {
856:   PetscScalar       *f;
857:   const PetscScalar *y;
858:   PetscInt          N,i;

860:   VecGetSize (Y,&N);
861:   VecGetArrayRead(Y,&y);
862:   VecGetArray(F,&f);
863:   f[0] = -2.0*y[0] + y[1];
864:   for (i = 1; i < N-1; i++) {
865:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
866:   }
867:   f[N-1] = y[N-2] - 2.0*y[N-1];
868:   VecRestoreArrayRead(Y,&y);
869:   VecRestoreArray(F,&f);
870:   return 0;
871: }

873: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
874: {
875:   PetscScalar       *f;
876:   const PetscScalar *y;
877:   PetscInt          N,i;

879:   VecGetSize (Y,&N);
880:   VecGetArrayRead(Y,&y);
881:   VecGetArray(F,&f);
882:   f[0] = -2.0*y[0] + y[1];
883:   for (i = 1; i < N-1; i++) {
884:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
885:   }
886:   f[N-1] = y[N-2] - 2.0*y[N-1];
887:   VecRestoreArrayRead(Y,&y);
888:   VecRestoreArray(F,&f);
889:   /* Left hand side = ydot - f(y) */
890:   VecAYPX(F,-1.0,Ydot);
891:   return 0;
892: }

894: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
895: {
896:   const PetscScalar *y;
897:   PetscScalar       value[3];
898:   PetscInt          N,i,col[3];

900:   VecGetSize (Y,&N);
901:   VecGetArrayRead(Y,&y);
902:   for (i = 0; i < N; i++) {
903:     if (i == 0) {
904:       value[0] = a+2;  col[0] = i;
905:       value[1] =  -1;  col[1] = i+1;
906:       value[2] =  0;   col[2] = i+2;
907:     } else if (i == N-1) {
908:       value[0] =  0;   col[0] = i-2;
909:       value[1] =  -1;  col[1] = i-1;
910:       value[2] = a+2;  col[2] = i;
911:     } else {
912:       value[0] = -1;   col[0] = i-1;
913:       value[1] = a+2;  col[1] = i;
914:       value[2] = -1;   col[2] = i+1;
915:     }
916:     MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
917:   }
918:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
919:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
920:   VecRestoreArrayRead(Y,&y);
921:   return 0;
922: }

924: /***************************************************************************/

926: /* Sets the initial solution for the IVP and sets up the function pointers*/
927: PetscErrorCode Initialize(Vec Y, void* s)
928: {
929:   char          *p = (char*) s;
930:   PetscScalar   *y;
931:   PetscReal     t0;
932:   PetscInt      N = GetSize((const char *)s);
933:   PetscBool     flg;

935:   VecZeroEntries(Y);
936:   VecGetArray(Y,&y);
937:   if (!strcmp(p,"hull1972a1")) {
938:     y[0] = 1.0;
939:     RHSFunction = RHSFunction_Hull1972A1;
940:     RHSJacobian = RHSJacobian_Hull1972A1;
941:     IFunction   = IFunction_Hull1972A1;
942:     IJacobian   = IJacobian_Hull1972A1;
943:   } else if (!strcmp(p,"hull1972a2")) {
944:     y[0] = 1.0;
945:     RHSFunction = RHSFunction_Hull1972A2;
946:     RHSJacobian = RHSJacobian_Hull1972A2;
947:     IFunction   = IFunction_Hull1972A2;
948:     IJacobian   = IJacobian_Hull1972A2;
949:   } else if (!strcmp(p,"hull1972a3")) {
950:     y[0] = 1.0;
951:     RHSFunction = RHSFunction_Hull1972A3;
952:     RHSJacobian = RHSJacobian_Hull1972A3;
953:     IFunction   = IFunction_Hull1972A3;
954:     IJacobian   = IJacobian_Hull1972A3;
955:   } else if (!strcmp(p,"hull1972a4")) {
956:     y[0] = 1.0;
957:     RHSFunction = RHSFunction_Hull1972A4;
958:     RHSJacobian = RHSJacobian_Hull1972A4;
959:     IFunction   = IFunction_Hull1972A4;
960:     IJacobian   = IJacobian_Hull1972A4;
961:   } else if (!strcmp(p,"hull1972a5")) {
962:     y[0] = 4.0;
963:     RHSFunction = RHSFunction_Hull1972A5;
964:     RHSJacobian = RHSJacobian_Hull1972A5;
965:     IFunction   = IFunction_Hull1972A5;
966:     IJacobian   = IJacobian_Hull1972A5;
967:   } else if (!strcmp(p,"hull1972b1")) {
968:     y[0] = 1.0;
969:     y[1] = 3.0;
970:     RHSFunction = RHSFunction_Hull1972B1;
971:     IFunction   = IFunction_Hull1972B1;
972:     IJacobian   = IJacobian_Hull1972B1;
973:   } else if (!strcmp(p,"hull1972b2")) {
974:     y[0] = 2.0;
975:     y[1] = 0.0;
976:     y[2] = 1.0;
977:     RHSFunction = RHSFunction_Hull1972B2;
978:     IFunction   = IFunction_Hull1972B2;
979:     IJacobian   = IJacobian_Hull1972B2;
980:   } else if (!strcmp(p,"hull1972b3")) {
981:     y[0] = 1.0;
982:     y[1] = 0.0;
983:     y[2] = 0.0;
984:     RHSFunction = RHSFunction_Hull1972B3;
985:     IFunction   = IFunction_Hull1972B3;
986:     IJacobian   = IJacobian_Hull1972B3;
987:   } else if (!strcmp(p,"hull1972b4")) {
988:     y[0] = 3.0;
989:     y[1] = 0.0;
990:     y[2] = 0.0;
991:     RHSFunction = RHSFunction_Hull1972B4;
992:     IFunction   = IFunction_Hull1972B4;
993:     IJacobian   = IJacobian_Hull1972B4;
994:   } else if (!strcmp(p,"hull1972b5")) {
995:     y[0] = 0.0;
996:     y[1] = 1.0;
997:     y[2] = 1.0;
998:     RHSFunction = RHSFunction_Hull1972B5;
999:     IFunction   = IFunction_Hull1972B5;
1000:     IJacobian   = IJacobian_Hull1972B5;
1001:   } else if (!strcmp(p,"kulik2013i")) {
1002:     t0=0.;
1003:     y[0] = PetscExpReal(PetscSinReal(t0*t0));
1004:     y[1] = PetscExpReal(5.*PetscSinReal(t0*t0));
1005:     y[2] = PetscSinReal(t0*t0)+1.0;
1006:     y[3] = PetscCosReal(t0*t0);
1007:     RHSFunction = RHSFunction_Kulikov2013I;
1008:     RHSJacobian = RHSJacobian_Kulikov2013I;
1009:     IFunction   = IFunction_Kulikov2013I;
1010:     IJacobian   = IJacobian_Kulikov2013I;
1011:   } else if (!strcmp(p,"hull1972c1")) {
1012:     y[0] = 1.0;
1013:     RHSFunction = RHSFunction_Hull1972C1;
1014:     IFunction   = IFunction_Hull1972C1;
1015:     IJacobian   = IJacobian_Hull1972C1;
1016:   } else if (!strcmp(p,"hull1972c2")) {
1017:     y[0] = 1.0;
1018:     RHSFunction = RHSFunction_Hull1972C2;
1019:     IFunction   = IFunction_Hull1972C2;
1020:     IJacobian   = IJacobian_Hull1972C2;
1021:   } else if ((!strcmp(p,"hull1972c3"))
1022:            ||(!strcmp(p,"hull1972c4"))) {
1023:     y[0] = 1.0;
1024:     RHSFunction = RHSFunction_Hull1972C34;
1025:     IFunction   = IFunction_Hull1972C34;
1026:     IJacobian   = IJacobian_Hull1972C34;
1027:   }
1028:   PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1030:   VecRestoreArray(Y,&y);
1031:   return 0;
1032: }

1034: /* Calculates the exact solution to problems that have one */
1035: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1036: {
1037:   char          *p = (char*) s;
1038:   PetscScalar   *y;

1040:   if (!strcmp(p,"hull1972a1")) {
1041:     VecGetArray(Y,&y);
1042:     y[0] = PetscExpReal(-t);
1043:     *flag = PETSC_TRUE;
1044:     VecRestoreArray(Y,&y);
1045:   } else if (!strcmp(p,"hull1972a2")) {
1046:     VecGetArray(Y,&y);
1047:     y[0] = 1.0/PetscSqrtReal(t+1);
1048:     *flag = PETSC_TRUE;
1049:     VecRestoreArray(Y,&y);
1050:   } else if (!strcmp(p,"hull1972a3")) {
1051:     VecGetArray(Y,&y);
1052:     y[0] = PetscExpReal(PetscSinReal(t));
1053:     *flag = PETSC_TRUE;
1054:     VecRestoreArray(Y,&y);
1055:   } else if (!strcmp(p,"hull1972a4")) {
1056:     VecGetArray(Y,&y);
1057:     y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1058:     *flag = PETSC_TRUE;
1059:     VecRestoreArray(Y,&y);
1060:   } else if (!strcmp(p,"kulik2013i")) {
1061:     VecGetArray(Y,&y);
1062:     y[0] = PetscExpReal(PetscSinReal(t*t));
1063:     y[1] = PetscExpReal(5.*PetscSinReal(t*t));
1064:     y[2] = PetscSinReal(t*t)+1.0;
1065:     y[3] = PetscCosReal(t*t);
1066:     *flag = PETSC_TRUE;
1067:     VecRestoreArray(Y,&y);
1068:   } else {
1069:     VecSet(Y,0);
1070:     *flag = PETSC_FALSE;
1071:   }
1072:   return 0;
1073: }

1075: /* Solves the specified ODE and computes the error if exact solution is available */
1076: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1077: {
1078:   TS              ts;               /* time-integrator                        */
1079:   Vec             Y;                /* Solution vector                        */
1080:   Vec             Yex;              /* Exact solution                         */
1081:   PetscInt        N;                /* Size of the system of equations        */
1082:   TSType          time_scheme;      /* Type of time-integration scheme        */
1083:   Mat             Jac = NULL;       /* Jacobian matrix                        */
1084:   Vec             Yerr;             /* Auxiliary solution vector              */
1085:   PetscReal       err_norm;         /* Estimated error norm                   */
1086:   PetscReal       final_time;       /* Actual final time from the integrator  */

1088:   N = GetSize((const char *)&ptype[0]);
1090:   VecCreate(PETSC_COMM_WORLD,&Y);
1091:   VecSetSizes(Y,N,PETSC_DECIDE);
1092:   VecSetUp(Y);
1093:   VecSet(Y,0);

1095:   /* Initialize the problem */
1096:   Initialize(Y,&ptype[0]);

1098:   /* Create and initialize the time-integrator                            */
1099:   TSCreate(PETSC_COMM_WORLD,&ts);
1100:   /* Default time integration options                                     */
1101:   TSSetType(ts,TSRK);
1102:   TSSetMaxSteps(ts,maxiter);
1103:   TSSetMaxTime(ts,tfinal);
1104:   TSSetTimeStep(ts,dt);
1105:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1106:   /* Read command line options for time integration                       */
1107:   TSSetFromOptions(ts);
1108:   /* Set solution vector                                                  */
1109:   TSSetSolution(ts,Y);
1110:   /* Specify left/right-hand side functions                               */
1111:   TSGetType(ts,&time_scheme);

1113:   if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1114:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1115:     TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1116:     MatCreate(PETSC_COMM_WORLD,&Jac);
1117:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1118:     MatSetFromOptions(Jac);
1119:     MatSetUp(Jac);
1120:     TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1121:   } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1122:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1123:     /* and its Jacobian function                                                 */
1124:     TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1125:     MatCreate(PETSC_COMM_WORLD,&Jac);
1126:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1127:     MatSetFromOptions(Jac);
1128:     MatSetUp(Jac);
1129:     TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1130:   }

1132:   /* Solve */
1133:   TSSolve(ts,Y);
1134:   TSGetTime(ts,&final_time);

1136:   /* Get the estimated error, if available */
1137:   VecDuplicate(Y,&Yerr);
1138:   VecZeroEntries(Yerr);
1139:   TSGetTimeError(ts,0,&Yerr);
1140:   VecNorm(Yerr,NORM_2,&err_norm);
1141:   VecDestroy(&Yerr);
1142:   PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);

1144:   /* Exact solution */
1145:   VecDuplicate(Y,&Yex);
1146:   if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1147:     PetscPrintf(PETSC_COMM_WORLD,"Note: There is a difference between the prescribed final time %g and the actual final time, %g.\n",(double)tfinal,(double)final_time);
1148:   }
1149:   ExactSolution(Yex,&ptype[0],final_time,exact_flag);

1151:   /* Calculate Error */
1152:   VecAYPX(Yex,-1.0,Y);
1153:   VecNorm(Yex,NORM_2,error);
1154:   *error = PetscSqrtReal(((*error)*(*error))/N);

1156:   /* Clean up and finalize */
1157:   MatDestroy(&Jac);
1158:   TSDestroy(&ts);
1159:   VecDestroy(&Yex);
1160:   VecDestroy(&Y);

1162:   return 0;
1163: }

1165: int main(int argc, char **argv)
1166: {
1167:   PetscErrorCode  ierr;                       /* Error code                                           */
1168:   char            ptype[256] = "hull1972a1";  /* Problem specification                                */
1169:   PetscInt        n_refine   = 1;             /* Number of refinement levels for convergence analysis */
1170:   PetscReal       refine_fac = 2.0;           /* Refinement factor for dt                             */
1171:   PetscReal       dt_initial = 0.01;          /* Initial default value of dt                          */
1172:   PetscReal       dt;
1173:   PetscReal       tfinal     = 20.0;          /* Final time for the time-integration                  */
1174:   PetscInt        maxiter    = 100000;        /* Maximum number of time-integration iterations        */
1175:   PetscReal       *error;                     /* Array to store the errors for convergence analysis   */
1176:   PetscMPIInt     size;                      /* No of processors                                     */
1177:   PetscBool       flag;                       /* Flag denoting availability of exact solution         */
1178:   PetscInt        r;

1180:   /* Initialize program */
1181:   PetscInitialize(&argc,&argv,(char*)0,help);

1183:   /* Check if running with only 1 proc */
1184:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

1187:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);
1188:   PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1189:   PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1190:   PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1191:   PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1192:   PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1193:   PetscOptionsEnd();

1195:   PetscMalloc1(n_refine,&error);
1196:   for (r = 0,dt = dt_initial; r < n_refine; r++) {
1197:     error[r] = 0;
1198:     if (r > 0) dt /= refine_fac;

1200:     PetscPrintf(PETSC_COMM_WORLD,"Solving ODE \"%s\" with dt %f, final time %f and system size %D.\n",ptype,(double)dt,(double)tfinal,GetSize(&ptype[0]));
1201:     SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1202:     if (flag) {
1203:       /* If exact solution available for the specified ODE */
1204:       if (r > 0) {
1205:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1206:         PetscPrintf(PETSC_COMM_WORLD,"Error           = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1207:       } else {
1208:         PetscPrintf(PETSC_COMM_WORLD,"Error           = %E.\n",error[r]);
1209:       }
1210:     }
1211:   }
1212:   PetscFree(error);
1213:   PetscFinalize();
1214:   return 0;
1215: }

1217: /*TEST

1219:     test:
1220:       suffix: 2
1221:       args: -ts_type glee -final_time 5 -ts_adapt_type none
1222:       timeoutfactor: 3
1223:       requires: !single

1225:     test:
1226:       suffix: 3
1227:       args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor  -ts_max_steps 50  -problem hull1972a3 -ts_adapt_glee_use_local 1
1228:       timeoutfactor: 3
1229:       requires: !single

1231:     test:
1232:       suffix: 4
1233:       args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor  -ts_max_steps 50  -problem hull1972a3  -ts_max_reject 100 -ts_adapt_glee_use_local 0
1234:       timeoutfactor: 3
1235:       requires: !single !__float128

1237: TEST*/