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

 14: */

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

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

 39:     Kulikov2013I -> "kulik2013i"

 41: */

 43: #include <petscts.h>

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

 51: /* Returns the size of the system of equations depending on problem specification */
 52: PetscErrorCode GetSize(const char *p, PetscInt *sz)
 53: {

 56:   if (!strcmp(p, "hull1972a1") || !strcmp(p, "hull1972a2") || !strcmp(p, "hull1972a3") || !strcmp(p, "hull1972a4") || !strcmp(p, "hull1972a5")) *sz = 1;
 57:   else if (!strcmp(p, "hull1972b1")) *sz = 2;
 58:   else if (!strcmp(p, "hull1972b2") || !strcmp(p, "hull1972b3") || !strcmp(p, "hull1972b4") || !strcmp(p, "hull1972b5")) *sz = 3;
 59:   else if (!strcmp(p, "kulik2013i")) *sz = 4;
 60:   else if (!strcmp(p, "hull1972c1") || !strcmp(p, "hull1972c2") || !strcmp(p, "hull1972c3")) *sz = 10;
 61:   else if (!strcmp(p, "hull1972c4")) *sz = 52;
 62:   else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Incorrect problem name %s", p);
 63:   return 0;
 64: }

 66: /****************************************************************/

 68: /* Problem specific functions */

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

 72: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
 73: {
 74:   PetscScalar       *f;
 75:   const PetscScalar *y;

 78:   VecGetArrayRead(Y, &y);
 79:   VecGetArray(F, &f);
 80:   f[0] = -y[0];
 81:   VecRestoreArrayRead(Y, &y);
 82:   VecRestoreArray(F, &f);
 83:   return 0;
 84: }

 86: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
 87: {
 88:   const PetscScalar *y;
 89:   PetscInt           row = 0, col = 0;
 90:   PetscScalar        value = -1.0;

 93:   VecGetArrayRead(Y, &y);
 94:   MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
 95:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 96:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 97:   VecRestoreArrayRead(Y, &y);
 98:   return 0;
 99: }

