Actual source code: ex31.c

petsc-3.5.4 2015-05-23
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: {
 82:   PetscScalar    *y,*f;

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

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

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

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

121:   VecGetArray(Y,&y);
122:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
123:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
124:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
125:   VecRestoreArray(Y,&y);
126:   return(0);
127: }

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

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

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

149: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
150: {
152:   PetscScalar    *y,*f;

155:   VecGetArray(Y,&y);
156:   VecGetArray(F,&f);
157:   f[0] = -0.5*y[0]*y[0]*y[0];
158:   VecRestoreArray(Y,&y);
159:   VecRestoreArray(F,&f);
160:   /* Left hand side = ydot - f(y) */
161:   VecAYPX(F,-1.0,Ydot);
162:   return(0);
163: }

167: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
168: {
170:   PetscScalar    *y;
171:   PetscInt       row = 0,col = 0;
172:   PetscScalar    value;

175:   VecGetArray(Y,&y);
176:   value = a + 0.5*3.0*y[0]*y[0];
177:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
178:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
179:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
180:   VecRestoreArray(Y,&y);
181:   return(0);
182: }

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

188: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
189: {
191:   PetscScalar    *y,*f;

194:   VecGetArray(Y,&y);
195:   VecGetArray(F,&f);
196:   f[0] = y[0]*PetscCosReal(t);
197:   VecRestoreArray(Y,&y);
198:   VecRestoreArray(F,&f);
199:   return(0);
200: }

204: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
205: {
207:   PetscScalar    *y,*f;

210:   VecGetArray(Y,&y);
211:   VecGetArray(F,&f);
212:   f[0] = y[0]*PetscCosReal(t);
213:   VecRestoreArray(Y,&y);
214:   VecRestoreArray(F,&f);
215:   /* Left hand side = ydot - f(y) */
216:   VecAYPX(F,-1.0,Ydot);
217:   return(0);
218: }

222: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
223: {
225:   PetscScalar    *y;
226:   PetscInt       row = 0,col = 0;
227:   PetscScalar    value = a - PetscCosReal(t);

230:   VecGetArray(Y,&y);
231:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
232:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
233:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
234:   VecRestoreArray(Y,&y);
235:   return(0);
236: }

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

242: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
243: {
245:   PetscScalar    *y,*f;

248:   VecGetArray(Y,&y);
249:   VecGetArray(F,&f);
250:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
251:   VecRestoreArray(Y,&y);
252:   VecRestoreArray(F,&f);
253:   return(0);
254: }

258: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
259: {
261:   PetscScalar    *y,*f;

264:   VecGetArray(Y,&y);
265:   VecGetArray(F,&f);
266:   f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
267:   VecRestoreArray(Y,&y);
268:   VecRestoreArray(F,&f);
269:   /* Left hand side = ydot - f(y) */
270:   VecAYPX(F,-1.0,Ydot);
271:   return(0);
272: }

276: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
277: {
279:   PetscScalar    *y;
280:   PetscInt       row = 0,col = 0;
281:   PetscScalar    value;

284:   VecGetArray(Y,&y);
285:   value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
286:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
287:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
288:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);

290:   VecRestoreArray(Y,&y);
291:   return(0);
292: }

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

298: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
299: {
301:   PetscScalar    *y,*f;

304:   VecGetArray(Y,&y);
305:   VecGetArray(F,&f);
306:   f[0] = (y[0]-t)/(y[0]+t);
307:   VecRestoreArray(Y,&y);
308:   VecRestoreArray(F,&f);
309:   return(0);
310: }

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

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

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

340:   VecGetArray(Y,&y);
341:   value = a - 2*t/((t+y[0])*(t+y[0]));
342:   MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
343:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
344:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
345:   VecRestoreArray(Y,&y);
346:   return(0);
347: }

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

353: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
354: {
356:   PetscScalar    *y,*f;

359:   VecGetArray(Y,&y);
360:   VecGetArray(F,&f);
361:   f[0] = 2.0*(y[0] - y[0]*y[1]);
362:   f[1] = -(y[1]-y[0]*y[1]);
363:   VecRestoreArray(Y,&y);
364:   VecRestoreArray(F,&f);
365:   return(0);
366: }

