Actual source code: ex31.c

petsc-3.8.4 2018-03-24
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>
 45: #if defined(PETSC_HAVE_STRING_H)
 46: #include <string.h>             /* strcmp */
 47: #endif

 49: /* Function declarations */
 50: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
 51: PetscErrorCode (*RHSJacobian) (TS,PetscReal,Vec,Mat,Mat,void*);
 52: PetscErrorCode (*IFunction)   (TS,PetscReal,Vec,Vec,Vec,void*);
 53: PetscErrorCode (*IJacobian)   (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

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

 77: /****************************************************************/

 79: /* Problem specific functions */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

290:   VecGetArrayRead(Y,&y);
291:   VecGetArray(F,&f);
292:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
293:   VecRestoreArrayRead(Y,&y);
294:   VecRestoreArray(F,&f);
295:   return(0);
296: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

461:   VecGetArrayRead(Y,&y);
462:   value[0][0] = a - 2.0*(1.0-y[1]);    value[0][1] = 2.0*y[0];
463:   value[1][0] = -y[1];                 value[1][1] = a + 1.0 - y[0];
464:   MatSetValues(A,2,&row[0],2,&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 B2 */

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

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

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

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

509: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
510: {
511:   PetscErrorCode    ierr;
512:   const PetscScalar *y;
513:   PetscInt          row[3] = {0,1,2};
514:   PetscScalar       value[3][3];

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

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

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

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

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

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

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

574:   VecGetArrayRead(Y,&y);
575:   value[0][0] = a + 1.0; value[0][1] = 0;             value[0][2] = 0;
576:   value[1][0] = -1.0;    value[1][1] = a + 2.0*y[1];  value[1][2] = 0;
577:   value[2][0] = 0;       value[2][1] = -2.0*y[1];     value[2][2] = a;
578:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
579:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
580:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
581:   VecRestoreArrayRead(Y,&y);
582:   return(0);
583: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1026: /***************************************************************************/

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

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

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

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

1181: /* Solves the specified ODE and computes the error if exact solution is available */
1182: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1183: {
1184:   PetscErrorCode  ierr;             /* Error code                             */
1185:   TS              ts;               /* time-integrator                        */
1186:   Vec             Y;                /* Solution vector                        */
1187:   Vec             Yex;              /* Exact solution                         */
1188:   PetscInt        N;                /* Size of the system of equations        */
1189:   TSType          time_scheme;      /* Type of time-integration scheme        */
1190:   Mat             Jac = NULL;       /* Jacobian matrix                        */
1191:   Vec             Yerr;             /* Auxiliary solution vector              */
1192:   PetscReal       err_norm;         /* Estimated error norm                   */

1195:   N = GetSize((const char *)&ptype[0]);
1196:   if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1197:   VecCreate(PETSC_COMM_WORLD,&Y);
1198:   VecSetSizes(Y,N,PETSC_DECIDE);
1199:   VecSetUp(Y);
1200:   VecSet(Y,0);

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

1205:   /* Create and initialize the time-integrator                            */
1206:   TSCreate(PETSC_COMM_WORLD,&ts);
1207:   /* Default time integration options                                     */
1208:   TSSetType(ts,TSRK);
1209:   TSSetMaxSteps(ts,maxiter);
1210:   TSSetMaxTime(ts,tfinal);
1211:   TSSetTimeStep(ts,dt);
1212:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1213:   /* Read command line options for time integration                       */
1214:   TSSetFromOptions(ts);
1215:   /* Set solution vector                                                  */
1216:   TSSetSolution(ts,Y);
1217:   /* Specify left/right-hand side functions                               */
1218:   TSGetType(ts,&time_scheme);

1220:   if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1221:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1222:     TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1223:     MatCreate(PETSC_COMM_WORLD,&Jac);
1224:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1225:     MatSetFromOptions(Jac);
1226:     MatSetUp(Jac);
1227:     TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1228:   } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1229:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1230:     /* and its Jacobian function                                                 */
1231:     TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1232:     MatCreate(PETSC_COMM_WORLD,&Jac);
1233:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1234:     MatSetFromOptions(Jac);
1235:     MatSetUp(Jac);
1236:     TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1237:   }

1239:   /* Solve */
1240:   TSSolve(ts,Y);

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

1250:   /* Exact solution */
1251:   VecDuplicate(Y,&Yex);
1252:   ExactSolution(Yex,&ptype[0],tfinal,exact_flag);

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

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

1265:   return(0);
1266: }

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

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

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

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

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

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