101: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
102: {
103:   const PetscScalar *y;
104:   PetscScalar       *f;

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

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

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

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

134: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
135: {
136:   const PetscScalar *y;
137:   PetscScalar       *f;

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

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

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

164: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
165: {
166:   PetscScalar       *f;
167:   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;

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

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

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

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

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

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

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

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

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

250:   VecGetArrayRead(Y, &y);
251:   MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
252:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
253:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
254:   VecRestoreArrayRead(Y, &y);
255:   return 0;
256: }

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

260: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
261: {
262:   PetscScalar       *f;
263:   const PetscScalar *y;

266:   VecGetArrayRead(Y, &y);
267:   VecGetArray(F, &f);
268:   f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
269:   VecRestoreArrayRead(Y, &y);
270:   VecRestoreArray(F, &f);
271:   return 0;
272: }

274: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
275: {
276:   const PetscScalar *y;
277:   PetscInt           row = 0, col = 0;
278:   PetscScalar        value;

281:   VecGetArrayRead(Y, &y);
282:   value = 0.25 * (1.0 - 0.05 * y[0]) - (0.25 * y[0]) * 0.05;
283:   MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
284:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
285:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
286:   VecRestoreArrayRead(Y, &y);
287:   return 0;
288: }

290: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
291: {
292:   PetscScalar       *f;
293:   const PetscScalar *y;

296:   VecGetArrayRead(Y, &y);
297:   VecGetArray(F, &f);
298:   f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
299:   VecRestoreArrayRead(Y, &y);
300:   VecRestoreArray(F, &f);
301:   /* Left hand side = ydot - f(y) */
302:   VecAYPX(F, -1.0, Ydot);
303:   return 0;
304: }

306: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
307: {
308:   const PetscScalar *y;
309:   PetscInt           row = 0, col = 0;
310:   PetscScalar        value;

313:   VecGetArrayRead(Y, &y);
314:   value = a - 0.25 * (1.0 - 0.05 * y[0]) + (0.25 * y[0]) * 0.05;
315:   MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
316:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
317:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
318:   VecRestoreArrayRead(Y, &y);
319:   return 0;
320: }

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

324: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
325: {
326:   PetscScalar       *f;
327:   const PetscScalar *y;

330:   VecGetArrayRead(Y, &y);
331:   VecGetArray(F, &f);
332:   f[0] = (y[0] - t) / (y[0] + t);
333:   VecRestoreArrayRead(Y, &y);
334:   VecRestoreArray(F, &f);
335:   return 0;
336: }

338: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
339: {
340:   const PetscScalar *y;
341:   PetscInt           row = 0, col = 0;
342:   PetscScalar        value;

345:   VecGetArrayRead(Y, &y);
346:   value = 2 * t / ((t + y[0]) * (t + y[0]));
347:   MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
348:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
349:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
350:   VecRestoreArrayRead(Y, &y);
351:   return 0;
352: }

354: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
355: {
356:   PetscScalar       *f;
357:   const PetscScalar *y;

360:   VecGetArrayRead(Y, &y);
361:   VecGetArray(F, &f);
362:   f[0] = (y[0] - t) / (y[0] + t);
363:   VecRestoreArrayRead(Y, &y);
364:   VecRestoreArray(F, &f);
365:   /* Left hand side = ydot - f(y) */
366:   VecAYPX(F, -1.0, Ydot);
367:   return 0;
368: }

370: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
371: {
372:   const PetscScalar *y;
373:   PetscInt           row = 0, col = 0;
374:   PetscScalar        value;

377:   VecGetArrayRead(Y, &y);
378:   value = a - 2 * t / ((t + y[0]) * (t + y[0]));
379:   MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
380:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
381:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
382:   VecRestoreArrayRead(Y, &y);
383:   return 0;
384: }

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

388: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
389: {
390:   PetscScalar       *f;
391:   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:   return 0;
401: }

403: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
404: {
405:   PetscScalar       *f;
406:   const PetscScalar *y;

409:   VecGetArrayRead(Y, &y);
410:   VecGetArray(F, &f);
411:   f[0] = 2.0 * (y[0] - y[0] * y[1]);
412:   f[1] = -(y[1] - y[0] * y[1]);
413:   VecRestoreArrayRead(Y, &y);
414:   VecRestoreArray(F, &f);
415:   /* Left hand side = ydot - f(y) */
416:   VecAYPX(F, -1.0, Ydot);
417:   return 0;
418: }

420: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
421: {
422:   const PetscScalar *y;
423:   PetscInt           row[2] = {0, 1};
424:   PetscScalar        value[2][2];

427:   VecGetArrayRead(Y, &y);
428:   value[0][0] = a - 2.0 * (1.0 - y[1]);
429:   value[0][1] = 2.0 * y[0];
430:   value[1][0] = -y[1];
431:   value[1][1] = a + 1.0 - y[0];
432:   MatSetValues(A, 2, &row[0], 2, &row[0], &value[0][0], INSERT_VALUES);
433:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
434:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
435:   VecRestoreArrayRead(Y, &y);
436:   return 0;
437: }

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

441: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
442: {
443:   PetscScalar       *f;
444:   const PetscScalar *y;

447:   VecGetArrayRead(Y, &y);
448:   VecGetArray(F, &f);
449:   f[0] = -y[0] + y[1];
450:   f[1] = y[0] - 2.0 * y[1] + y[2];
451:   f[2] = y[1] - y[2];
452:   VecRestoreArrayRead(Y, &y);
453:   VecRestoreArray(F, &f);
454:   return 0;
455: }

457: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
458: {
459:   PetscScalar       *f;
460:   const PetscScalar *y;

463:   VecGetArrayRead(Y, &y);
464:   VecGetArray(F, &f);
465:   f[0] = -y[0] + y[1];
466:   f[1] = y[0] - 2.0 * y[1] + y[2];
467:   f[2] = y[1] - y[2];
468:   VecRestoreArrayRead(Y, &y);
469:   VecRestoreArray(F, &f);
470:   /* Left hand side = ydot - f(y) */
471:   VecAYPX(F, -1.0, Ydot);
472:   return 0;
473: }

475: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
476: {
477:   const PetscScalar *y;
478:   PetscInt           row[3] = {0, 1, 2};
479:   PetscScalar        value[3][3];

482:   VecGetArrayRead(Y, &y);
483:   value[0][0] = a + 1.0;
484:   value[0][1] = -1.0;
485:   value[0][2] = 0;
486:   value[1][0] = -1.0;
487:   value[1][1] = a + 2.0;
488:   value[1][2] = -1.0;
489:   value[2][0] = 0;
490:   value[2][1] = -1.0;
491:   value[2][2] = a + 1.0;
492:   MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
493:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
494:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
495:   VecRestoreArrayRead(Y, &y);
496:   return 0;
497: }

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

501: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
502: {
503:   PetscScalar       *f;
504:   const PetscScalar *y;

507:   VecGetArrayRead(Y, &y);
508:   VecGetArray(F, &f);
509:   f[0] = -y[0];
510:   f[1] = y[0] - y[1] * y[1];
511:   f[2] = y[1] * y[1];
512:   VecRestoreArrayRead(Y, &y);
513:   VecRestoreArray(F, &f);
514:   return 0;
515: }

517: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
518: {
519:   PetscScalar       *f;
520:   const PetscScalar *y;

523:   VecGetArrayRead(Y, &y);
524:   VecGetArray(F, &f);
525:   f[0] = -y[0];
526:   f[1] = y[0] - y[1] * y[1];
527:   f[2] = y[1] * y[1];
528:   VecRestoreArrayRead(Y, &y);
529:   VecRestoreArray(F, &f);
530:   /* Left hand side = ydot - f(y) */
531:   VecAYPX(F, -1.0, Ydot);
532:   return 0;
533: }

535: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
536: {
537:   const PetscScalar *y;
538:   PetscInt           row[3] = {0, 1, 2};
539:   PetscScalar        value[3][3];

542:   VecGetArrayRead(Y, &y);
543:   value[0][0] = a + 1.0;
544:   value[0][1] = 0;
545:   value[0][2] = 0;
546:   value[1][0] = -1.0;
547:   value[1][1] = a + 2.0 * y[1];
548:   value[1][2] = 0;
549:   value[2][0] = 0;
550:   value[2][1] = -2.0 * y[1];
551:   value[2][2] = a;
552:   MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
553:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
554:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
555:   VecRestoreArrayRead(Y, &y);
556:   return 0;
557: }

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

561: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
562: {
563:   PetscScalar       *f;
564:   const PetscScalar *y;

567:   VecGetArrayRead(Y, &y);
568:   VecGetArray(F, &f);
569:   f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
570:   f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
571:   f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
572:   VecRestoreArrayRead(Y, &y);
573:   VecRestoreArray(F, &f);
574:   return 0;
575: }

577: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
578: {
579:   PetscScalar       *f;
580:   const PetscScalar *y;

583:   VecGetArrayRead(Y, &y);
584:   VecGetArray(F, &f);
585:   f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
586:   f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
587:   f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
588:   VecRestoreArrayRead(Y, &y);
589:   VecRestoreArray(F, &f);
590:   /* Left hand side = ydot - f(y) */
591:   VecAYPX(F, -1.0, Ydot);
592:   return 0;
593: }

595: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
596: {
597:   const PetscScalar *y;
598:   PetscInt           row[3] = {0, 1, 2};
599:   PetscScalar        value[3][3], fac, fac2;

602:   VecGetArrayRead(Y, &y);
603:   fac         = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -1.5);
604:   fac2        = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -0.5);
605:   value[0][0] = a + (y[1] * y[1] * y[2]) * fac;
606:   value[0][1] = 1.0 - (y[0] * y[1] * y[2]) * fac;
607:   value[0][2] = y[0] * fac2;
608:   value[1][0] = -1.0 - y[0] * y[1] * y[2] * fac;
609:   value[1][1] = a + y[0] * y[0] * y[2] * fac;
610:   value[1][2] = y[1] * fac2;
611:   value[2][0] = -y[1] * y[1] * fac;
612:   value[2][1] = y[0] * y[1] * fac;
613:   value[2][2] = a;
614:   MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
615:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
616:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
617:   VecRestoreArrayRead(Y, &y);
618:   return 0;
619: }

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

