Actual source code: ex31.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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: {
 56:   if      ((!strcmp(p,"hull1972a1"))
 57:          ||(!strcmp(p,"hull1972a2"))
 58:          ||(!strcmp(p,"hull1972a3"))
 59:          ||(!strcmp(p,"hull1972a4"))
 60:          ||(!strcmp(p,"hull1972a5")) )  PetscFunctionReturn(1);
 61:   else if  (!strcmp(p,"hull1972b1")  )  PetscFunctionReturn(2);
 62:   else if ((!strcmp(p,"hull1972b2"))
 63:          ||(!strcmp(p,"hull1972b3"))
 64:          ||(!strcmp(p,"hull1972b4"))
 65:          ||(!strcmp(p,"hull1972b5")) )  PetscFunctionReturn(3);
 66:   else if ((!strcmp(p,"kulik2013i")) )  PetscFunctionReturn(4);
 67:   else if ((!strcmp(p,"hull1972c1"))
 68:          ||(!strcmp(p,"hull1972c2"))
 69:          ||(!strcmp(p,"hull1972c3")) )  PetscFunctionReturn(10);
 70:   else if  (!strcmp(p,"hull1972c4")  )  PetscFunctionReturn(51);
 71:   else                                  PetscFunctionReturn(-1);
 72: }

 74: /****************************************************************/

 76: /* Problem specific functions */

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

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

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

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

103:   VecGetArrayRead(Y,&y);
104:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
105:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
106:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
107:   VecRestoreArrayRead(Y,&y);
108:   return(0);
109: }

111: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
112: {
113:   PetscErrorCode    ierr;
114:   const PetscScalar *y;
115:   PetscScalar       *f;

118:   VecGetArrayRead(Y,&y);
119:   VecGetArray(F,&f);
120:   f[0] = -y[0];
121:   VecRestoreArrayRead(Y,&y);
122:   VecRestoreArray(F,&f);
123:   /* Left hand side = ydot - f(y) */
124:   VecAYPX(F,-1.0,Ydot);
125:   return(0);
126: }

128: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
129: {
130:   PetscErrorCode    ierr;
131:   const PetscScalar *y;
132:   PetscInt          row = 0,col = 0;
133:   PetscScalar       value = a - 1.0;

136:   VecGetArrayRead(Y,&y);
137:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
138:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
139:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
140:   VecRestoreArrayRead(Y,&y);
141:   return(0);
142: }

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

146: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
147: {
148:   PetscErrorCode    ierr;
149:   const PetscScalar *y;
150:   PetscScalar       *f;  

153:   VecGetArrayRead(Y,&y);
154:   VecGetArray(F,&f);
155:   f[0] = -0.5*y[0]*y[0]*y[0];
156:   VecRestoreArrayRead(Y,&y);
157:   VecRestoreArray(F,&f);
158:   return(0);
159: }

161: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
162: {
163:   PetscErrorCode    ierr;
164:   const PetscScalar *y;
165:   PetscInt          row = 0,col = 0;
166:   PetscScalar       value;

169:   VecGetArrayRead(Y,&y);
170:   value = -0.5*3.0*y[0]*y[0];
171:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
172:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
173:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
174:   VecRestoreArrayRead(Y,&y);
175:   return(0);
176: }

178: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
179: {
180:   PetscErrorCode    ierr;
181:   PetscScalar       *f;
182:   const PetscScalar *y;

185:   VecGetArrayRead(Y,&y);
186:   VecGetArray(F,&f);
187:   f[0] = -0.5*y[0]*y[0]*y[0];
188:   VecRestoreArrayRead(Y,&y);
189:   VecRestoreArray(F,&f);
190:   /* Left hand side = ydot - f(y) */
191:   VecAYPX(F,-1.0,Ydot);
192:   return(0);
193: }