370: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
371: {
373:   PetscScalar    *y,*f;

376:   VecGetArray(Y,&y);
377:   VecGetArray(F,&f);
378:   f[0] = 2.0*(y[0] - y[0]*y[1]);
379:   f[1] = -(y[1]-y[0]*y[1]);
380:   VecRestoreArray(Y,&y);
381:   VecRestoreArray(F,&f);
382:   /* Left hand side = ydot - f(y) */
383:   VecAYPX(F,-1.0,Ydot);
384:   return(0);
385: }

389: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
390: {
392:   PetscScalar    *y;
393:   PetscInt       row[2] = {0,1};
394:   PetscScalar    value[2][2];

397:   VecGetArray(Y,&y);
398:   value[0][0] = a - 2.0*(1.0-y[1]);    value[0][1] = 2.0*y[0];
399:   value[1][0] = -y[1];                 value[1][1] = a + 1.0 - y[0];
400:   MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
401:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
402:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
403:   VecRestoreArray(Y,&y);
404:   return(0);
405: }

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

411: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
412: {
414:   PetscScalar    *y,*f;

417:   VecGetArray(Y,&y);
418:   VecGetArray(F,&f);
419:   f[0] = -y[0] + y[1];
420:   f[1] = y[0] - 2.0*y[1] + y[2];
421:   f[2] = y[1] - y[2];
422:   VecRestoreArray(Y,&y);
423:   VecRestoreArray(F,&f);
424:   return(0);
425: }

429: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
430: {
431:   PetscErrorCode  ierr;
432:   PetscScalar    *y,*f;

435:   VecGetArray(Y,&y);
436:   VecGetArray(F,&f);
437:   f[0] = -y[0] + y[1];
438:   f[1] = y[0] - 2.0*y[1] + y[2];
439:   f[2] = y[1] - y[2];
440:   VecRestoreArray(Y,&y);
441:   VecRestoreArray(F,&f);
442:   /* Left hand side = ydot - f(y) */
443:   VecAYPX(F,-1.0,Ydot);
444:   return(0);
445: }

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

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

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

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

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

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

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

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

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

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

533: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
534: {
536:   PetscScalar    *y,*f;

539:   VecGetArray(Y,&y);
540:   VecGetArray(F,&f);
541:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
542:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
543:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
544:   VecRestoreArray(Y,&y);
545:   VecRestoreArray(F,&f);
546:   return(0);
547: }

551: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
552: {
553:   PetscErrorCode  ierr;
554:   PetscScalar    *y,*f;

557:   VecGetArray(Y,&y);
558:   VecGetArray(F,&f);
559:   f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
560:   f[1] =  y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
561:   f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
562:   VecRestoreArray(Y,&y);
563:   VecRestoreArray(F,&f);
564:   /* Left hand side = ydot - f(y) */
565:   VecAYPX(F,-1.0,Ydot);
566:   return(0);
567: }

571: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
572: {
574:   PetscScalar    *y;
575:   PetscInt       row[3] = {0,1,2};
576:   PetscScalar    value[3][3],fac,fac2;

579:   VecGetArray(Y,&y);
580:   fac  = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
581:   fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
582:   value[0][0] = a + (y[1]*y[1]*y[2])*fac;
583:   value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
584:   value[0][2] = y[0]*fac2;
585:   value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
586:   value[1][1] = a + y[0]*y[0]*y[2]*fac;
587:   value[1][2] = y[1]*fac2;
588:   value[2][0] = -y[1]*y[1]*fac;
589:   value[2][1] = y[0]*y[1]*fac;
590:   value[2][2] = a;
591:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
592:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
593:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
594:   VecRestoreArray(Y,&y);
595:   return(0);
596: }

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

602: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
603: {
605:   PetscScalar    *y,*f;

608:   VecGetArray(Y,&y);
609:   VecGetArray(F,&f);
610:   f[0] = y[1]*y[2];
611:   f[1] = -y[0]*y[2];
612:   f[2] = -0.51*y[0]*y[1];
613:   VecRestoreArray(Y,&y);
614:   VecRestoreArray(F,&f);
615:   return(0);
616: }

