Actual source code: ex31.c

petsc-master 2015-08-27
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"
 36: */

 38: #include <petscts.h>
 39: #if defined(PETSC_HAVE_STRING_H)
 40: #include <string.h>             /* strcmp */
 41: #endif

 43: /* Function declarations */
 44: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
 45: PetscErrorCode (*IFunction)   (TS,PetscReal,Vec,Vec,Vec,void*);
 46: PetscErrorCode (*IJacobian)   (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

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

 71: /****************************************************************/

 73: /* Problem specific functions */

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

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

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

 96: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
 97: {
 98:   PetscErrorCode    ierr;
 99:   const PetscScalar *y;
100:   PetscScalar       *f;

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

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

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

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

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

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

152: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
153: {
154:   PetscErrorCode    ierr;
155:   PetscScalar       *f;
156:   const PetscScalar *y;

159:   VecGetArrayRead(Y,&y);
160:   VecGetArray(F,&f);
161:   f[0] = -0.5*y[0]*y[0]*y[0];
162:   VecRestoreArrayRead(Y,&y);
163:   VecRestoreArray(F,&f);
164:   /* Left hand side = ydot - f(y) */
165:   VecAYPX(F,-1.0,Ydot);
166:   return(0);
167: }

171: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
172: {
173:   PetscErrorCode    ierr;
174:   const PetscScalar *y;
175:   PetscInt          row = 0,col = 0;
176:   PetscScalar       value;

179:   VecGetArrayRead(Y,&y);
180:   value = a + 0.5*3.0*y[0]*y[0];
181:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
182:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
183:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
184:   VecRestoreArrayRead(Y,&y);
185:   return(0);
186: }

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

192: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
193: {
194:   PetscErrorCode    ierr;
195:   const PetscScalar *y;
196:   PetscScalar       *f;

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

209: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
210: {
211:   PetscErrorCode    ierr;
212:   PetscScalar       *f;
213:   const PetscScalar *y;

216:   VecGetArrayRead(Y,&y);
217:   VecGetArray(F,&f);
218:   f[0] = y[0]*PetscCosReal(t);
219:   VecRestoreArrayRead(Y,&y);
220:   VecRestoreArray(F,&f);
221:   /* Left hand side = ydot - f(y) */
222:   VecAYPX(F,-1.0,Ydot);
223:   return(0);
224: }

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

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

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

248: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, 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] = (0.25*y[0])*(1.0-0.05*y[0]);
258:   VecRestoreArrayRead(Y,&y);
259:   VecRestoreArray(F,&f);
260:   return(0);
261: }

265: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
266: {
267:   PetscErrorCode    ierr;
268:   PetscScalar       *f;
269:   const PetscScalar *y;

272:   VecGetArrayRead(Y,&y);
273:   VecGetArray(F,&f);
274:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
275:   VecRestoreArrayRead(Y,&y);
276:   VecRestoreArray(F,&f);
277:   /* Left hand side = ydot - f(y) */
278:   VecAYPX(F,-1.0,Ydot);
279:   return(0);
280: }

284: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
285: {
286:   PetscErrorCode    ierr;
287:   const PetscScalar *y;
288:   PetscInt          row = 0,col = 0;
289:   PetscScalar       value;

292:   VecGetArrayRead(Y,&y);
293:   value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
294:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
295:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
296:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
297:   VecRestoreArrayRead(Y,&y);
298:   return(0);
299: }

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

305: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
306: {
307:   PetscErrorCode    ierr;
308:   PetscScalar       *f;
309:   const PetscScalar *y;

312:   VecGetArrayRead(Y,&y);
313:   VecGetArray(F,&f);
314:   f[0] = (y[0]-t)/(y[0]+t);
315:   VecRestoreArrayRead(Y,&y);
316:   VecRestoreArray(F,&f);
317:   return(0);
318: }

322: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
323: {
324:   PetscErrorCode    ierr;
325:   PetscScalar       *f;
326:   const PetscScalar *y;

329:   VecGetArrayRead(Y,&y);
330:   VecGetArray(F,&f);
331:   f[0] = (y[0]-t)/(y[0]+t);
332:   VecRestoreArrayRead(Y,&y);
333:   VecRestoreArray(F,&f);
334:   /* Left hand side = ydot - f(y) */
335:   VecAYPX(F,-1.0,Ydot);
336:   return(0);
337: }

341: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
342: {
343:   PetscErrorCode    ierr;
344:   const PetscScalar *y;
345:   PetscInt          row = 0,col = 0;
346:   PetscScalar       value;

349:   VecGetArrayRead(Y,&y);
350:   value = a - 2*t/((t+y[0])*(t+y[0]));
351:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
352:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
353:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
354:   VecRestoreArrayRead(Y,&y);
355:   return(0);
356: }

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

362: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
363: {
364:   PetscErrorCode    ierr;
365:   PetscScalar       *f;
366:   const PetscScalar *y;

369:   VecGetArrayRead(Y,&y);
370:   VecGetArray(F,&f);
371:   f[0] = 2.0*(y[0] - y[0]*y[1]);
372:   f[1] = -(y[1]-y[0]*y[1]);
373:   VecRestoreArrayRead(Y,&y);
374:   VecRestoreArray(F,&f);
375:   return(0);
376: }

380: PetscErrorCode IFunction_Hull1972B1(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] = 2.0*(y[0] - y[0]*y[1]);
390:   f[1] = -(y[1]-y[0]*y[1]);
391:   VecRestoreArrayRead(Y,&y);
392:   VecRestoreArray(F,&f);
393:   /* Left hand side = ydot - f(y) */
394:   VecAYPX(F,-1.0,Ydot);
395:   return(0);
396: }

