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: PetscInt GetSize(const char *p)
 53: {
 54:   if      ((!strcmp(p,"hull1972a1"))
 55:          ||(!strcmp(p,"hull1972a2"))
 56:          ||(!strcmp(p,"hull1972a3"))
 57:          ||(!strcmp(p,"hull1972a4"))
 58:          ||(!strcmp(p,"hull1972a5"))) return 1;
 59:   else if  (!strcmp(p,"hull1972b1")) return 2;
 60:   else if ((!strcmp(p,"hull1972b2"))
 61:          ||(!strcmp(p,"hull1972b3"))
 62:          ||(!strcmp(p,"hull1972b4"))
 63:          ||(!strcmp(p,"hull1972b5"))) return 3;
 64:   else if ((!strcmp(p,"kulik2013i"))) return 4;
 65:   else if ((!strcmp(p,"hull1972c1"))
 66:          ||(!strcmp(p,"hull1972c2"))
 67:          ||(!strcmp(p,"hull1972c3"))) return 10;
 68:   else if  (!strcmp(p,"hull1972c4")) return 51;
 69:   else PetscFunctionReturn(-1);
 70: }

 72: /****************************************************************/

 74: /* Problem specific functions */

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

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

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

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

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

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

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

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

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

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

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

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

149: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
150: {
151:   const PetscScalar *y;
152:   PetscInt          row = 0,col = 0;
153:   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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

923: /***************************************************************************/

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

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

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

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

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

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

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

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

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

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

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

1143:   /* Exact solution */
1144:   VecDuplicate(Y,&Yex);
1145:   if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1146:     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);
1147:   }
1148:   ExactSolution(Yex,&ptype[0],final_time,exact_flag);

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

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

1161:   return 0;
1162: }

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

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

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

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

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

1199:     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]));
1200:     SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1201:     if (flag) {
1202:       /* If exact solution available for the specified ODE */
1203:       if (r > 0) {
1204:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1205:         PetscPrintf(PETSC_COMM_WORLD,"Error           = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1206:       } else {
1207:         PetscPrintf(PETSC_COMM_WORLD,"Error           = %E.\n",error[r]);
1208:       }
1209:     }
1210:   }
1211:   PetscFree(error);
1212:   PetscFinalize();
1213:   return 0;
1214: }

1216: /*TEST

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

1224:     test:
1225:       suffix: 3
1226:       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
1227:       timeoutfactor: 3
1228:       requires: !single

1230:     test:
1231:       suffix: 4
1232:       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
1233:       timeoutfactor: 3
1234:       requires: !single !__float128

1236: TEST*/