620: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
621: {
623:   PetscScalar    *y,*f;

626:   VecGetArray(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:   VecRestoreArray(Y,&y);
632:   VecRestoreArray(F,&f);
633:   /* Left hand side = ydot - f(y) */
634:   VecAYPX(F,-1.0,Ydot);
635:   return(0);
636: }

640: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
641: {
643:   PetscScalar    *y;
644:   PetscInt       row[3] = {0,1,2};
645:   PetscScalar    value[3][3];

648:   VecGetArray(Y,&y);
649:   value[0][0] = a;          value[0][1] = -y[2];      value[0][2] = -y[1];
650:   value[1][0] = y[2];       value[1][1] = a;          value[1][2] = y[0];
651:   value[2][0] = 0.51*y[1];  value[2][1] = 0.51*y[0];  value[2][2] = a;
652:   MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
653:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
654:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
655:   VecRestoreArray(Y,&y);
656:   return(0);
657: }

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

663: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
664: {
666:   PetscScalar    *y,*f;
667:   PetscInt       N,i;

670:   VecGetSize (Y,&N);
671:   VecGetArray(Y,&y);
672:   VecGetArray(F,&f);
673:   f[0] = -y[0];
674:   for (i = 1; i < N-1; i++) {
675:     f[i] = y[i-1] - y[i];
676:   }
677:   f[N-1] = y[N-2];
678:   VecRestoreArray(Y,&y);
679:   VecRestoreArray(F,&f);
680:   return(0);
681: }

685: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
686: {
688:   PetscScalar    *y,*f;
689:   PetscInt       N,i;

692:   VecGetSize (Y,&N);
693:   VecGetArray(Y,&y);
694:   VecGetArray(F,&f);
695:   f[0] = -y[0];
696:   for (i = 1; i < N-1; i++) {
697:     f[i] = y[i-1] - y[i];
698:   }
699:   f[N-1] = y[N-2];
700:   VecRestoreArray(Y,&y);
701:   VecRestoreArray(F,&f);
702:   /* Left hand side = ydot - f(y) */
703:   VecAYPX(F,-1.0,Ydot);
704:   return(0);
705: }

709: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
710: {
711:   PetscErrorCode  ierr;
712:   PetscScalar    *y;
713:   PetscInt        N,i,col[2];
714:   PetscScalar     value[2];

717:   VecGetSize (Y,&N);
718:   VecGetArray(Y,&y);
719:   i = 0;
720:   value[0] = a+1; col[0] = 0;
721:   value[1] =  0;  col[1] = 1;
722:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
723:   for (i = 0; i < N; i++) {
724:     value[0] =  -1; col[0] = i-1;
725:     value[1] = a+1; col[1] = i;
726:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
727:   }
728:   i = N-1;
729:   value[0] = -1;  col[0] = N-2;
730:   value[1] = a;   col[1] = N-1;
731:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
732:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
733:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
734:   VecRestoreArray(Y,&y);
735:   return(0);
736: }

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

742: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
743: {
744:   PetscErrorCode  ierr;
745:   PetscScalar    *y,*f;
746:   PetscInt        N,i;

749:   VecGetSize (Y,&N);
750:   VecGetArray(Y,&y);
751:   VecGetArray(F,&f);
752:   f[0] = -y[0];
753:   for (i = 1; i < N-1; i++) {
754:     f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
755:   }
756:   f[N-1] = (PetscReal)(N-1)*y[N-2];
757:   VecRestoreArray(Y,&y);
758:   VecRestoreArray(F,&f);
759:   return(0);
760: }

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

771:   VecGetSize (Y,&N);
772:   VecGetArray(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:   VecRestoreArray(Y,&y);
780:   VecRestoreArray(F,&f);
781:   /* Left hand side = ydot - f(y) */
782:   VecAYPX(F,-1.0,Ydot);
783:   return(0);
784: }

788: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
789: {
790:   PetscErrorCode  ierr;
791:   PetscScalar    *y;
792:   PetscInt        N,i,col[2];
793:   PetscScalar     value[2];

796:   VecGetSize (Y,&N);
797:   VecGetArray(Y,&y);
798:   i = 0;
799:   value[0] = a+1;                 col[0] = 0;
800:   value[1] = 0;                   col[1] = 1;
801:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
802:   for (i = 0; i < N; i++) {
803:     value[0] = -(PetscReal) i;      col[0] = i-1;
804:     value[1] = a+(PetscReal)(i+1);  col[1] = i;
805:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
806:   }
807:   i = N-1;
808:   value[0] = -(PetscReal) (N-1);  col[0] = N-2;
809:   value[1] = a;                   col[1] = N-1;
810:   MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
811:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
812:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
813:   VecRestoreArray(Y,&y);
814:   return(0);
815: }

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

821: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
822: {
824:   PetscScalar    *y,*f;
825:   PetscInt       N,i;

828:   VecGetSize (Y,&N);
829:   VecGetArray(Y,&y);
830:   VecGetArray(F,&f);
831:   f[0] = -2.0*y[0] + y[1];
832:   for (i = 1; i < N-1; i++) {
833:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
834:   }
835:   f[N-1] = y[N-2] - 2.0*y[N-1];
836:   VecRestoreArray(Y,&y);
837:   VecRestoreArray(F,&f);
838:   return(0);
839: }

843: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
844: {
846:   PetscScalar    *y,*f;
847:   PetscInt       N,i;

850:   VecGetSize (Y,&N);
851:   VecGetArray(Y,&y);
852:   VecGetArray(F,&f);
853:   f[0] = -2.0*y[0] + y[1];
854:   for (i = 1; i < N-1; i++) {
855:     f[i] = y[i-1] - 2.0*y[i] + y[i+1];
856:   }
857:   f[N-1] = y[N-2] - 2.0*y[N-1];
858:   VecRestoreArray(Y,&y);
859:   VecRestoreArray(F,&f);
860:   /* Left hand side = ydot - f(y) */
861:   VecAYPX(F,-1.0,Ydot);
862:   return(0);
863: }

867: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
868: {
870:   PetscScalar    *y,value[3];
871:   PetscInt       N,i,col[3];

874:   VecGetSize (Y,&N);
875:   VecGetArray(Y,&y);
876:   for (i = 0; i < N; i++) {
877:     if (i == 0) {
878:       value[0] = a+2;  col[0] = i;
879:       value[1] =  -1;  col[1] = i+1;
880:       value[2] =  0;   col[2] = i+2;
881:     } else if (i == N-1) {
882:       value[0] =  0;   col[0] = i-2;
883:       value[1] =  -1;  col[1] = i-1;
884:       value[2] = a+2;  col[2] = i;
885:     } else {
886:       value[0] = -1;   col[0] = i-1;
887:       value[1] = a+2;  col[1] = i;
888:       value[2] = -1;   col[2] = i+1;
889:     }
890:     MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
891:   }
892:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
893:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
894:   VecRestoreArray(Y,&y);
895:   return(0);
896: }

898: /***************************************************************************/

902: /* Sets the initial solution for the IVP and sets up the function pointers*/
903: PetscErrorCode Initialize(Vec Y, void* s)
904: {
906:   char          *p = (char*) s;
907:   PetscScalar   *y;
908:   PetscInt      N = GetSize((const char *)s);
909:   PetscBool     flg;

912:   VecZeroEntries(Y);
913:   VecGetArray(Y,&y);
914:   if (!strcmp(p,"hull1972a1")) {
915:     y[0] = 1.0;
916:     RHSFunction = RHSFunction_Hull1972A1;
917:     IFunction   = IFunction_Hull1972A1;
918:     IJacobian   = IJacobian_Hull1972A1;
919:   } else if (!strcmp(p,"hull1972a2")) {
920:     y[0] = 1.0;
921:     RHSFunction = RHSFunction_Hull1972A2;
922:     IFunction   = IFunction_Hull1972A2;
923:     IJacobian   = IJacobian_Hull1972A2;
924:   } else if (!strcmp(p,"hull1972a3")) {
925:     y[0] = 1.0;
926:     RHSFunction = RHSFunction_Hull1972A3;
927:     IFunction   = IFunction_Hull1972A3;
928:     IJacobian   = IJacobian_Hull1972A3;
929:   } else if (!strcmp(p,"hull1972a4")) {
930:     y[0] = 1.0;
931:     RHSFunction = RHSFunction_Hull1972A4;
932:     IFunction   = IFunction_Hull1972A4;
933:     IJacobian   = IJacobian_Hull1972A4;
934:   } else if (!strcmp(p,"hull1972a5")) {
935:     y[0] = 4.0;
936:     RHSFunction = RHSFunction_Hull1972A5;
937:     IFunction   = IFunction_Hull1972A5;
938:     IJacobian   = IJacobian_Hull1972A5;
939:   } else if (!strcmp(p,"hull1972b1")) {
940:     y[0] = 1.0;
941:     y[1] = 3.0;
942:     RHSFunction = RHSFunction_Hull1972B1;
943:     IFunction   = IFunction_Hull1972B1;
944:     IJacobian   = IJacobian_Hull1972B1;
945:   } else if (!strcmp(p,"hull1972b2")) {
946:     y[0] = 2.0;
947:     y[1] = 0.0;
948:     y[2] = 1.0;
949:     RHSFunction = RHSFunction_Hull1972B2;
950:     IFunction   = IFunction_Hull1972B2;
951:     IJacobian   = IJacobian_Hull1972B2;
952:   } else if (!strcmp(p,"hull1972b3")) {
953:     y[0] = 1.0;
954:     y[1] = 0.0;
955:     y[2] = 0.0;
956:     RHSFunction = RHSFunction_Hull1972B3;
957:     IFunction   = IFunction_Hull1972B3;
958:     IJacobian   = IJacobian_Hull1972B3;
959:   } else if (!strcmp(p,"hull1972b4")) {
960:     y[0] = 3.0;
961:     y[1] = 0.0;
962:     y[2] = 0.0;
963:     RHSFunction = RHSFunction_Hull1972B4;
964:     IFunction   = IFunction_Hull1972B4;
965:     IJacobian   = IJacobian_Hull1972B4;
966:   } else if (!strcmp(p,"hull1972b5")) {
967:     y[0] = 0.0;
968:     y[1] = 1.0;
969:     y[2] = 1.0;
970:     RHSFunction = RHSFunction_Hull1972B5;
971:     IFunction   = IFunction_Hull1972B5;
972:     IJacobian   = IJacobian_Hull1972B5;
973:   } else if (!strcmp(p,"hull1972c1")) {
974:     y[0] = 1.0;
975:     RHSFunction = RHSFunction_Hull1972C1;
976:     IFunction   = IFunction_Hull1972C1;
977:     IJacobian   = IJacobian_Hull1972C1;
978:   } else if (!strcmp(p,"hull1972c2")) {
979:     y[0] = 1.0;
980:     RHSFunction = RHSFunction_Hull1972C2;
981:     IFunction   = IFunction_Hull1972C2;
982:     IJacobian   = IJacobian_Hull1972C2;
983:   } else if ((!strcmp(p,"hull1972c3"))
984:            ||(!strcmp(p,"hull1972c4"))){
985:     y[0] = 1.0;
986:     RHSFunction = RHSFunction_Hull1972C34;
987:     IFunction   = IFunction_Hull1972C34;
988:     IJacobian   = IJacobian_Hull1972C34;
989:   }
990:   PetscOptionsGetScalarArray(NULL,"-yinit",y,&N,&flg);
991:   if ((N != GetSize((const char*)s)) && flg) {
992:     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));
993:   }
994:   VecRestoreArray(Y,&y);
995:   return(0);
996: }