623: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
624: {
625:   PetscScalar       *f;
626:   const PetscScalar *y;

629:   VecGetArrayRead(Y, &y);
630:   VecGetArray(F, &f);
631:   f[0] = y[1] * y[2];
632:   f[1] = -y[0] * y[2];
633:   f[2] = -0.51 * y[0] * y[1];
634:   VecRestoreArrayRead(Y, &y);
635:   VecRestoreArray(F, &f);
636:   return 0;
637: }

639: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
640: {
641:   PetscScalar       *f;
642:   const PetscScalar *y;

645:   VecGetArrayRead(Y, &y);
646:   VecGetArray(F, &f);
647:   f[0] = y[1] * y[2];
648:   f[1] = -y[0] * y[2];
649:   f[2] = -0.51 * y[0] * y[1];
650:   VecRestoreArrayRead(Y, &y);
651:   VecRestoreArray(F, &f);
652:   /* Left hand side = ydot - f(y) */
653:   VecAYPX(F, -1.0, Ydot);
654:   return 0;
655: }

657: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
658: {
659:   const PetscScalar *y;
660:   PetscInt           row[3] = {0, 1, 2};
661:   PetscScalar        value[3][3];

664:   VecGetArrayRead(Y, &y);
665:   value[0][0] = a;
666:   value[0][1] = -y[2];
667:   value[0][2] = -y[1];
668:   value[1][0] = y[2];
669:   value[1][1] = a;
670:   value[1][2] = y[0];
671:   value[2][0] = 0.51 * y[1];
672:   value[2][1] = 0.51 * y[0];
673:   value[2][2] = a;
674:   MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
675:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
676:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
677:   VecRestoreArrayRead(Y, &y);
678:   return 0;
679: }

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

683: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
684: {
685:   PetscScalar       *f;
686:   const PetscScalar *y;

689:   VecGetArrayRead(Y, &y);
690:   VecGetArray(F, &f);
691:   f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
692:   f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
693:   f[2] = 2. * t * y[3];
694:   f[3] = -2. * t * PetscLogScalar(y[0]);
695:   VecRestoreArrayRead(Y, &y);
696:   VecRestoreArray(F, &f);
697:   return 0;
698: }

