Actual source code: ex31.c

      1: static char help[] = "Solves the ordinary differential equations (IVPs) using explicit and implicit time-integration methods.\n";

      3: /*

      5:   Concepts:   TS
      6:   Useful command line parameters:
      7:   -problem <hull1972a1>: choose which problem to solve (see references
      8:                       for complete listing of problems).
      9:   -ts_type <euler>: specify time-integrator
     10:   -ts_adapt_type <basic>: specify time-step adapting (none,basic,advanced)
     11:   -refinement_levels <1>: number of refinement levels for convergence analysis
     12:   -refinement_factor <2.0>: factor to refine time step size by for convergence analysis
     13:   -dt <0.01>: specify time step (initial time step for convergence analysis)

     15: */

     17: /*
     18: List of cases and their names in the code:-
     19:   From Hull, T.E., Enright, W.H., Fellen, B.M., and Sedgwick, A.E.,
     20:       "Comparing Numerical Methods for Ordinary Differential
     21:        Equations", SIAM J. Numer. Anal., 9(4), 1972, pp. 603 - 635
     22:     A1 -> "hull1972a1" (exact solution available)
     23:     A2 -> "hull1972a2" (exact solution available)
     24:     A3 -> "hull1972a3" (exact solution available)
     25:     A4 -> "hull1972a4" (exact solution available)
     26:     A5 -> "hull1972a5"
     27:     B1 -> "hull1972b1"
     28:     B2 -> "hull1972b2"
     29:     B3 -> "hull1972b3"
     30:     B4 -> "hull1972b4"
     31:     B5 -> "hull1972b5"
     32:     C1 -> "hull1972c1"
     33:     C2 -> "hull1972c2"
     34:     C3 -> "hull1972c3"
     35:     C4 -> "hull1972c4"

     37:  From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
     38:        https://arxiv.org/abs/1503.05166, 2016

     40:     Kulikov2013I -> "kulik2013i"

     42: */

     44: #include <petscts.h>

     46: /* Function declarations */
     47: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
     48: PetscErrorCode (*RHSJacobian) (TS,PetscReal,Vec,Mat,Mat,void*);
     49: PetscErrorCode (*IFunction)   (TS,PetscReal,Vec,Vec,Vec,void*);
     50: PetscErrorCode (*IJacobian)   (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

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

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

     76: /* Problem specific functions */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1249:   /* Exact solution */
    1250:   VecDuplicate(Y,&Yex);
    1251:   if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
    1252:     PetscPrintf(PETSC_COMM_WORLD,"Note: There is a difference between the prescribed final time %g and the actual final time, %g.\n",(double)tfinal,(double)final_time);
    1253:   }
    1254:   ExactSolution(Yex,&ptype[0],final_time,exact_flag);

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

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

    1267:   return(0);
    1268: }

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

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

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

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

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

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

    1322: /*TEST

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

    1330:     test:
    1331:       suffix: 3
    1332:       args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor  -ts_max_steps 50  -problem hull1972a3 -ts_adapt_glee_use_local 1
    1333:       timeoutfactor: 3
    1334:       requires: !single

    1336:     test:
    1337:       suffix: 4
    1338:       args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor  -ts_max_steps 50  -problem hull1972a3  -ts_max_reject 100 -ts_adapt_glee_use_local 0
    1339:       timeoutfactor: 3
    1340:       requires: !single !__float128

    1342: TEST*/