1000: /* Calculates the exact solution to problems that have one */
1001: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1002: {
1004:   char          *p = (char*) s;
1005:   PetscScalar   *y;

1008:   if (!strcmp(p,"hull1972a1")) {
1009:     VecGetArray(Y,&y);
1010:     y[0] = PetscExpReal(-t);
1011:     *flag = PETSC_TRUE;
1012:     VecRestoreArray(Y,&y);
1013:   } else if (!strcmp(p,"hull1972a2")) {
1014:     VecGetArray(Y,&y);
1015:     y[0] = 1.0/PetscSqrtReal(t+1);
1016:     *flag = PETSC_TRUE;
1017:     VecRestoreArray(Y,&y);
1018:   } else if (!strcmp(p,"hull1972a3")) {
1019:     VecGetArray(Y,&y);
1020:     y[0] = PetscExpReal(PetscSinReal(t));
1021:     *flag = PETSC_TRUE;
1022:     VecRestoreArray(Y,&y);
1023:   } else if (!strcmp(p,"hull1972a4")) {
1024:     VecGetArray(Y,&y);
1025:     y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1026:     *flag = PETSC_TRUE;
1027:     VecRestoreArray(Y,&y);
1028:   } else {
1029:     VecSet(Y,0);
1030:     *flag = PETSC_FALSE;
1031:   }
1032:   return(0);
1033: }