700: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
701: {
702:   const PetscScalar *y;
703:   PetscInt           row[4] = {0, 1, 2, 3};
704:   PetscScalar        value[4][4];
705:   PetscScalar        m1, m2;
707:   VecGetArrayRead(Y, &y);
708:   m1          = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.));
709:   m2          = 2. * t * PetscPowScalar(y[1], 1. / 5.);
710:   value[0][0] = 0.;
711:   value[0][1] = m1;
712:   value[0][2] = 0.;
713:   value[0][3] = m2;
714:   m1          = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
715:   m2          = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.));
716:   value[1][0] = 0.;
717:   value[1][1] = 0.;
718:   value[1][2] = m1;
719:   value[1][3] = m2;
720:   value[2][0] = 0.;
721:   value[2][1] = 0.;
722:   value[2][2] = 0.;
723:   value[2][3] = 2 * t;
724:   value[3][0] = -2. * t / y[0];
725:   value[3][1] = 0.;
726:   value[3][2] = 0.;
727:   value[3][3] = 0.;
728:   MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES);
729:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
730:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
731:   VecRestoreArrayRead(Y, &y);
732:   return 0;
733: }

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

741:   VecGetArrayRead(Y, &y);
742:   VecGetArray(F, &f);
743:   f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
744:   f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
745:   f[2] = 2. * t * y[3];
746:   f[3] = -2. * t * PetscLogScalar(y[0]);
747:   VecRestoreArrayRead(Y, &y);
748:   VecRestoreArray(F, &f);
749:   /* Left hand side = ydot - f(y) */
750:   VecAYPX(F, -1.0, Ydot);
751:   return 0;
752: }

754: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
755: {
756:   const PetscScalar *y;
757:   PetscInt           row[4] = {0, 1, 2, 3};
758:   PetscScalar        value[4][4];
759:   PetscScalar        m1, m2;

762:   VecGetArrayRead(Y, &y);
763:   m1          = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.));
764:   m2          = 2. * t * PetscPowScalar(y[1], 1. / 5.);
765:   value[0][0] = a;
766:   value[0][1] = m1;
767:   value[0][2] = 0.;
768:   value[0][3] = m2;
769:   m1          = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
770:   m2          = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.));
771:   value[1][0] = 0.;
772:   value[1][1] = a;
773:   value[1][2] = m1;
774:   value[1][3] = m2;
775:   value[2][0] = 0.;
776:   value[2][1] = 0.;
777:   value[2][2] = a;
778:   value[2][3] = 2 * t;
779:   value[3][0] = -2. * t / y[0];
780:   value[3][1] = 0.;
781:   value[3][2] = 0.;
782:   value[3][3] = a;
783:   MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES);
784:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
785:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
786:   VecRestoreArrayRead(Y, &y);
787:   return 0;
788: }

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

792: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
793: {
794:   PetscScalar       *f;
795:   const PetscScalar *y;
796:   PetscInt           N, i;

799:   VecGetSize(Y, &N);
800:   VecGetArrayRead(Y, &y);
801:   VecGetArray(F, &f);
802:   f[0] = -y[0];
803:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - y[i];
804:   f[N - 1] = y[N - 2];
805:   VecRestoreArrayRead(Y, &y);
806:   VecRestoreArray(F, &f);
807:   return 0;
808: }

810: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
811: {
812:   PetscScalar       *f;
813:   const PetscScalar *y;
814:   PetscInt           N, i;

817:   VecGetSize(Y, &N);
818:   VecGetArrayRead(Y, &y);
819:   VecGetArray(F, &f);
820:   f[0] = -y[0];
821:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - y[i];
822:   f[N - 1] = y[N - 2];
823:   VecRestoreArrayRead(Y, &y);
824:   VecRestoreArray(F, &f);
825:   /* Left hand side = ydot - f(y) */
826:   VecAYPX(F, -1.0, Ydot);
827:   return 0;
828: }

830: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
831: {
832:   const PetscScalar *y;
833:   PetscInt           N, i, col[2];
834:   PetscScalar        value[2];

837:   VecGetSize(Y, &N);
838:   VecGetArrayRead(Y, &y);
839:   i        = 0;
840:   value[0] = a + 1;
841:   col[0]   = 0;
842:   value[1] = 0;
843:   col[1]   = 1;
844:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
845:   for (i = 0; i < N; i++) {
846:     value[0] = -1;
847:     col[0]   = i - 1;
848:     value[1] = a + 1;
849:     col[1]   = i;
850:     MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
851:   }
852:   i        = N - 1;
853:   value[0] = -1;
854:   col[0]   = N - 2;
855:   value[1] = a;
856:   col[1]   = N - 1;
857:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
858:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
859:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
860:   VecRestoreArrayRead(Y, &y);
861:   return 0;
862: }

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

866: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
867: {
868:   const PetscScalar *y;
869:   PetscScalar       *f;
870:   PetscInt           N, i;

873:   VecGetSize(Y, &N);
874:   VecGetArrayRead(Y, &y);
875:   VecGetArray(F, &f);
876:   f[0] = -y[0];
877:   for (i = 1; i < N - 1; i++) f[i] = (PetscReal)i * y[i - 1] - (PetscReal)(i + 1) * y[i];
878:   f[N - 1] = (PetscReal)(N - 1) * y[N - 2];
879:   VecRestoreArrayRead(Y, &y);
880:   VecRestoreArray(F, &f);
881:   return 0;
882: }