400: PetscErrorCode IJacobian_Hull1972B1(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[2] = {0,1};
405:   PetscScalar       value[2][2];

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

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

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

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

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

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

462: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
463: {
464:   PetscErrorCode    ierr;
465:   const PetscScalar *y;
466:   PetscInt          row[3] = {0,1,2};
467:   PetscScalar       value[3][3];

470:   VecGetArrayRead(Y,&y);
471:   value[0][0] = a + 1.0;  value[0][1] = -1.0;     value[0][2] = 0;
472:   value[1][0] = -1.0;     value[1][1] = a + 2.0;  value[1][2] = -1.0;
473:   value[2][0] = 0;        value[2][1] = -1.0;     value[2][2] = a + 1.0;
474:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
475:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
476:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
477:   VecRestoreArrayRead(Y,&y);
478:   return(0);
479: }

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

485: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
486: {
487:   PetscErrorCode    ierr;
488:   PetscScalar       *f;
489:   const PetscScalar *y;

492:   VecGetArrayRead(Y,&y);
493:   VecGetArray(F,&f);
494:   f[0] = -y[0];
495:   f[1] = y[0] - y[1]*y[1];
496:   f[2] = y[1]*y[1];
497:   VecRestoreArrayRead(Y,&y);
498:   VecRestoreArray(F,&f);
499:   return(0);
500: }

504: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
505: {
506:   PetscErrorCode    ierr;
507:   PetscScalar       *f;
508:   const PetscScalar *y;

511:   VecGetArrayRead(Y,&y);
512:   VecGetArray(F,&f);
513:   f[0] = -y[0];
514:   f[1] = y[0] - y[1]*y[1];
515:   f[2] = y[1]*y[1];
516:   VecRestoreArrayRead(Y,&y);
517:   VecRestoreArray(F,&f);
518:   /* Left hand side = ydot - f(y) */
519:   VecAYPX(F,-1.0,Ydot);
520:   return(0);
521: }

525: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
526: {
527:   PetscErrorCode    ierr;
528:   const PetscScalar *y;
529:   PetscInt          row[3] = {0,1,2};
530:   PetscScalar       value[3][3];

533:   VecGetArrayRead(Y,&y);
534:   value[0][0] = a + 1.0; value[0][1] = 0;             value[0][2] = 0;
535:   value[1][0] = -1.0;    value[1][1] = a + 2.0*y[1];  value[1][2] = 0;
536:   value[2][0] = 0;       value[2][1] = -2.0*y[1];     value[2][2] = a;
537:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
538:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
539:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
540:   VecRestoreArrayRead(Y,&y);
541:   return(0);
542: }

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

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

555:   VecGetArrayRead(Y,&y);
556:   VecGetArray(F,&f);
557:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
558:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
559:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
560:   VecRestoreArrayRead(Y,&y);
561:   VecRestoreArray(F,&f);
562:   return(0);
563: }