1037: /* Solves the specified ODE and computes the error if exact solution is available */
1038: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1039: {
1040:   PetscErrorCode  ierr;             /* Error code                             */
1041:   TS              ts;               /* time-integrator                        */
1042:   Vec             Y;                /* Solution vector                        */
1043:   Vec             Yex;              /* Exact solution                         */
1044:   PetscInt        N;                /* Size of the system of equations        */
1045:   TSType          time_scheme;      /* Type of time-integration scheme        */
1046:   Mat             Jac = NULL;       /* Jacobian matrix                        */

1049:   N = GetSize((const char *)&ptype[0]);
1050:   if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1051:   VecCreate(PETSC_COMM_WORLD,&Y);
1052:   VecSetSizes(Y,N,PETSC_DECIDE);
1053:   VecSetUp(Y);
1054:   VecSet(Y,0);

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

1059:   /* Create and initialize the time-integrator                            */
1060:   TSCreate(PETSC_COMM_WORLD,&ts);
1061:   /* Default time integration options                                     */
1062:   TSSetType(ts,TSEULER);
1063:   TSSetDuration(ts,maxiter,tfinal);
1064:   TSSetInitialTimeStep(ts,0.0,dt);
1065:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1066:   /* Read command line options for time integration                       */
1067:   TSSetFromOptions(ts);
1068:   /* Set solution vector                                                  */
1069:   TSSetSolution(ts,Y);
1070:   /* Specify left/right-hand side functions                               */
1071:   TSGetType(ts,&time_scheme);
1072:   if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP))) {
1073:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1074:     TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1075:   } else if ((!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSARKIMEX))) {
1076:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1077:     /* and its Jacobian function                                                 */
1078:     TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1079:     MatCreate(PETSC_COMM_WORLD,&Jac);
1080:     MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1081:     MatSetFromOptions(Jac);
1082:     MatSetUp(Jac);
1083:     TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1084:   }