884: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
885: {
886:   PetscScalar       *f;
887:   const PetscScalar *y;
888:   PetscInt           N, i;

891:   VecGetSize(Y, &N);
892:   VecGetArrayRead(Y, &y);
893:   VecGetArray(F, &f);
894:   f[0] = -y[0];
895:   for (i = 1; i < N - 1; i++) f[i] = (PetscReal)i * y[i - 1] - (PetscReal)(i + 1) * y[i];
896:   f[N - 1] = (PetscReal)(N - 1) * y[N - 2];
897:   VecRestoreArrayRead(Y, &y);
898:   VecRestoreArray(F, &f);
899:   /* Left hand side = ydot - f(y) */
900:   VecAYPX(F, -1.0, Ydot);
901:   return 0;
902: }

904: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
905: {
906:   const PetscScalar *y;
907:   PetscInt           N, i, col[2];
908:   PetscScalar        value[2];

911:   VecGetSize(Y, &N);
912:   VecGetArrayRead(Y, &y);
913:   i        = 0;
914:   value[0] = a + 1;
915:   col[0]   = 0;
916:   value[1] = 0;
917:   col[1]   = 1;
918:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
919:   for (i = 0; i < N; i++) {
920:     value[0] = -(PetscReal)i;
921:     col[0]   = i - 1;
922:     value[1] = a + (PetscReal)(i + 1);
923:     col[1]   = i;
924:     MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
925:   }
926:   i        = N - 1;
927:   value[0] = -(PetscReal)(N - 1);
928:   col[0]   = N - 2;
929:   value[1] = a;
930:   col[1]   = N - 1;
931:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
932:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
933:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
934:   VecRestoreArrayRead(Y, &y);
935:   return 0;
936: }

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

940: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
941: {
942:   PetscScalar       *f;
943:   const PetscScalar *y;
944:   PetscInt           N, i;

947:   VecGetSize(Y, &N);
948:   VecGetArrayRead(Y, &y);
949:   VecGetArray(F, &f);
950:   f[0] = -2.0 * y[0] + y[1];
951:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - 2.0 * y[i] + y[i + 1];
952:   f[N - 1] = y[N - 2] - 2.0 * y[N - 1];
953:   VecRestoreArrayRead(Y, &y);
954:   VecRestoreArray(F, &f);
955:   return 0;
956: }

958: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
959: {
960:   PetscScalar       *f;
961:   const PetscScalar *y;
962:   PetscInt           N, i;

965:   VecGetSize(Y, &N);
966:   VecGetArrayRead(Y, &y);
967:   VecGetArray(F, &f);
968:   f[0] = -2.0 * y[0] + y[1];
969:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - 2.0 * y[i] + y[i + 1];
970:   f[N - 1] = y[N - 2] - 2.0 * y[N - 1];
971:   VecRestoreArrayRead(Y, &y);
972:   VecRestoreArray(F, &f);
973:   /* Left hand side = ydot - f(y) */
974:   VecAYPX(F, -1.0, Ydot);
975:   return 0;
976: }

978: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
979: {
980:   const PetscScalar *y;
981:   PetscScalar        value[3];
982:   PetscInt           N, i, col[3];

985:   VecGetSize(Y, &N);
986:   VecGetArrayRead(Y, &y);
987:   for (i = 0; i < N; i++) {
988:     if (i == 0) {
989:       value[0] = a + 2;
990:       col[0]   = i;
991:       value[1] = -1;
992:       col[1]   = i + 1;
993:       value[2] = 0;
994:       col[2]   = i + 2;
995:     } else if (i == N - 1) {
996:       value[0] = 0;
997:       col[0]   = i - 2;
998:       value[1] = -1;
999:       col[1]   = i - 1;
1000:       value[2] = a + 2;
1001:       col[2]   = i;
1002:     } else {
1003:       value[0] = -1;
1004:       col[0]   = i - 1;
1005:       value[1] = a + 2;
1006:       col[1]   = i;
1007:       value[2] = -1;
1008:       col[2]   = i + 1;
1009:     }
1010:     MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
1011:   }
1012:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1013:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1014:   VecRestoreArrayRead(Y, &y);
1015:   return 0;
1016: }

1018: /***************************************************************************/