195: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
196: {
197:   PetscErrorCode    ierr;
198:   const PetscScalar *y;
199:   PetscInt          row = 0,col = 0;
200:   PetscScalar       value;

203:   VecGetArrayRead(Y,&y);
204:   value = a + 0.5*3.0*y[0]*y[0];
205:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
206:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
207:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
208:   VecRestoreArrayRead(Y,&y);
209:   return(0);
210: }

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

214: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
215: {
216:   PetscErrorCode    ierr;
217:   const PetscScalar *y;
218:   PetscScalar       *f;

221:   VecGetArrayRead(Y,&y);
222:   VecGetArray(F,&f);
223:   f[0] = y[0]*PetscCosReal(t);
224:   VecRestoreArrayRead(Y,&y);
225:   VecRestoreArray(F,&f);
226:   return(0);
227: }

229: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
230: {
231:   PetscErrorCode    ierr;
232:   const PetscScalar *y;
233:   PetscInt          row = 0,col = 0;
234:   PetscScalar       value = PetscCosReal(t);

237:   VecGetArrayRead(Y,&y);
238:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
239:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
240:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
241:   VecRestoreArrayRead(Y,&y);
242:   return(0);
243: }

245: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
246: {
247:   PetscErrorCode    ierr;
248:   PetscScalar       *f;
249:   const PetscScalar *y;

252:   VecGetArrayRead(Y,&y);
253:   VecGetArray(F,&f);
254:   f[0] = y[0]*PetscCosReal(t);
255:   VecRestoreArrayRead(Y,&y);
256:   VecRestoreArray(F,&f);
257:   /* Left hand side = ydot - f(y) */
258:   VecAYPX(F,-1.0,Ydot);
259:   return(0);
260: }

262: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
263: {
264:   PetscErrorCode    ierr;
265:   const PetscScalar *y;
266:   PetscInt          row = 0,col = 0;
267:   PetscScalar       value = a - PetscCosReal(t);

270:   VecGetArrayRead(Y,&y);
271:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
272:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
273:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
274:   VecRestoreArrayRead(Y,&y);
275:   return(0);
276: }

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

280: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
281: {
282:   PetscErrorCode    ierr;
283:   PetscScalar       *f;
284:   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:   return(0);
293: }

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

303:   VecGetArrayRead(Y,&y);
304:   value = 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: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
313: {
314:   PetscErrorCode    ierr;
315:   PetscScalar       *f;
316:   const PetscScalar *y;

319:   VecGetArrayRead(Y,&y);
320:   VecGetArray(F,&f);
321:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
322:   VecRestoreArrayRead(Y,&y);
323:   VecRestoreArray(F,&f);
324:   /* Left hand side = ydot - f(y) */
325:   VecAYPX(F,-1.0,Ydot);
326:   return(0);
327: }

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

337:   VecGetArrayRead(Y,&y);
338:   value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
339:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
340:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
341:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
342:   VecRestoreArrayRead(Y,&y);
343:   return(0);
344: }

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

348: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
349: {
350:   PetscErrorCode    ierr;
351:   PetscScalar       *f;
352:   const PetscScalar *y;

355:   VecGetArrayRead(Y,&y);
356:   VecGetArray(F,&f);
357:   f[0] = (y[0]-t)/(y[0]+t);
358:   VecRestoreArrayRead(Y,&y);
359:   VecRestoreArray(F,&f);
360:   return(0);
361: }

363: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
364: {
365:   PetscErrorCode    ierr;
366:   const PetscScalar *y;
367:   PetscInt          row = 0,col = 0;
368:   PetscScalar       value;

371:   VecGetArrayRead(Y,&y);
372:   value = 2*t/((t+y[0])*(t+y[0]));
373:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
374:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
375:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
376:   VecRestoreArrayRead(Y,&y);
377:   return(0);
378: }

380: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
381: {
382:   PetscErrorCode    ierr;
383:   PetscScalar       *f;
384:   const PetscScalar *y;

387:   VecGetArrayRead(Y,&y);
388:   VecGetArray(F,&f);
389:   f[0] = (y[0]-t)/(y[0]+t);
390:   VecRestoreArrayRead(Y,&y);
391:   VecRestoreArray(F,&f);
392:   /* Left hand side = ydot - f(y) */
393:   VecAYPX(F,-1.0,Ydot);
394:   return(0);
395: }