1086:   /* Solve */
1087:   TSSolve(ts,Y);

1089:   /* Exact solution */
1090:   VecDuplicate(Y,&Yex);
1091:   ExactSolution(Yex,&ptype[0],tfinal,exact_flag);

1093:   /* Calculate Error */
1094:   VecAYPX(Yex,-1.0,Y);
1095:   VecNorm(Yex,NORM_2,error);
1096:   *error = PetscSqrtReal(((*error)*(*error))/N);

1098:   /* Clean up and finalize */
1099:   MatDestroy(&Jac);
1100:   TSDestroy(&ts);
1101:   VecDestroy(&Yex);
1102:   VecDestroy(&Y);

1104:   return(0);
1105: }

1109: int main(int argc, char **argv)
1110: {
1111:   PetscErrorCode  ierr;                       /* Error code                                           */
1112:   char            ptype[256] = "hull1972a1";  /* Problem specification                                */
1113:   PetscInt        n_refine   = 1;             /* Number of refinement levels for convergence analysis */
1114:   PetscReal       refine_fac = 2.0;           /* Refinement factor for dt                             */
1115:   PetscReal       dt_initial = 0.01;          /* Initial default value of dt                          */
1116:   PetscReal       dt;
1117:   PetscReal       tfinal     = 20.0;          /* Final time for the time-integration                  */
1118:   PetscInt        maxiter    = 100000;        /* Maximum number of time-integration iterations        */
1119:   PetscReal       *error;                     /* Array to store the errors for convergence analysis   */
1120:   PetscMPIInt     size;                      /* No of processors                                     */
1121:   PetscBool       flag;                       /* Flag denoting availability of exact solution         */
1122:   PetscInt        r;

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

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

1131:   PetscOptionsString("-problem","Problem specification","<hull1972a1>",
1132:                             ptype,ptype,sizeof(ptype),NULL);
1133:   PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis",
1134:                          "<1>",n_refine,&n_refine,NULL);
1135:   PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",
1136:                           refine_fac,&refine_fac,NULL);
1137:   PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)",
1138:                           "<0.01>",dt_initial,&dt_initial,NULL);
1139:   PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",
1140:                           tfinal,&tfinal,NULL);

1142:   PetscMalloc1(n_refine,&error);
1143:   for (r = 0,dt = dt_initial; r < n_refine; r++) {
1144:     error[r] = 0;
1145:     if (r > 0) dt /= refine_fac;

1147:     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]));
1148:     SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1149:     if (flag) {
1150:       /* If exact solution available for the specified ODE */
1151:       if (r > 0) {
1152:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1153:         PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f\n.",(double)error[r],(double)conv_rate);
1154:       } else {
1155:         PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1156:       }
1157:     }
1158:   }
1159:   PetscFree(error);

1161:   /* Exit */
1162:   PetscFinalize();
1163:   return(0);
1164: }