1020: /* Sets the initial solution for the IVP and sets up the function pointers*/
1021: PetscErrorCode Initialize(Vec Y, void *s)
1022: {
1023:   char        *p = (char *)s;
1024:   PetscScalar *y;
1025:   PetscReal    t0;
1026:   PetscInt     N;
1027:   PetscBool    flg;

1030:   GetSize((const char *)s, &N);
1031:   VecZeroEntries(Y);
1032:   VecGetArray(Y, &y);
1033:   if (!strcmp(p, "hull1972a1")) {
1034:     y[0]        = 1.0;
1035:     RHSFunction = RHSFunction_Hull1972A1;
1036:     RHSJacobian = RHSJacobian_Hull1972A1;
1037:     IFunction   = IFunction_Hull1972A1;
1038:     IJacobian   = IJacobian_Hull1972A1;
1039:   } else if (!strcmp(p, "hull1972a2")) {
1040:     y[0]        = 1.0;
1041:     RHSFunction = RHSFunction_Hull1972A2;
1042:     RHSJacobian = RHSJacobian_Hull1972A2;
1043:     IFunction   = IFunction_Hull1972A2;
1044:     IJacobian   = IJacobian_Hull1972A2;
1045:   } else if (!strcmp(p, "hull1972a3")) {
1046:     y[0]        = 1.0;
1047:     RHSFunction = RHSFunction_Hull1972A3;
1048:     RHSJacobian = RHSJacobian_Hull1972A3;
1049:     IFunction   = IFunction_Hull1972A3;
1050:     IJacobian   = IJacobian_Hull1972A3;
1051:   } else if (!strcmp(p, "hull1972a4")) {
1052:     y[0]        = 1.0;
1053:     RHSFunction = RHSFunction_Hull1972A4;
1054:     RHSJacobian = RHSJacobian_Hull1972A4;
1055:     IFunction   = IFunction_Hull1972A4;
1056:     IJacobian   = IJacobian_Hull1972A4;
1057:   } else if (!strcmp(p, "hull1972a5")) {
1058:     y[0]        = 4.0;
1059:     RHSFunction = RHSFunction_Hull1972A5;
1060:     RHSJacobian = RHSJacobian_Hull1972A5;
1061:     IFunction   = IFunction_Hull1972A5;
1062:     IJacobian   = IJacobian_Hull1972A5;
1063:   } else if (!strcmp(p, "hull1972b1")) {
1064:     y[0]        = 1.0;
1065:     y[1]        = 3.0;
1066:     RHSFunction = RHSFunction_Hull1972B1;
1067:     IFunction   = IFunction_Hull1972B1;
1068:     IJacobian   = IJacobian_Hull1972B1;
1069:   } else if (!strcmp(p, "hull1972b2")) {
1070:     y[0]        = 2.0;
1071:     y[1]        = 0.0;
1072:     y[2]        = 1.0;
1073:     RHSFunction = RHSFunction_Hull1972B2;
1074:     IFunction   = IFunction_Hull1972B2;
1075:     IJacobian   = IJacobian_Hull1972B2;
1076:   } else if (!strcmp(p, "hull1972b3")) {
1077:     y[0]        = 1.0;
1078:     y[1]        = 0.0;
1079:     y[2]        = 0.0;
1080:     RHSFunction = RHSFunction_Hull1972B3;
1081:     IFunction   = IFunction_Hull1972B3;
1082:     IJacobian   = IJacobian_Hull1972B3;
1083:   } else if (!strcmp(p, "hull1972b4")) {
1084:     y[0]        = 3.0;
1085:     y[1]        = 0.0;
1086:     y[2]        = 0.0;
1087:     RHSFunction = RHSFunction_Hull1972B4;
1088:     IFunction   = IFunction_Hull1972B4;
1089:     IJacobian   = IJacobian_Hull1972B4;
1090:   } else if (!strcmp(p, "hull1972b5")) {
1091:     y[0]        = 0.0;
1092:     y[1]        = 1.0;
1093:     y[2]        = 1.0;
1094:     RHSFunction = RHSFunction_Hull1972B5;
1095:     IFunction   = IFunction_Hull1972B5;
1096:     IJacobian   = IJacobian_Hull1972B5;
1097:   } else if (!strcmp(p, "kulik2013i")) {
1098:     t0          = 0.;
1099:     y[0]        = PetscExpReal(PetscSinReal(t0 * t0));
1100:     y[1]        = PetscExpReal(5. * PetscSinReal(t0 * t0));
1101:     y[2]        = PetscSinReal(t0 * t0) + 1.0;
1102:     y[3]        = PetscCosReal(t0 * t0);
1103:     RHSFunction = RHSFunction_Kulikov2013I;
1104:     RHSJacobian = RHSJacobian_Kulikov2013I;
1105:     IFunction   = IFunction_Kulikov2013I;
1106:     IJacobian   = IJacobian_Kulikov2013I;
1107:   } else if (!strcmp(p, "hull1972c1")) {
1108:     y[0]        = 1.0;
1109:     RHSFunction = RHSFunction_Hull1972C1;
1110:     IFunction   = IFunction_Hull1972C1;
1111:     IJacobian   = IJacobian_Hull1972C1;
1112:   } else if (!strcmp(p, "hull1972c2")) {
1113:     y[0]        = 1.0;
1114:     RHSFunction = RHSFunction_Hull1972C2;
1115:     IFunction   = IFunction_Hull1972C2;
1116:     IJacobian   = IJacobian_Hull1972C2;
1117:   } else if (!strcmp(p, "hull1972c3") || !strcmp(p, "hull1972c4")) {
1118:     y[0]        = 1.0;
1119:     RHSFunction = RHSFunction_Hull1972C34;
1120:     IFunction   = IFunction_Hull1972C34;
1121:     IJacobian   = IJacobian_Hull1972C34;
1122:   }
1123:   PetscOptionsGetScalarArray(NULL, NULL, "-yinit", y, &N, &flg);
1124:   VecRestoreArray(Y, &y);
1125:   return 0;
1126: }