567: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
568: {
569:   PetscErrorCode    ierr;
570:   PetscScalar       *f;
571:   const PetscScalar *y;

574:   VecGetArrayRead(Y,&y);
575:   VecGetArray(F,&f);
576:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
577:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
578:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
579:   VecRestoreArrayRead(Y,&y);
580:   VecRestoreArray(F,&f);
581:   /* Left hand side = ydot - f(y) */
582:   VecAYPX(F,-1.0,Ydot);
583:   return(0);
584: }

588: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
589: {
590:   PetscErrorCode    ierr;
591:   const PetscScalar *y;
592:   PetscInt          row[3] = {0,1,2};
593:   PetscScalar       value[3][3],fac,fac2;

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

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

619: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
620: {
621:   PetscErrorCode    ierr;
622:   PetscScalar       *f;
623:   const PetscScalar *y;

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

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

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

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

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

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

682: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
683: {
684:   PetscErrorCode    ierr;
685:   PetscScalar       *f;
686:   const PetscScalar *y;
687:   PetscInt          N,i;

690:   VecGetSize (Y,&N);
691:   VecGetArrayRead(Y,&y);
692:   VecGetArray(F,&f);
693:   f[0] = -y[0];
694:   for (i = 1; i < N-1; i++) {
695:     f[i] = y[i-1] - y[i];
696:   }
697:   f[N-1] = y[N-2];
698:   VecRestoreArrayRead(Y,&y);
699:   VecRestoreArray(F,&f);
700:   return(0);
701: }

705: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
706: {
707:   PetscErrorCode    ierr;
708:   PetscScalar       *f;
709:   const PetscScalar *y;
710:   PetscInt          N,i;

713:   VecGetSize (Y,&N);
714:   VecGetArrayRead(Y,&y);
715:   VecGetArray(F,&f);
716:   f[0] = -y[0];
717:   for (i = 1; i < N-1; i++) {
718:     f[i] = y[i-1] - y[i];
719:   }
720:   f[N-1] = y[N-2];
721:   VecRestoreArrayRead(Y,&y);
722:   VecRestoreArray(F,&f);
723:   /* Left hand side = ydot - f(y) */
724:   VecAYPX(F,-1.0,Ydot);
725:   return(0);
726: }

730: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
731: {
732:   PetscErrorCode    ierr;
733:   const PetscScalar *y;
734:   PetscInt          N,i,col[2];
735:   PetscScalar       value[2];

738:   VecGetSize (Y,&N);
739:   VecGetArrayRead(Y,&y);
740:   i = 0;
741:   value[0] = a+1; col[0] = 0;
742:   value[1] =  0;  col[1] = 1;
743:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
744:   for (i = 0; i < N; i++) {
745:     value[0] =  -1; col[0] = i-1;
746:     value[1] = a+1; col[1] = i;
747:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
748:   }
749:   i = N-1;
750:   value[0] = -1;  col[0] = N-2;
751:   value[1] = a;   col[1] = N-1;
752:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
753:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
754:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
755:   VecRestoreArrayRead(Y,&y);
756:   return(0);
757: }

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

763: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
764: {
765:   PetscErrorCode    ierr;
766:   const PetscScalar *y;
767:   PetscScalar       *f;
768:   PetscInt          N,i;

771:   VecGetSize (Y,&N);
772:   VecGetArrayRead(Y,&y);
773:   VecGetArray(F,&f);
774:   f[0] = -y[0];
775:   for (i = 1; i < N-1; i++) {
776:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
777:   }
778:   f[N-1] = (PetscReal)(N-1)*y[N-2];
779:   VecRestoreArrayRead(Y,&y);
780:   VecRestoreArray(F,&f);
781:   return(0);
782: }

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

794:   VecGetSize (Y,&N);
795:   VecGetArrayRead(Y,&y);
796:   VecGetArray(F,&f);
797:   f[0] = -y[0];
798:   for (i = 1; i < N-1; i++) {
799:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
800:   }
801:   f[N-1] = (PetscReal)(N-1)*y[N-2];
802:   VecRestoreArrayRead(Y,&y);
803:   VecRestoreArray(F,&f);
804:   /* Left hand side = ydot - f(y) */
805:   VecAYPX(F,-1.0,Ydot);
806:   return(0);
807: }

811: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
812: {
813:   PetscErrorCode    ierr;
814:   const PetscScalar *y;
815:   PetscInt          N,i,col[2];
816:   PetscScalar       value[2];

819:   VecGetSize (Y,&N);
820:   VecGetArrayRead(Y,&y);
821:   i = 0;
822:   value[0] = a+1;                 col[0] = 0;
823:   value[1] = 0;                   col[1] = 1;
824:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
825:   for (i = 0; i < N; i++) {
826:     value[0] = -(PetscReal) i;      col[0] = i-1;
827:     value[1] = a+(PetscReal)(i+1);  col[1] = i;
828:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
829:   }
830:   i = N-1;
831:   value[0] = -(PetscReal) (N-1);  col[0] = N-2;
832:   value[1] = a;                   col[1] = N-1;
833:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
834:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
835:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
836:   VecRestoreArrayRead(Y,&y);
837:   return(0);
838: }

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

844: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
845: {
846:   PetscErrorCode    ierr;
847:   PetscScalar       *f;
848:   const PetscScalar *y;
849:   PetscInt          N,i;

852:   VecGetSize (Y,&N);
853:   VecGetArrayRead(Y,&y);
854:   VecGetArray(F,&f);
855:   f[0] = -2.0*y[0] + y[1];
856:   for (i = 1; i < N-1; i++) {
857:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
858:   }
859:   f[N-1] = y[N-2] - 2.0*y[N-1];
860:   VecRestoreArrayRead(Y,&y);
861:   VecRestoreArray(F,&f);
862:   return(0);
863: }

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

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

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

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

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

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

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

1026: /* Calculates the exact solution to problems that have one */
1027: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1028: {
1030:   char          *p = (char*) s;
1031:   PetscScalar   *y;

1034:   if (!strcmp(p,"hull1972a1")) {
1035:     VecGetArray(Y,&y);
1036:     y[0] = PetscExpReal(-t);
1037:     *flag = PETSC_TRUE;
1038:     VecRestoreArray(Y,&y);
1039:   } else if (!strcmp(p,"hull1972a2")) {
1040:     VecGetArray(Y,&y);
1041:     y[0] = 1.0/PetscSqrtReal(t+1);
1042:     *flag = PETSC_TRUE;
1043:     VecRestoreArray(Y,&y);
1044:   } else if (!strcmp(p,"hull1972a3")) {
1045:     VecGetArray(Y,&y);
1046:     y[0] = PetscExpReal(PetscSinReal(t));
1047:     *flag = PETSC_TRUE;
1048:     VecRestoreArray(Y,&y);
1049:   } else if (!strcmp(p,"hull1972a4")) {
1050:     VecGetArray(Y,&y);
1051:     y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1052:     *flag = PETSC_TRUE;
1053:     VecRestoreArray(Y,&y);
1054:   } else {
1055:     VecSet(Y,0);
1056:     *flag = PETSC_FALSE;
1057:   }
1058:   return(0);
1059: }

1063: /* Solves the specified ODE and computes the error if exact solution is available */
1064: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1065: {
1066:   PetscErrorCode  ierr;             /* Error code                             */
1067:   TS              ts;               /* time-integrator                        */
1068:   Vec             Y;                /* Solution vector                        */
1069:   Vec             Yex;              /* Exact solution                         */
1070:   PetscInt        N;                /* Size of the system of equations        */
1071:   TSType          time_scheme;      /* Type of time-integration scheme        */
1072:   Mat             Jac = NULL;       /* Jacobian matrix                        */

1075:   N = GetSize((const char *)&ptype[0]);
1076:   if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1077:   VecCreate(PETSC_COMM_WORLD,&Y);
1078:   VecSetSizes(Y,N,PETSC_DECIDE);
1079:   VecSetUp(Y);
1080:   VecSet(Y,0);

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

1085:   /* Create and initialize the time-integrator                            */
1086:   TSCreate(PETSC_COMM_WORLD,&ts);
1087:   /* Default time integration options                                     */
1088:   TSSetType(ts,TSEULER);
1089:   TSSetDuration(ts,maxiter,tfinal);
1090:   TSSetInitialTimeStep(ts,0.0,dt);
1091:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1092:   /* Read command line options for time integration                       */
1093:   TSSetFromOptions(ts);
1094:   /* Set solution vector                                                  */
1095:   TSSetSolution(ts,Y);
1096:   /* Specify left/right-hand side functions                               */
1097:   TSGetType(ts,&time_scheme);
1098:   if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP))) {
1099:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1100:     TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1101:   } else if ((!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSARKIMEX))) {
1102:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1103:     /* and its Jacobian function                                                 */
1104:     TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1105:     MatCreate(PETSC_COMM_WORLD,&Jac);
1106:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1107:     MatSetFromOptions(Jac);
1108:     MatSetUp(Jac);
1109:     TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1110:   }