397: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
398: {
399:   PetscErrorCode    ierr;
400:   const PetscScalar *y;
401:   PetscInt          row = 0,col = 0;
402:   PetscScalar       value;

405:   VecGetArrayRead(Y,&y);
406:   value = a - 2*t/((t+y[0])*(t+y[0]));
407:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
408:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
409:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
410:   VecRestoreArrayRead(Y,&y);
411:   return(0);
412: }

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

416: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
417: {
418:   PetscErrorCode    ierr;
419:   PetscScalar       *f;
420:   const PetscScalar *y;

423:   VecGetArrayRead(Y,&y);
424:   VecGetArray(F,&f);
425:   f[0] = 2.0*(y[0] - y[0]*y[1]);
426:   f[1] = -(y[1]-y[0]*y[1]);
427:   VecRestoreArrayRead(Y,&y);
428:   VecRestoreArray(F,&f);
429:   return(0);
430: }

432: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
433: {
434:   PetscErrorCode    ierr;
435:   PetscScalar       *f;
436:   const PetscScalar *y;

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

450: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
451: {
452:   PetscErrorCode    ierr;
453:   const PetscScalar *y;
454:   PetscInt          row[2] = {0,1};
455:   PetscScalar       value[2][2];

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

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

470: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
471: {
472:   PetscErrorCode    ierr;
473:   PetscScalar       *f;
474:   const PetscScalar *y;

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

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

494:   VecGetArrayRead(Y,&y);
495:   VecGetArray(F,&f);
496:   f[0] = -y[0] + y[1];
497:   f[1] = y[0] - 2.0*y[1] + y[2];
498:   f[2] = y[1] - y[2];
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_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
507: {
508:   PetscErrorCode    ierr;
509:   const PetscScalar *y;
510:   PetscInt          row[3] = {0,1,2};
511:   PetscScalar       value[3][3];

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

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

527: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
528: {
529:   PetscErrorCode    ierr;
530:   PetscScalar       *f;
531:   const PetscScalar *y;

534:   VecGetArrayRead(Y,&y);
535:   VecGetArray(F,&f);
536:   f[0] = -y[0];
537:   f[1] = y[0] - y[1]*y[1];
538:   f[2] = y[1]*y[1];
539:   VecRestoreArrayRead(Y,&y);
540:   VecRestoreArray(F,&f);
541:   return(0);
542: }

544: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
545: {
546:   PetscErrorCode    ierr;
547:   PetscScalar       *f;
548:   const PetscScalar *y;

551:   VecGetArrayRead(Y,&y);
552:   VecGetArray(F,&f);
553:   f[0] = -y[0];
554:   f[1] = y[0] - y[1]*y[1];
555:   f[2] = y[1]*y[1];
556:   VecRestoreArrayRead(Y,&y);
557:   VecRestoreArray(F,&f);
558:   /* Left hand side = ydot - f(y) */
559:   VecAYPX(F,-1.0,Ydot);
560:   return(0);
561: }

563: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
564: {
565:   PetscErrorCode    ierr;
566:   const PetscScalar *y;
567:   PetscInt          row[3] = {0,1,2};
568:   PetscScalar       value[3][3];

571:   VecGetArrayRead(Y,&y);
572:   value[0][0] = a + 1.0; value[0][1] = 0;             value[0][2] = 0;
573:   value[1][0] = -1.0;    value[1][1] = a + 2.0*y[1];  value[1][2] = 0;
574:   value[2][0] = 0;       value[2][1] = -2.0*y[1];     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 B4 */

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

591:   VecGetArrayRead(Y,&y);
592:   VecGetArray(F,&f);
593:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
594:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
595:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
596:   VecRestoreArrayRead(Y,&y);
597:   VecRestoreArray(F,&f);
598:   return(0);
599: }

601: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
602: {
603:   PetscErrorCode    ierr;
604:   PetscScalar       *f;
605:   const PetscScalar *y;

608:   VecGetArrayRead(Y,&y);
609:   VecGetArray(F,&f);
610:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
611:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
612:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
613:   VecRestoreArrayRead(Y,&y);
614:   VecRestoreArray(F,&f);
615:   /* Left hand side = ydot - f(y) */
616:   VecAYPX(F,-1.0,Ydot);
617:   return(0);
618: }

620: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
621: {
622:   PetscErrorCode    ierr;
623:   const PetscScalar *y;
624:   PetscInt          row[3] = {0,1,2};
625:   PetscScalar       value[3][3],fac,fac2;

628:   VecGetArrayRead(Y,&y);
629:   fac  = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
630:   fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
631:   value[0][0] = a + (y[1]*y[1]*y[2])*fac;
632:   value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
633:   value[0][2] = y[0]*fac2;
634:   value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
635:   value[1][1] = a + y[0]*y[0]*y[2]*fac;
636:   value[1][2] = y[1]*fac2;
637:   value[2][0] = -y[1]*y[1]*fac;
638:   value[2][1] = y[0]*y[1]*fac;
639:   value[2][2] = a;
640:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
641:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
642:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
643:   VecRestoreArrayRead(Y,&y);
644:   return(0);
645: }

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

649: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
650: {
651:   PetscErrorCode    ierr;
652:   PetscScalar       *f;
653:   const PetscScalar *y;

656:   VecGetArrayRead(Y,&y);
657:   VecGetArray(F,&f);
658:   f[0] = y[1]*y[2];
659:   f[1] = -y[0]*y[2];
660:   f[2] = -0.51*y[0]*y[1];
661:   VecRestoreArrayRead(Y,&y);
662:   VecRestoreArray(F,&f);
663:   return(0);
664: }

666: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
667: {
668:   PetscErrorCode    ierr;
669:   PetscScalar       *f;
670:   const PetscScalar *y;

673:   VecGetArrayRead(Y,&y);
674:   VecGetArray(F,&f);
675:   f[0] = y[1]*y[2];
676:   f[1] = -y[0]*y[2];
677:   f[2] = -0.51*y[0]*y[1];
678:   VecRestoreArrayRead(Y,&y);
679:   VecRestoreArray(F,&f);
680:   /* Left hand side = ydot - f(y) */
681:   VecAYPX(F,-1.0,Ydot);
682:   return(0);
683: }

685: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
686: {
687:   PetscErrorCode    ierr;
688:   const PetscScalar *y;
689:   PetscInt          row[3] = {0,1,2};
690:   PetscScalar       value[3][3];

693:   VecGetArrayRead(Y,&y);
694:   value[0][0] = a;          value[0][1] = -y[2];      value[0][2] = -y[1];
695:   value[1][0] = y[2];       value[1][1] = a;          value[1][2] = y[0];
696:   value[2][0] = 0.51*y[1];  value[2][1] = 0.51*y[0];  value[2][2] = a;
697:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
698:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
699:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
700:   VecRestoreArrayRead(Y,&y);
701:   return(0);
702: }


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

707: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
708: {
709:   PetscErrorCode    ierr;
710:   PetscScalar       *f;
711:   const PetscScalar *y;

714:   VecGetArrayRead(Y,&y);
715:   VecGetArray(F,&f);
716:   f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
717:   f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
718:   f[2] = 2.*t*y[3];
719:   f[3] = -2.*t*PetscLogScalar(y[0]);
720:   VecRestoreArrayRead(Y,&y);
721:   VecRestoreArray(F,&f);
722:   return(0);
723: }

725: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
726: {
727:   PetscErrorCode    ierr;
728:   const PetscScalar *y;
729:   PetscInt          row[4] = {0,1,2,3};
730:   PetscScalar       value[4][4];
731:   PetscScalar       m1,m2;
733:   VecGetArrayRead(Y,&y);
734:   m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
735:   m2=2.*t*PetscPowScalar(y[1],1./5.);
736:   value[0][0] = 0. ;        value[0][1] = m1; value[0][2] = 0.;  value[0][3] = m2;
737:   m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
738:   m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
739:   value[1][0] = 0.;        value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
740:   value[2][0] = 0.;        value[2][1] = 0.;  value[2][2] = 0.; value[2][3] = 2*t;
741:   value[3][0] = -2.*t/y[0];value[3][1] = 0.;  value[3][2] = 0.; value[3][3] = 0.;
742:   MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
743:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
744:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
745:   VecRestoreArrayRead(Y,&y);
746:   return(0);
747: }

749: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
750: {
751:   PetscErrorCode    ierr;
752:   PetscScalar       *f;
753:   const PetscScalar *y;

756:   VecGetArrayRead(Y,&y);
757:   VecGetArray(F,&f);
758:   f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
759:   f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
760:   f[2] = 2.*t*y[3];
761:   f[3] = -2.*t*PetscLogScalar(y[0]);
762:   VecRestoreArrayRead(Y,&y);
763:   VecRestoreArray(F,&f);
764:   /* Left hand side = ydot - f(y) */
765:   VecAYPX(F,-1.0,Ydot);
766:   return(0);
767: }

769: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
770: {
771:   PetscErrorCode    ierr;
772:   const PetscScalar *y;
773:   PetscInt          row[4] = {0,1,2,3};
774:   PetscScalar       value[4][4];
775:   PetscScalar       m1,m2;

778:   VecGetArrayRead(Y,&y);
779:   m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
780:   m2=2.*t*PetscPowScalar(y[1],1./5.);
781:   value[0][0] = a ;        value[0][1] = m1;  value[0][2] = 0.; value[0][3] = m2;
782:   m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
783:   m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
784:   value[1][0] = 0.;        value[1][1] = a ;  value[1][2] = m1; value[1][3] = m2;
785:   value[2][0] = 0.;        value[2][1] = 0.;  value[2][2] = a;  value[2][3] = 2*t;
786:   value[3][0] = -2.*t/y[0];value[3][1] = 0.;  value[3][2] = 0.; value[3][3] = a;
787:   MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
788:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
789:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
790:   VecRestoreArrayRead(Y,&y);
791:   return(0);
792: }


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

797: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
798: {
799:   PetscErrorCode    ierr;
800:   PetscScalar       *f;
801:   const PetscScalar *y;
802:   PetscInt          N,i;

805:   VecGetSize (Y,&N);
806:   VecGetArrayRead(Y,&y);
807:   VecGetArray(F,&f);
808:   f[0] = -y[0];
809:   for (i = 1; i < N-1; i++) {
810:     f[i] = y[i-1] - y[i];
811:   }
812:   f[N-1] = y[N-2];
813:   VecRestoreArrayRead(Y,&y);
814:   VecRestoreArray(F,&f);
815:   return(0);
816: }

818: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
819: {
820:   PetscErrorCode    ierr;
821:   PetscScalar       *f;
822:   const PetscScalar *y;
823:   PetscInt          N,i;

826:   VecGetSize (Y,&N);
827:   VecGetArrayRead(Y,&y);
828:   VecGetArray(F,&f);
829:   f[0] = -y[0];
830:   for (i = 1; i < N-1; i++) {
831:     f[i] = y[i-1] - y[i];
832:   }
833:   f[N-1] = y[N-2];
834:   VecRestoreArrayRead(Y,&y);
835:   VecRestoreArray(F,&f);
836:   /* Left hand side = ydot - f(y) */
837:   VecAYPX(F,-1.0,Ydot);
838:   return(0);
839: }

841: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
842: {
843:   PetscErrorCode    ierr;
844:   const PetscScalar *y;
845:   PetscInt          N,i,col[2];
846:   PetscScalar       value[2];

849:   VecGetSize (Y,&N);
850:   VecGetArrayRead(Y,&y);
851:   i = 0;
852:   value[0] = a+1; col[0] = 0;
853:   value[1] =  0;  col[1] = 1;
854:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
855:   for (i = 0; i < N; i++) {
856:     value[0] =  -1; col[0] = i-1;
857:     value[1] = a+1; col[1] = i;
858:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
859:   }
860:   i = N-1;
861:   value[0] = -1;  col[0] = N-2;
862:   value[1] = a;   col[1] = N-1;
863:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
864:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
865:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
866:   VecRestoreArrayRead(Y,&y);
867:   return(0);
868: }

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

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

880:   VecGetSize (Y,&N);
881:   VecGetArrayRead(Y,&y);
882:   VecGetArray(F,&f);
883:   f[0] = -y[0];
884:   for (i = 1; i < N-1; i++) {
885:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
886:   }
887:   f[N-1] = (PetscReal)(N-1)*y[N-2];
888:   VecRestoreArrayRead(Y,&y);
889:   VecRestoreArray(F,&f);
890:   return(0);
891: }

893: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
894: {
895:   PetscErrorCode    ierr;
896:   PetscScalar       *f;
897:   const PetscScalar *y;
898:   PetscInt          N,i;

901:   VecGetSize (Y,&N);
902:   VecGetArrayRead(Y,&y);
903:   VecGetArray(F,&f);
904:   f[0] = -y[0];
905:   for (i = 1; i < N-1; i++) {
906:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
907:   }
908:   f[N-1] = (PetscReal)(N-1)*y[N-2];
909:   VecRestoreArrayRead(Y,&y);
910:   VecRestoreArray(F,&f);
911:   /* Left hand side = ydot - f(y) */
912:   VecAYPX(F,-1.0,Ydot);
913:   return(0);
914: }

916: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
917: {
918:   PetscErrorCode    ierr;
919:   const PetscScalar *y;
920:   PetscInt          N,i,col[2];
921:   PetscScalar       value[2];

924:   VecGetSize (Y,&N);
925:   VecGetArrayRead(Y,&y);
926:   i = 0;
927:   value[0] = a+1;                 col[0] = 0;
928:   value[1] = 0;                   col[1] = 1;
929:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
930:   for (i = 0; i < N; i++) {
931:     value[0] = -(PetscReal) i;      col[0] = i-1;
932:     value[1] = a+(PetscReal)(i+1);  col[1] = i;
933:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
934:   }
935:   i = N-1;
936:   value[0] = -(PetscReal) (N-1);  col[0] = N-2;
937:   value[1] = a;                   col[1] = N-1;
938:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
939:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
940:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
941:   VecRestoreArrayRead(Y,&y);
942:   return(0);
943: }

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

947: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
948: {
949:   PetscErrorCode    ierr;
950:   PetscScalar       *f;
951:   const PetscScalar *y;
952:   PetscInt          N,i;

955:   VecGetSize (Y,&N);
956:   VecGetArrayRead(Y,&y);
957:   VecGetArray(F,&f);
958:   f[0] = -2.0*y[0] + y[1];
959:   for (i = 1; i < N-1; i++) {
960:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
961:   }
962:   f[N-1] = y[N-2] - 2.0*y[N-1];
963:   VecRestoreArrayRead(Y,&y);
964:   VecRestoreArray(F,&f);
965:   return(0);
966: }

968: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
969: {
970:   PetscErrorCode    ierr;
971:   PetscScalar       *f;
972:   const PetscScalar *y;
973:   PetscInt          N,i;

976:   VecGetSize (Y,&N);
977:   VecGetArrayRead(Y,&y);
978:   VecGetArray(F,&f);
979:   f[0] = -2.0*y[0] + y[1];
980:   for (i = 1; i < N-1; i++) {
981:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
982:   }
983:   f[N-1] = y[N-2] - 2.0*y[N-1];
984:   VecRestoreArrayRead(Y,&y);
985:   VecRestoreArray(F,&f);
986:   /* Left hand side = ydot - f(y) */
987:   VecAYPX(F,-1.0,Ydot);
988:   return(0);
989: }

991: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
992: {
993:   PetscErrorCode    ierr;
994:   const PetscScalar *y;
995:   PetscScalar       value[3];  
996:   PetscInt          N,i,col[3];

999:   VecGetSize (Y,&N);
1000:   VecGetArrayRead(Y,&y);
1001:   for (i = 0; i < N; i++) {
1002:     if (i == 0) {
1003:       value[0] = a+2;  col[0] = i;
1004:       value[1] =  -1;  col[1] = i+1;
1005:       value[2] =  0;   col[2] = i+2;
1006:     } else if (i == N-1) {
1007:       value[0] =  0;   col[0] = i-2;
1008:       value[1] =  -1;  col[1] = i-1;
1009:       value[2] = a+2;  col[2] = i;
1010:     } else {
1011:       value[0] = -1;   col[0] = i-1;
1012:       value[1] = a+2;  col[1] = i;
1013:       value[2] = -1;   col[2] = i+1;
1014:     }
1015:     MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
1016:   }
1017:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1018:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
1019:   VecRestoreArrayRead(Y,&y);
1020:   return(0);
1021: }

1023: /***************************************************************************/

1025: /* Sets the initial solution for the IVP and sets up the function pointers*/
1026: PetscErrorCode Initialize(Vec Y, void* s)
1027: {
1029:   char          *p = (char*) s;
1030:   PetscScalar   *y;
1031:   PetscReal     t0;
1032:   PetscInt      N = GetSize((const char *)s);
1033:   PetscBool     flg;

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

1135: /* Calculates the exact solution to problems that have one */
1136: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1137: {
1139:   char          *p = (char*) s;
1140:   PetscScalar   *y;

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

1178: /* Solves the specified ODE and computes the error if exact solution is available */
1179: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1180: {
1181:   PetscErrorCode  ierr;             /* Error code                             */
1182:   TS              ts;               /* time-integrator                        */
1183:   Vec             Y;                /* Solution vector                        */
1184:   Vec             Yex;              /* Exact solution                         */
1185:   PetscInt        N;                /* Size of the system of equations        */
1186:   TSType          time_scheme;      /* Type of time-integration scheme        */
1187:   Mat             Jac = NULL;       /* Jacobian matrix                        */
1188:   Vec             Yerr;             /* Auxiliary solution vector              */
1189:   PetscReal       err_norm;         /* Estimated error norm                   */
1190:   PetscReal       final_time;       /* Actual final time from the integrator  */
1191:   
1193:   N = GetSize((const char *)&ptype[0]);
1194:   if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1195:   VecCreate(PETSC_COMM_WORLD,&Y);
1196:   VecSetSizes(Y,N,PETSC_DECIDE);
1197:   VecSetUp(Y);
1198:   VecSet(Y,0);

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

1203:   /* Create and initialize the time-integrator                            */
1204:   TSCreate(PETSC_COMM_WORLD,&ts);
1205:   /* Default time integration options                                     */
1206:   TSSetType(ts,TSRK);
1207:   TSSetMaxSteps(ts,maxiter);
1208:   TSSetMaxTime(ts,tfinal);
1209:   TSSetTimeStep(ts,dt);
1210:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1211:   /* Read command line options for time integration                       */
1212:   TSSetFromOptions(ts);
1213:   /* Set solution vector                                                  */
1214:   TSSetSolution(ts,Y);
1215:   /* Specify left/right-hand side functions                               */
1216:   TSGetType(ts,&time_scheme);
1217:   
1218:   if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1219:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1220:     TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1221:     MatCreate(PETSC_COMM_WORLD,&Jac);
1222:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1223:     MatSetFromOptions(Jac);
1224:     MatSetUp(Jac);
1225:     TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1226:   } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1227:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1228:     /* and its Jacobian function                                                 */
1229:     TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1230:     MatCreate(PETSC_COMM_WORLD,&Jac);
1231:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1232:     MatSetFromOptions(Jac);
1233:     MatSetUp(Jac);
1234:     TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1235:   }

1237:   /* Solve */
1238:   TSSolve(ts,Y);
1239:   TSGetTime(ts,&final_time);

1241:   /* Get the estimated error, if available */
1242:   VecDuplicate(Y,&Yerr);
1243:   VecZeroEntries(Yerr);
1244:   TSGetTimeError(ts,0,&Yerr);
1245:   VecNorm(Yerr,NORM_2,&err_norm);
1246:   VecDestroy(&Yerr);
1247:   PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);

1249:   /* Exact solution */
1250:   VecDuplicate(Y,&Yex);
1251:   if(PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1252:     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);
1253:   }
1254:   ExactSolution(Yex,&ptype[0],final_time,exact_flag);

1256:   /* Calculate Error */
1257:   VecAYPX(Yex,-1.0,Y);
1258:   VecNorm(Yex,NORM_2,error);
1259:   *error = PetscSqrtReal(((*error)*(*error))/N);

1261:   /* Clean up and finalize */
1262:   MatDestroy(&Jac);
1263:   TSDestroy(&ts);
1264:   VecDestroy(&Yex);
1265:   VecDestroy(&Y);

1267:   return(0);
1268: }

1270: int main(int argc, char **argv)
1271: {
1272:   PetscErrorCode  ierr;                       /* Error code                                           */
1273:   char            ptype[256] = "hull1972a1";  /* Problem specification                                */
1274:   PetscInt        n_refine   = 1;             /* Number of refinement levels for convergence analysis */
1275:   PetscReal       refine_fac = 2.0;           /* Refinement factor for dt                             */
1276:   PetscReal       dt_initial = 0.01;          /* Initial default value of dt                          */
1277:   PetscReal       dt;
1278:   PetscReal       tfinal     = 20.0;          /* Final time for the time-integration                  */
1279:   PetscInt        maxiter    = 100000;        /* Maximum number of time-integration iterations        */
1280:   PetscReal       *error;                     /* Array to store the errors for convergence analysis   */
1281:   PetscMPIInt     size;                      /* No of processors                                     */
1282:   PetscBool       flag;                       /* Flag denoting availability of exact solution         */
1283:   PetscInt        r;

1285:   /* Initialize program */
1286:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

1288:   /* Check if running with only 1 proc */
1289:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
1290:   if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

1292:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);
1293:   PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1294:   PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1295:   PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1296:   PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1297:   PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1298:   PetscOptionsEnd();

1300:   PetscMalloc1(n_refine,&error);
1301:   for (r = 0,dt = dt_initial; r < n_refine; r++) {
1302:     error[r] = 0;
1303:     if (r > 0) dt /= refine_fac;

1305:     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]));
1306:     SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1307:     if (flag) {
1308:       /* If exact solution available for the specified ODE */
1309:       if (r > 0) {
1310:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1311:         PetscPrintf(PETSC_COMM_WORLD,"Error           = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1312:       } else {
1313:         PetscPrintf(PETSC_COMM_WORLD,"Error           = %E.\n",error[r]);
1314:       }
1315:     }
1316:   }
1317:   PetscFree(error);
1318:   PetscFinalize();
1319:   return ierr;
1320: }

1322: /*TEST

1324:     test:
1325:       suffix: 2
1326:       args: -ts_type glee -final_time 5 -ts_adapt_type none
1327:       timeoutfactor: 3
1328:       requires: !single

1330:     test:
1331:       suffix: 3
1332:       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
1333:       timeoutfactor: 3
1334:       requires: !single

1336:     test:
1337:       suffix: 4
1338:       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
1339:       timeoutfactor: 3
1340:       requires: !single !__float128

1342: TEST*/