1128: /* Calculates the exact solution to problems that have one */
1129: PetscErrorCode ExactSolution(Vec Y, void *s, PetscReal t, PetscBool *flag)
1130: {
1131:   char        *p = (char *)s;
1132:   PetscScalar *y;

1135:   if (!strcmp(p, "hull1972a1")) {
1136:     VecGetArray(Y, &y);
1137:     y[0]  = PetscExpReal(-t);
1138:     *flag = PETSC_TRUE;
1139:     VecRestoreArray(Y, &y);
1140:   } else if (!strcmp(p, "hull1972a2")) {
1141:     VecGetArray(Y, &y);
1142:     y[0]  = 1.0 / PetscSqrtReal(t + 1);
1143:     *flag = PETSC_TRUE;
1144:     VecRestoreArray(Y, &y);
1145:   } else if (!strcmp(p, "hull1972a3")) {
1146:     VecGetArray(Y, &y);
1147:     y[0]  = PetscExpReal(PetscSinReal(t));
1148:     *flag = PETSC_TRUE;
1149:     VecRestoreArray(Y, &y);
1150:   } else if (!strcmp(p, "hull1972a4")) {
1151:     VecGetArray(Y, &y);
1152:     y[0]  = 20.0 / (1 + 19.0 * PetscExpReal(-t / 4.0));
1153:     *flag = PETSC_TRUE;
1154:     VecRestoreArray(Y, &y);
1155:   } else if (!strcmp(p, "kulik2013i")) {
1156:     VecGetArray(Y, &y);
1157:     y[0]  = PetscExpReal(PetscSinReal(t * t));
1158:     y[1]  = PetscExpReal(5. * PetscSinReal(t * t));
1159:     y[2]  = PetscSinReal(t * t) + 1.0;
1160:     y[3]  = PetscCosReal(t * t);
1161:     *flag = PETSC_TRUE;
1162:     VecRestoreArray(Y, &y);
1163:   } else {
1164:     VecSet(Y, 0);
1165:     *flag = PETSC_FALSE;
1166:   }
1167:   return 0;
1168: }

1170: /* Solves the specified ODE and computes the error if exact solution is available */
1171: PetscErrorCode SolveODE(char *ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1172: {
1173:   TS        ts;          /* time-integrator                        */
1174:   Vec       Y;           /* Solution vector                        */
1175:   Vec       Yex;         /* Exact solution                         */
1176:   PetscInt  N;           /* Size of the system of equations        */
1177:   TSType    time_scheme; /* Type of time-integration scheme        */
1178:   Mat       Jac = NULL;  /* Jacobian matrix                        */
1179:   Vec       Yerr;        /* Auxiliary solution vector              */
1180:   PetscReal err_norm;    /* Estimated error norm                   */
1181:   PetscReal final_time;  /* Actual final time from the integrator  */

1184:   GetSize((const char *)&ptype[0], &N);
1186:   VecCreate(PETSC_COMM_WORLD, &Y);
1187:   VecSetSizes(Y, N, PETSC_DECIDE);
1188:   VecSetUp(Y);
1189:   VecSet(Y, 0);

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

1194:   /* Create and initialize the time-integrator                            */
1195:   TSCreate(PETSC_COMM_WORLD, &ts);
1196:   /* Default time integration options                                     */
1197:   TSSetType(ts, TSRK);
1198:   TSSetMaxSteps(ts, maxiter);
1199:   TSSetMaxTime(ts, tfinal);
1200:   TSSetTimeStep(ts, dt);
1201:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1202:   /* Read command line options for time integration                       */
1203:   TSSetFromOptions(ts);
1204:   /* Set solution vector                                                  */
1205:   TSSetSolution(ts, Y);
1206:   /* Specify left/right-hand side functions                               */
1207:   TSGetType(ts, &time_scheme);

1209:   if ((!strcmp(time_scheme, TSEULER)) || (!strcmp(time_scheme, TSRK)) || (!strcmp(time_scheme, TSSSP) || (!strcmp(time_scheme, TSGLEE)))) {
1210:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1211:     TSSetRHSFunction(ts, NULL, RHSFunction, &ptype[0]);
1212:     MatCreate(PETSC_COMM_WORLD, &Jac);
1213:     MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N);
1214:     MatSetFromOptions(Jac);
1215:     MatSetUp(Jac);
1216:     TSSetRHSJacobian(ts, Jac, Jac, RHSJacobian, &ptype[0]);
1217:   } else if ((!strcmp(time_scheme, TSTHETA)) || (!strcmp(time_scheme, TSBEULER)) || (!strcmp(time_scheme, TSCN)) || (!strcmp(time_scheme, TSALPHA)) || (!strcmp(time_scheme, TSARKIMEX))) {
1218:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1219:     /* and its Jacobian function                                                 */
1220:     TSSetIFunction(ts, NULL, IFunction, &ptype[0]);
1221:     MatCreate(PETSC_COMM_WORLD, &Jac);
1222:     MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N);
1223:     MatSetFromOptions(Jac);
1224:     MatSetUp(Jac);
1225:     TSSetIJacobian(ts, Jac, Jac, IJacobian, &ptype[0]);
1226:   }