1112:   /* Solve */
1113:   TSSolve(ts,Y);

1115:   /* Exact solution */
1116:   VecDuplicate(Y,&Yex);
1117:   ExactSolution(Yex,&ptype[0],tfinal,exact_flag);

1119:   /* Calculate Error */
1120:   VecAYPX(Yex,-1.0,Y);
1121:   VecNorm(Yex,NORM_2,error);
1122:   *error = PetscSqrtReal(((*error)*(*error))/N);

1124:   /* Clean up and finalize */
1125:   MatDestroy(&Jac);
1126:   TSDestroy(&ts);
1127:   VecDestroy(&Yex);
1128:   VecDestroy(&Y);

1130:   return(0);
1131: }

1135: int main(int argc, char **argv)
1136: {
1137:   PetscErrorCode  ierr;                       /* Error code                                           */
1138:   char            ptype[256] = "hull1972a1";  /* Problem specification                                */
1139:   PetscInt        n_refine   = 1;             /* Number of refinement levels for convergence analysis */
1140:   PetscReal       refine_fac = 2.0;           /* Refinement factor for dt                             */
1141:   PetscReal       dt_initial = 0.01;          /* Initial default value of dt                          */
1142:   PetscReal       dt;
1143:   PetscReal       tfinal     = 20.0;          /* Final time for the time-integration                  */
1144:   PetscInt        maxiter    = 100000;        /* Maximum number of time-integration iterations        */
1145:   PetscReal       *error;                     /* Array to store the errors for convergence analysis   */
1146:   PetscMPIInt     size;                      /* No of processors                                     */
1147:   PetscBool       flag;                       /* Flag denoting availability of exact solution         */
1148:   PetscInt        r;

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

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

1157:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
1158:   PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1159:   PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1160:   PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1161:   PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1162:   PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1163:   PetscOptionsEnd();

1165:   PetscMalloc1(n_refine,&error);
1166:   for (r = 0,dt = dt_initial; r < n_refine; r++) {
1167:     error[r] = 0;
1168:     if (r > 0) dt /= refine_fac;

1170:     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]));
1171:     SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1172:     if (flag) {
1173:       /* If exact solution available for the specified ODE */
1174:       if (r > 0) {
1175:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1176:         PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f\n.",(double)error[r],(double)conv_rate);
1177:       } else {
1178:         PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1179:       }
1180:     }
1181:   }
1182:   PetscFree(error);

1184:   /* Exit */
1185:   PetscFinalize();
1186:   return(0);
1187: }