1228:   /* Solve */
1229:   TSSolve(ts, Y);
1230:   TSGetTime(ts, &final_time);

1232:   /* Get the estimated error, if available */
1233:   VecDuplicate(Y, &Yerr);
1234:   VecZeroEntries(Yerr);
1235:   TSGetTimeError(ts, 0, &Yerr);
1236:   VecNorm(Yerr, NORM_2, &err_norm);
1237:   VecDestroy(&Yerr);
1238:   PetscPrintf(PETSC_COMM_WORLD, "Estimated Error = %e.\n", (double)err_norm);

1240:   /* Exact solution */
1241:   VecDuplicate(Y, &Yex);
1242:   if (PetscAbsReal(final_time - tfinal) > 2. * PETSC_MACHINE_EPSILON) {
1243:     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);
1244:   }
1245:   ExactSolution(Yex, &ptype[0], final_time, exact_flag);

1247:   /* Calculate Error */
1248:   VecAYPX(Yex, -1.0, Y);
1249:   VecNorm(Yex, NORM_2, error);
1250:   *error = PetscSqrtReal(((*error) * (*error)) / N);

1252:   /* Clean up and finalize */
1253:   MatDestroy(&Jac);
1254:   TSDestroy(&ts);
1255:   VecDestroy(&Yex);
1256:   VecDestroy(&Y);

1258:   return 0;
1259: }

1261: int main(int argc, char **argv)
1262: {
1263:   char        ptype[256] = "hull1972a1"; /* Problem specification                                */
1264:   PetscInt    n_refine   = 1;            /* Number of refinement levels for convergence analysis */
1265:   PetscReal   refine_fac = 2.0;          /* Refinement factor for dt                             */
1266:   PetscReal   dt_initial = 0.01;         /* Initial default value of dt                          */
1267:   PetscReal   dt;
1268:   PetscReal   tfinal  = 20.0;   /* Final time for the time-integration                  */
1269:   PetscInt    maxiter = 100000; /* Maximum number of time-integration iterations        */
1270:   PetscReal  *error;            /* Array to store the errors for convergence analysis   */
1271:   PetscMPIInt size;             /* No of processors                                     */
1272:   PetscBool   flag;             /* Flag denoting availability of exact solution         */
1273:   PetscInt    r, N;

1275:   /* Initialize program */
1277:   GetSize(&ptype[0], &N);
1278:   PetscInitialize(&argc, &argv, (char *)0, help);

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

1284:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL);
1285:   PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL);
1286:   PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL);
1287:   PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL);
1288:   PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL);
1289:   PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL);
1290:   PetscOptionsEnd();

1292:   PetscMalloc1(n_refine, &error);
1293:   for (r = 0, dt = dt_initial; r < n_refine; r++) {
1294:     error[r] = 0;
1295:     if (r > 0) dt /= refine_fac;

1297:     PetscPrintf(PETSC_COMM_WORLD, "Solving ODE \"%s\" with dt %f, final time %f and system size %" PetscInt_FMT ".\n", ptype, (double)dt, (double)tfinal, N);
1298:     SolveODE(&ptype[0], dt, tfinal, maxiter, &error[r], &flag);
1299:     if (flag) {
1300:       /* If exact solution available for the specified ODE */
1301:       if (r > 0) {
1302:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r - 1])) / (-PetscLogReal(refine_fac));
1303:         PetscPrintf(PETSC_COMM_WORLD, "Error           = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate);
1304:       } else {
1305:         PetscPrintf(PETSC_COMM_WORLD, "Error           = %E.\n", (double)error[r]);
1306:       }
1307:     }
1308:   }
1309:   PetscFree(error);
1310:   PetscFinalize();
1311:   return 0;
1312: }

1314: /*TEST

1316:     test:
1317:       suffix: 2
1318:       args: -ts_type glee -final_time 5 -ts_adapt_type none
1319:       timeoutfactor: 3
1320:       requires: !single

1322:     test:
1323:       suffix: 3
1324:       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
1325:       timeoutfactor: 3
1326:       requires: !single

1328:     test:
1329:       suffix: 4
1330:       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
1331:       timeoutfactor: 3
1332:       requires: !single !__float128

1334: TEST*/