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

     14: */

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

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

     39:     Kulikov2013I -> "kulik2013i"

     41: */

     43: #include <petscts.h>

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

     51: /* Returns the size of the system of equations depending on problem specification */
     52: PetscErrorCode GetSize(const char *p, PetscInt *sz)
     53: {
     54:   PetscFunctionBeginUser;
     55:   if (!strcmp(p, "hull1972a1") || !strcmp(p, "hull1972a2") || !strcmp(p, "hull1972a3") || !strcmp(p, "hull1972a4") || !strcmp(p, "hull1972a5")) *sz = 1;
     56:   else if (!strcmp(p, "hull1972b1")) *sz = 2;
     57:   else if (!strcmp(p, "hull1972b2") || !strcmp(p, "hull1972b3") || !strcmp(p, "hull1972b4") || !strcmp(p, "hull1972b5")) *sz = 3;
     58:   else if (!strcmp(p, "kulik2013i")) *sz = 4;
     59:   else if (!strcmp(p, "hull1972c1") || !strcmp(p, "hull1972c2") || !strcmp(p, "hull1972c3")) *sz = 10;
     60:   else if (!strcmp(p, "hull1972c4")) *sz = 52;
     61:   else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Incorrect problem name %s", p);
     62:   PetscFunctionReturn(PETSC_SUCCESS);
     63: }

     65: /****************************************************************/

     67: /* Problem specific functions */

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

     71: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
     72: {
     73:   PetscScalar       *f;
     74:   const PetscScalar *y;

     76:   PetscFunctionBeginUser;
     77:   PetscCall(VecGetArrayRead(Y, &y));
     78:   PetscCall(VecGetArray(F, &f));
     79:   f[0] = -y[0];
     80:   PetscCall(VecRestoreArrayRead(Y, &y));
     81:   PetscCall(VecRestoreArray(F, &f));
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
     86: {
     87:   const PetscScalar *y;
     88:   PetscInt           row = 0, col = 0;
     89:   PetscScalar        value = -1.0;

     91:   PetscFunctionBeginUser;
     92:   PetscCall(VecGetArrayRead(Y, &y));
     93:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
     94:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     95:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     96:   PetscCall(VecRestoreArrayRead(Y, &y));
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    101: {
    102:   const PetscScalar *y;
    103:   PetscScalar       *f;

    105:   PetscFunctionBeginUser;
    106:   PetscCall(VecGetArrayRead(Y, &y));
    107:   PetscCall(VecGetArray(F, &f));
    108:   f[0] = -y[0];
    109:   PetscCall(VecRestoreArrayRead(Y, &y));
    110:   PetscCall(VecRestoreArray(F, &f));
    111:   /* Left hand side = ydot - f(y) */
    112:   PetscCall(VecAYPX(F, -1.0, Ydot));
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

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

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

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

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

    138:   PetscFunctionBeginUser;
    139:   PetscCall(VecGetArrayRead(Y, &y));
    140:   PetscCall(VecGetArray(F, &f));
    141:   f[0] = -0.5 * y[0] * y[0] * y[0];
    142:   PetscCall(VecRestoreArrayRead(Y, &y));
    143:   PetscCall(VecRestoreArray(F, &f));
    144:   PetscFunctionReturn(PETSC_SUCCESS);
    145: }

    147: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
    148: {
    149:   const PetscScalar *y;
    150:   PetscInt           row = 0, col = 0;
    151:   PetscScalar        value;

    153:   PetscFunctionBeginUser;
    154:   PetscCall(VecGetArrayRead(Y, &y));
    155:   value = -0.5 * 3.0 * y[0] * y[0];
    156:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    157:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    158:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    159:   PetscCall(VecRestoreArrayRead(Y, &y));
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    164: {
    165:   PetscScalar       *f;
    166:   const PetscScalar *y;

    168:   PetscFunctionBeginUser;
    169:   PetscCall(VecGetArrayRead(Y, &y));
    170:   PetscCall(VecGetArray(F, &f));
    171:   f[0] = -0.5 * y[0] * y[0] * y[0];
    172:   PetscCall(VecRestoreArrayRead(Y, &y));
    173:   PetscCall(VecRestoreArray(F, &f));
    174:   /* Left hand side = ydot - f(y) */
    175:   PetscCall(VecAYPX(F, -1.0, Ydot));
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    180: {
    181:   const PetscScalar *y;
    182:   PetscInt           row = 0, col = 0;
    183:   PetscScalar        value;

    185:   PetscFunctionBeginUser;
    186:   PetscCall(VecGetArrayRead(Y, &y));
    187:   value = a + 0.5 * 3.0 * y[0] * y[0];
    188:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    189:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    190:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    191:   PetscCall(VecRestoreArrayRead(Y, &y));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

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

    197: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    198: {
    199:   const PetscScalar *y;
    200:   PetscScalar       *f;

    202:   PetscFunctionBeginUser;
    203:   PetscCall(VecGetArrayRead(Y, &y));
    204:   PetscCall(VecGetArray(F, &f));
    205:   f[0] = y[0] * PetscCosReal(t);
    206:   PetscCall(VecRestoreArrayRead(Y, &y));
    207:   PetscCall(VecRestoreArray(F, &f));
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
    212: {
    213:   const PetscScalar *y;
    214:   PetscInt           row = 0, col = 0;
    215:   PetscScalar        value = PetscCosReal(t);

    217:   PetscFunctionBeginUser;
    218:   PetscCall(VecGetArrayRead(Y, &y));
    219:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    220:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    221:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    222:   PetscCall(VecRestoreArrayRead(Y, &y));
    223:   PetscFunctionReturn(PETSC_SUCCESS);
    224: }

    226: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    227: {
    228:   PetscScalar       *f;
    229:   const PetscScalar *y;

    231:   PetscFunctionBeginUser;
    232:   PetscCall(VecGetArrayRead(Y, &y));
    233:   PetscCall(VecGetArray(F, &f));
    234:   f[0] = y[0] * PetscCosReal(t);
    235:   PetscCall(VecRestoreArrayRead(Y, &y));
    236:   PetscCall(VecRestoreArray(F, &f));
    237:   /* Left hand side = ydot - f(y) */
    238:   PetscCall(VecAYPX(F, -1.0, Ydot));
    239:   PetscFunctionReturn(PETSC_SUCCESS);
    240: }

    242: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    243: {
    244:   const PetscScalar *y;
    245:   PetscInt           row = 0, col = 0;
    246:   PetscScalar        value = a - PetscCosReal(t);

    248:   PetscFunctionBeginUser;
    249:   PetscCall(VecGetArrayRead(Y, &y));
    250:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    251:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    252:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    253:   PetscCall(VecRestoreArrayRead(Y, &y));
    254:   PetscFunctionReturn(PETSC_SUCCESS);
    255: }

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

    259: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    260: {
    261:   PetscScalar       *f;
    262:   const PetscScalar *y;

    264:   PetscFunctionBeginUser;
    265:   PetscCall(VecGetArrayRead(Y, &y));
    266:   PetscCall(VecGetArray(F, &f));
    267:   f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
    268:   PetscCall(VecRestoreArrayRead(Y, &y));
    269:   PetscCall(VecRestoreArray(F, &f));
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
    274: {
    275:   const PetscScalar *y;
    276:   PetscInt           row = 0, col = 0;
    277:   PetscScalar        value;

    279:   PetscFunctionBeginUser;
    280:   PetscCall(VecGetArrayRead(Y, &y));
    281:   value = 0.25 * (1.0 - 0.05 * y[0]) - (0.25 * y[0]) * 0.05;
    282:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    283:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    284:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    285:   PetscCall(VecRestoreArrayRead(Y, &y));
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    289: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    290: {
    291:   PetscScalar       *f;
    292:   const PetscScalar *y;

    294:   PetscFunctionBeginUser;
    295:   PetscCall(VecGetArrayRead(Y, &y));
    296:   PetscCall(VecGetArray(F, &f));
    297:   f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
    298:   PetscCall(VecRestoreArrayRead(Y, &y));
    299:   PetscCall(VecRestoreArray(F, &f));
    300:   /* Left hand side = ydot - f(y) */
    301:   PetscCall(VecAYPX(F, -1.0, Ydot));
    302:   PetscFunctionReturn(PETSC_SUCCESS);
    303: }

    305: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    306: {
    307:   const PetscScalar *y;
    308:   PetscInt           row = 0, col = 0;
    309:   PetscScalar        value;

    311:   PetscFunctionBeginUser;
    312:   PetscCall(VecGetArrayRead(Y, &y));
    313:   value = a - 0.25 * (1.0 - 0.05 * y[0]) + (0.25 * y[0]) * 0.05;
    314:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    315:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    316:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    317:   PetscCall(VecRestoreArrayRead(Y, &y));
    318:   PetscFunctionReturn(PETSC_SUCCESS);
    319: }

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

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

    328:   PetscFunctionBeginUser;
    329:   PetscCall(VecGetArrayRead(Y, &y));
    330:   PetscCall(VecGetArray(F, &f));
    331:   f[0] = (y[0] - t) / (y[0] + t);
    332:   PetscCall(VecRestoreArrayRead(Y, &y));
    333:   PetscCall(VecRestoreArray(F, &f));
    334:   PetscFunctionReturn(PETSC_SUCCESS);
    335: }

    337: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
    338: {
    339:   const PetscScalar *y;
    340:   PetscInt           row = 0, col = 0;
    341:   PetscScalar        value;

    343:   PetscFunctionBeginUser;
    344:   PetscCall(VecGetArrayRead(Y, &y));
    345:   value = 2 * t / ((t + y[0]) * (t + y[0]));
    346:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    347:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    348:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    349:   PetscCall(VecRestoreArrayRead(Y, &y));
    350:   PetscFunctionReturn(PETSC_SUCCESS);
    351: }

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

    358:   PetscFunctionBeginUser;
    359:   PetscCall(VecGetArrayRead(Y, &y));
    360:   PetscCall(VecGetArray(F, &f));
    361:   f[0] = (y[0] - t) / (y[0] + t);
    362:   PetscCall(VecRestoreArrayRead(Y, &y));
    363:   PetscCall(VecRestoreArray(F, &f));
    364:   /* Left hand side = ydot - f(y) */
    365:   PetscCall(VecAYPX(F, -1.0, Ydot));
    366:   PetscFunctionReturn(PETSC_SUCCESS);
    367: }

    369: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    370: {
    371:   const PetscScalar *y;
    372:   PetscInt           row = 0, col = 0;
    373:   PetscScalar        value;

    375:   PetscFunctionBeginUser;
    376:   PetscCall(VecGetArrayRead(Y, &y));
    377:   value = a - 2 * t / ((t + y[0]) * (t + y[0]));
    378:   PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
    379:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    380:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    381:   PetscCall(VecRestoreArrayRead(Y, &y));
    382:   PetscFunctionReturn(PETSC_SUCCESS);
    383: }

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

    387: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    388: {
    389:   PetscScalar       *f;
    390:   const PetscScalar *y;

    392:   PetscFunctionBeginUser;
    393:   PetscCall(VecGetArrayRead(Y, &y));
    394:   PetscCall(VecGetArray(F, &f));
    395:   f[0] = 2.0 * (y[0] - y[0] * y[1]);
    396:   f[1] = -(y[1] - y[0] * y[1]);
    397:   PetscCall(VecRestoreArrayRead(Y, &y));
    398:   PetscCall(VecRestoreArray(F, &f));
    399:   PetscFunctionReturn(PETSC_SUCCESS);
    400: }

    402: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    403: {
    404:   PetscScalar       *f;
    405:   const PetscScalar *y;

    407:   PetscFunctionBeginUser;
    408:   PetscCall(VecGetArrayRead(Y, &y));
    409:   PetscCall(VecGetArray(F, &f));
    410:   f[0] = 2.0 * (y[0] - y[0] * y[1]);
    411:   f[1] = -(y[1] - y[0] * y[1]);
    412:   PetscCall(VecRestoreArrayRead(Y, &y));
    413:   PetscCall(VecRestoreArray(F, &f));
    414:   /* Left hand side = ydot - f(y) */
    415:   PetscCall(VecAYPX(F, -1.0, Ydot));
    416:   PetscFunctionReturn(PETSC_SUCCESS);
    417: }

    419: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    420: {
    421:   const PetscScalar *y;
    422:   PetscInt           row[2] = {0, 1};
    423:   PetscScalar        value[2][2];

    425:   PetscFunctionBeginUser;
    426:   PetscCall(VecGetArrayRead(Y, &y));
    427:   value[0][0] = a - 2.0 * (1.0 - y[1]);
    428:   value[0][1] = 2.0 * y[0];
    429:   value[1][0] = -y[1];
    430:   value[1][1] = a + 1.0 - y[0];
    431:   PetscCall(MatSetValues(A, 2, &row[0], 2, &row[0], &value[0][0], INSERT_VALUES));
    432:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    433:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    434:   PetscCall(VecRestoreArrayRead(Y, &y));
    435:   PetscFunctionReturn(PETSC_SUCCESS);
    436: }

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

    440: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    441: {
    442:   PetscScalar       *f;
    443:   const PetscScalar *y;

    445:   PetscFunctionBeginUser;
    446:   PetscCall(VecGetArrayRead(Y, &y));
    447:   PetscCall(VecGetArray(F, &f));
    448:   f[0] = -y[0] + y[1];
    449:   f[1] = y[0] - 2.0 * y[1] + y[2];
    450:   f[2] = y[1] - y[2];
    451:   PetscCall(VecRestoreArrayRead(Y, &y));
    452:   PetscCall(VecRestoreArray(F, &f));
    453:   PetscFunctionReturn(PETSC_SUCCESS);
    454: }

    456: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    457: {
    458:   PetscScalar       *f;
    459:   const PetscScalar *y;

    461:   PetscFunctionBeginUser;
    462:   PetscCall(VecGetArrayRead(Y, &y));
    463:   PetscCall(VecGetArray(F, &f));
    464:   f[0] = -y[0] + y[1];
    465:   f[1] = y[0] - 2.0 * y[1] + y[2];
    466:   f[2] = y[1] - y[2];
    467:   PetscCall(VecRestoreArrayRead(Y, &y));
    468:   PetscCall(VecRestoreArray(F, &f));
    469:   /* Left hand side = ydot - f(y) */
    470:   PetscCall(VecAYPX(F, -1.0, Ydot));
    471:   PetscFunctionReturn(PETSC_SUCCESS);
    472: }

    474: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    475: {
    476:   const PetscScalar *y;
    477:   PetscInt           row[3] = {0, 1, 2};
    478:   PetscScalar        value[3][3];

    480:   PetscFunctionBeginUser;
    481:   PetscCall(VecGetArrayRead(Y, &y));
    482:   value[0][0] = a + 1.0;
    483:   value[0][1] = -1.0;
    484:   value[0][2] = 0;
    485:   value[1][0] = -1.0;
    486:   value[1][1] = a + 2.0;
    487:   value[1][2] = -1.0;
    488:   value[2][0] = 0;
    489:   value[2][1] = -1.0;
    490:   value[2][2] = a + 1.0;
    491:   PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
    492:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    493:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    494:   PetscCall(VecRestoreArrayRead(Y, &y));
    495:   PetscFunctionReturn(PETSC_SUCCESS);
    496: }

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

    500: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    501: {
    502:   PetscScalar       *f;
    503:   const PetscScalar *y;

    505:   PetscFunctionBeginUser;
    506:   PetscCall(VecGetArrayRead(Y, &y));
    507:   PetscCall(VecGetArray(F, &f));
    508:   f[0] = -y[0];
    509:   f[1] = y[0] - y[1] * y[1];
    510:   f[2] = y[1] * y[1];
    511:   PetscCall(VecRestoreArrayRead(Y, &y));
    512:   PetscCall(VecRestoreArray(F, &f));
    513:   PetscFunctionReturn(PETSC_SUCCESS);
    514: }

    516: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    517: {
    518:   PetscScalar       *f;
    519:   const PetscScalar *y;

    521:   PetscFunctionBeginUser;
    522:   PetscCall(VecGetArrayRead(Y, &y));
    523:   PetscCall(VecGetArray(F, &f));
    524:   f[0] = -y[0];
    525:   f[1] = y[0] - y[1] * y[1];
    526:   f[2] = y[1] * y[1];
    527:   PetscCall(VecRestoreArrayRead(Y, &y));
    528:   PetscCall(VecRestoreArray(F, &f));
    529:   /* Left hand side = ydot - f(y) */
    530:   PetscCall(VecAYPX(F, -1.0, Ydot));
    531:   PetscFunctionReturn(PETSC_SUCCESS);
    532: }

    534: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    535: {
    536:   const PetscScalar *y;
    537:   PetscInt           row[3] = {0, 1, 2};
    538:   PetscScalar        value[3][3];

    540:   PetscFunctionBeginUser;
    541:   PetscCall(VecGetArrayRead(Y, &y));
    542:   value[0][0] = a + 1.0;
    543:   value[0][1] = 0;
    544:   value[0][2] = 0;
    545:   value[1][0] = -1.0;
    546:   value[1][1] = a + 2.0 * y[1];
    547:   value[1][2] = 0;
    548:   value[2][0] = 0;
    549:   value[2][1] = -2.0 * y[1];
    550:   value[2][2] = a;
    551:   PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
    552:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    553:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    554:   PetscCall(VecRestoreArrayRead(Y, &y));
    555:   PetscFunctionReturn(PETSC_SUCCESS);
    556: }

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

    560: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    561: {
    562:   PetscScalar       *f;
    563:   const PetscScalar *y;

    565:   PetscFunctionBeginUser;
    566:   PetscCall(VecGetArrayRead(Y, &y));
    567:   PetscCall(VecGetArray(F, &f));
    568:   f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
    569:   f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
    570:   f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
    571:   PetscCall(VecRestoreArrayRead(Y, &y));
    572:   PetscCall(VecRestoreArray(F, &f));
    573:   PetscFunctionReturn(PETSC_SUCCESS);
    574: }

    576: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    577: {
    578:   PetscScalar       *f;
    579:   const PetscScalar *y;

    581:   PetscFunctionBeginUser;
    582:   PetscCall(VecGetArrayRead(Y, &y));
    583:   PetscCall(VecGetArray(F, &f));
    584:   f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
    585:   f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
    586:   f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
    587:   PetscCall(VecRestoreArrayRead(Y, &y));
    588:   PetscCall(VecRestoreArray(F, &f));
    589:   /* Left hand side = ydot - f(y) */
    590:   PetscCall(VecAYPX(F, -1.0, Ydot));
    591:   PetscFunctionReturn(PETSC_SUCCESS);
    592: }

    594: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    595: {
    596:   const PetscScalar *y;
    597:   PetscInt           row[3] = {0, 1, 2};
    598:   PetscScalar        value[3][3], fac, fac2;

    600:   PetscFunctionBeginUser;
    601:   PetscCall(VecGetArrayRead(Y, &y));
    602:   fac         = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -1.5);
    603:   fac2        = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -0.5);
    604:   value[0][0] = a + (y[1] * y[1] * y[2]) * fac;
    605:   value[0][1] = 1.0 - (y[0] * y[1] * y[2]) * fac;
    606:   value[0][2] = y[0] * fac2;
    607:   value[1][0] = -1.0 - y[0] * y[1] * y[2] * fac;
    608:   value[1][1] = a + y[0] * y[0] * y[2] * fac;
    609:   value[1][2] = y[1] * fac2;
    610:   value[2][0] = -y[1] * y[1] * fac;
    611:   value[2][1] = y[0] * y[1] * fac;
    612:   value[2][2] = a;
    613:   PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
    614:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    615:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    616:   PetscCall(VecRestoreArrayRead(Y, &y));
    617:   PetscFunctionReturn(PETSC_SUCCESS);
    618: }

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

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

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

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

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

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

    662:   PetscFunctionBeginUser;
    663:   PetscCall(VecGetArrayRead(Y, &y));
    664:   value[0][0] = a;
    665:   value[0][1] = -y[2];
    666:   value[0][2] = -y[1];
    667:   value[1][0] = y[2];
    668:   value[1][1] = a;
    669:   value[1][2] = y[0];
    670:   value[2][0] = 0.51 * y[1];
    671:   value[2][1] = 0.51 * y[0];
    672:   value[2][2] = a;
    673:   PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
    674:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    675:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    676:   PetscCall(VecRestoreArrayRead(Y, &y));
    677:   PetscFunctionReturn(PETSC_SUCCESS);
    678: }

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

    682: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    683: {
    684:   PetscScalar       *f;
    685:   const PetscScalar *y;

    687:   PetscFunctionBeginUser;
    688:   PetscCall(VecGetArrayRead(Y, &y));
    689:   PetscCall(VecGetArray(F, &f));
    690:   f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
    691:   f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
    692:   f[2] = 2. * t * y[3];
    693:   f[3] = -2. * t * PetscLogScalar(y[0]);
    694:   PetscCall(VecRestoreArrayRead(Y, &y));
    695:   PetscCall(VecRestoreArray(F, &f));
    696:   PetscFunctionReturn(PETSC_SUCCESS);
    697: }

    699: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
    700: {
    701:   const PetscScalar *y;
    702:   PetscInt           row[4] = {0, 1, 2, 3};
    703:   PetscScalar        value[4][4];
    704:   PetscScalar        m1, m2;

    706:   PetscFunctionBeginUser;
    707:   PetscCall(VecGetArrayRead(Y, &y));
    708:   m1          = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.));
    709:   m2          = 2. * t * PetscPowScalar(y[1], 1. / 5.);
    710:   value[0][0] = 0.;
    711:   value[0][1] = m1;
    712:   value[0][2] = 0.;
    713:   value[0][3] = m2;
    714:   m1          = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
    715:   m2          = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.));
    716:   value[1][0] = 0.;
    717:   value[1][1] = 0.;
    718:   value[1][2] = m1;
    719:   value[1][3] = m2;
    720:   value[2][0] = 0.;
    721:   value[2][1] = 0.;
    722:   value[2][2] = 0.;
    723:   value[2][3] = 2 * t;
    724:   value[3][0] = -2. * t / y[0];
    725:   value[3][1] = 0.;
    726:   value[3][2] = 0.;
    727:   value[3][3] = 0.;
    728:   PetscCall(MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES));
    729:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    730:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    731:   PetscCall(VecRestoreArrayRead(Y, &y));
    732:   PetscFunctionReturn(PETSC_SUCCESS);
    733: }

    735: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    736: {
    737:   PetscScalar       *f;
    738:   const PetscScalar *y;

    740:   PetscFunctionBeginUser;
    741:   PetscCall(VecGetArrayRead(Y, &y));
    742:   PetscCall(VecGetArray(F, &f));
    743:   f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
    744:   f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
    745:   f[2] = 2. * t * y[3];
    746:   f[3] = -2. * t * PetscLogScalar(y[0]);
    747:   PetscCall(VecRestoreArrayRead(Y, &y));
    748:   PetscCall(VecRestoreArray(F, &f));
    749:   /* Left hand side = ydot - f(y) */
    750:   PetscCall(VecAYPX(F, -1.0, Ydot));
    751:   PetscFunctionReturn(PETSC_SUCCESS);
    752: }

    754: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    755: {
    756:   const PetscScalar *y;
    757:   PetscInt           row[4] = {0, 1, 2, 3};
    758:   PetscScalar        value[4][4];
    759:   PetscScalar        m1, m2;

    761:   PetscFunctionBeginUser;
    762:   PetscCall(VecGetArrayRead(Y, &y));
    763:   m1          = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.));
    764:   m2          = 2. * t * PetscPowScalar(y[1], 1. / 5.);
    765:   value[0][0] = a;
    766:   value[0][1] = m1;
    767:   value[0][2] = 0.;
    768:   value[0][3] = m2;
    769:   m1          = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
    770:   m2          = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.));
    771:   value[1][0] = 0.;
    772:   value[1][1] = a;
    773:   value[1][2] = m1;
    774:   value[1][3] = m2;
    775:   value[2][0] = 0.;
    776:   value[2][1] = 0.;
    777:   value[2][2] = a;
    778:   value[2][3] = 2 * t;
    779:   value[3][0] = -2. * t / y[0];
    780:   value[3][1] = 0.;
    781:   value[3][2] = 0.;
    782:   value[3][3] = a;
    783:   PetscCall(MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES));
    784:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    785:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    786:   PetscCall(VecRestoreArrayRead(Y, &y));
    787:   PetscFunctionReturn(PETSC_SUCCESS);
    788: }

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

    792: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    793: {
    794:   PetscScalar       *f;
    795:   const PetscScalar *y;
    796:   PetscInt           N, i;

    798:   PetscFunctionBeginUser;
    799:   PetscCall(VecGetSize(Y, &N));
    800:   PetscCall(VecGetArrayRead(Y, &y));
    801:   PetscCall(VecGetArray(F, &f));
    802:   f[0] = -y[0];
    803:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - y[i];
    804:   f[N - 1] = y[N - 2];
    805:   PetscCall(VecRestoreArrayRead(Y, &y));
    806:   PetscCall(VecRestoreArray(F, &f));
    807:   PetscFunctionReturn(PETSC_SUCCESS);
    808: }

    810: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    811: {
    812:   PetscScalar       *f;
    813:   const PetscScalar *y;
    814:   PetscInt           N, i;

    816:   PetscFunctionBeginUser;
    817:   PetscCall(VecGetSize(Y, &N));
    818:   PetscCall(VecGetArrayRead(Y, &y));
    819:   PetscCall(VecGetArray(F, &f));
    820:   f[0] = -y[0];
    821:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - y[i];
    822:   f[N - 1] = y[N - 2];
    823:   PetscCall(VecRestoreArrayRead(Y, &y));
    824:   PetscCall(VecRestoreArray(F, &f));
    825:   /* Left hand side = ydot - f(y) */
    826:   PetscCall(VecAYPX(F, -1.0, Ydot));
    827:   PetscFunctionReturn(PETSC_SUCCESS);
    828: }

    830: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    831: {
    832:   const PetscScalar *y;
    833:   PetscInt           N, i, col[2];
    834:   PetscScalar        value[2];

    836:   PetscFunctionBeginUser;
    837:   PetscCall(VecGetSize(Y, &N));
    838:   PetscCall(VecGetArrayRead(Y, &y));
    839:   i        = 0;
    840:   value[0] = a + 1;
    841:   col[0]   = 0;
    842:   value[1] = 0;
    843:   col[1]   = 1;
    844:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
    845:   for (i = 0; i < N; i++) {
    846:     value[0] = -1;
    847:     col[0]   = i - 1;
    848:     value[1] = a + 1;
    849:     col[1]   = i;
    850:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
    851:   }
    852:   i        = N - 1;
    853:   value[0] = -1;
    854:   col[0]   = N - 2;
    855:   value[1] = a;
    856:   col[1]   = N - 1;
    857:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
    858:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    859:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    860:   PetscCall(VecRestoreArrayRead(Y, &y));
    861:   PetscFunctionReturn(PETSC_SUCCESS);
    862: }

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

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

    872:   PetscFunctionBeginUser;
    873:   PetscCall(VecGetSize(Y, &N));
    874:   PetscCall(VecGetArrayRead(Y, &y));
    875:   PetscCall(VecGetArray(F, &f));
    876:   f[0] = -y[0];
    877:   for (i = 1; i < N - 1; i++) f[i] = (PetscReal)i * y[i - 1] - (PetscReal)(i + 1) * y[i];
    878:   f[N - 1] = (PetscReal)(N - 1) * y[N - 2];
    879:   PetscCall(VecRestoreArrayRead(Y, &y));
    880:   PetscCall(VecRestoreArray(F, &f));
    881:   PetscFunctionReturn(PETSC_SUCCESS);
    882: }

    884: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    885: {
    886:   PetscScalar       *f;
    887:   const PetscScalar *y;
    888:   PetscInt           N, i;

    890:   PetscFunctionBeginUser;
    891:   PetscCall(VecGetSize(Y, &N));
    892:   PetscCall(VecGetArrayRead(Y, &y));
    893:   PetscCall(VecGetArray(F, &f));
    894:   f[0] = -y[0];
    895:   for (i = 1; i < N - 1; i++) f[i] = (PetscReal)i * y[i - 1] - (PetscReal)(i + 1) * y[i];
    896:   f[N - 1] = (PetscReal)(N - 1) * y[N - 2];
    897:   PetscCall(VecRestoreArrayRead(Y, &y));
    898:   PetscCall(VecRestoreArray(F, &f));
    899:   /* Left hand side = ydot - f(y) */
    900:   PetscCall(VecAYPX(F, -1.0, Ydot));
    901:   PetscFunctionReturn(PETSC_SUCCESS);
    902: }

    904: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    905: {
    906:   const PetscScalar *y;
    907:   PetscInt           N, i, col[2];
    908:   PetscScalar        value[2];

    910:   PetscFunctionBeginUser;
    911:   PetscCall(VecGetSize(Y, &N));
    912:   PetscCall(VecGetArrayRead(Y, &y));
    913:   i        = 0;
    914:   value[0] = a + 1;
    915:   col[0]   = 0;
    916:   value[1] = 0;
    917:   col[1]   = 1;
    918:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
    919:   for (i = 0; i < N; i++) {
    920:     value[0] = -(PetscReal)i;
    921:     col[0]   = i - 1;
    922:     value[1] = a + (PetscReal)(i + 1);
    923:     col[1]   = i;
    924:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
    925:   }
    926:   i        = N - 1;
    927:   value[0] = -(PetscReal)(N - 1);
    928:   col[0]   = N - 2;
    929:   value[1] = a;
    930:   col[1]   = N - 1;
    931:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
    932:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    933:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    934:   PetscCall(VecRestoreArrayRead(Y, &y));
    935:   PetscFunctionReturn(PETSC_SUCCESS);
    936: }

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

    940: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
    941: {
    942:   PetscScalar       *f;
    943:   const PetscScalar *y;
    944:   PetscInt           N, i;

    946:   PetscFunctionBeginUser;
    947:   PetscCall(VecGetSize(Y, &N));
    948:   PetscCall(VecGetArrayRead(Y, &y));
    949:   PetscCall(VecGetArray(F, &f));
    950:   f[0] = -2.0 * y[0] + y[1];
    951:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - 2.0 * y[i] + y[i + 1];
    952:   f[N - 1] = y[N - 2] - 2.0 * y[N - 1];
    953:   PetscCall(VecRestoreArrayRead(Y, &y));
    954:   PetscCall(VecRestoreArray(F, &f));
    955:   PetscFunctionReturn(PETSC_SUCCESS);
    956: }

    958: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
    959: {
    960:   PetscScalar       *f;
    961:   const PetscScalar *y;
    962:   PetscInt           N, i;

    964:   PetscFunctionBeginUser;
    965:   PetscCall(VecGetSize(Y, &N));
    966:   PetscCall(VecGetArrayRead(Y, &y));
    967:   PetscCall(VecGetArray(F, &f));
    968:   f[0] = -2.0 * y[0] + y[1];
    969:   for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - 2.0 * y[i] + y[i + 1];
    970:   f[N - 1] = y[N - 2] - 2.0 * y[N - 1];
    971:   PetscCall(VecRestoreArrayRead(Y, &y));
    972:   PetscCall(VecRestoreArray(F, &f));
    973:   /* Left hand side = ydot - f(y) */
    974:   PetscCall(VecAYPX(F, -1.0, Ydot));
    975:   PetscFunctionReturn(PETSC_SUCCESS);
    976: }

    978: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
    979: {
    980:   const PetscScalar *y;
    981:   PetscScalar        value[3];
    982:   PetscInt           N, i, col[3];

    984:   PetscFunctionBeginUser;
    985:   PetscCall(VecGetSize(Y, &N));
    986:   PetscCall(VecGetArrayRead(Y, &y));
    987:   for (i = 0; i < N; i++) {
    988:     if (i == 0) {
    989:       value[0] = a + 2;
    990:       col[0]   = i;
    991:       value[1] = -1;
    992:       col[1]   = i + 1;
    993:       value[2] = 0;
    994:       col[2]   = i + 2;
    995:     } else if (i == N - 1) {
    996:       value[0] = 0;
    997:       col[0]   = i - 2;
    998:       value[1] = -1;
    999:       col[1]   = i - 1;
    1000:       value[2] = a + 2;
    1001:       col[2]   = i;
    1002:     } else {
    1003:       value[0] = -1;
    1004:       col[0]   = i - 1;
    1005:       value[1] = a + 2;
    1006:       col[1]   = i;
    1007:       value[2] = -1;
    1008:       col[2]   = i + 1;
    1009:     }
    1010:     PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
    1011:   }
    1012:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    1013:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    1014:   PetscCall(VecRestoreArrayRead(Y, &y));
    1015:   PetscFunctionReturn(PETSC_SUCCESS);
    1016: }

    1018: /***************************************************************************/

    1020: /* Sets the initial solution for the IVP and sets up the function pointers*/
    1021: PetscErrorCode Initialize(Vec Y, void *s)
    1022: {
    1023:   char        *p = (char *)s;
    1024:   PetscScalar *y;
    1025:   PetscReal    t0;
    1026:   PetscInt     N;
    1027:   PetscBool    flg;

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

    1128: /* Calculates the exact solution to problems that have one */
    1129: PetscErrorCode ExactSolution(Vec Y, void *s, PetscReal t, PetscBool *flag)
    1130: {
    1131:   char        *p = (char *)s;
    1132:   PetscScalar *y;

    1134:   PetscFunctionBeginUser;
    1135:   if (!strcmp(p, "hull1972a1")) {
    1136:     PetscCall(VecGetArray(Y, &y));
    1137:     y[0]  = PetscExpReal(-t);
    1138:     *flag = PETSC_TRUE;
    1139:     PetscCall(VecRestoreArray(Y, &y));
    1140:   } else if (!strcmp(p, "hull1972a2")) {
    1141:     PetscCall(VecGetArray(Y, &y));
    1142:     y[0]  = 1.0 / PetscSqrtReal(t + 1);
    1143:     *flag = PETSC_TRUE;
    1144:     PetscCall(VecRestoreArray(Y, &y));
    1145:   } else if (!strcmp(p, "hull1972a3")) {
    1146:     PetscCall(VecGetArray(Y, &y));
    1147:     y[0]  = PetscExpReal(PetscSinReal(t));
    1148:     *flag = PETSC_TRUE;
    1149:     PetscCall(VecRestoreArray(Y, &y));
    1150:   } else if (!strcmp(p, "hull1972a4")) {
    1151:     PetscCall(VecGetArray(Y, &y));
    1152:     y[0]  = 20.0 / (1 + 19.0 * PetscExpReal(-t / 4.0));
    1153:     *flag = PETSC_TRUE;
    1154:     PetscCall(VecRestoreArray(Y, &y));
    1155:   } else if (!strcmp(p, "kulik2013i")) {
    1156:     PetscCall(VecGetArray(Y, &y));
    1157:     y[0]  = PetscExpReal(PetscSinReal(t * t));
    1158:     y[1]  = PetscExpReal(5. * PetscSinReal(t * t));
    1159:     y[2]  = PetscSinReal(t * t) + 1.0;
    1160:     y[3]  = PetscCosReal(t * t);
    1161:     *flag = PETSC_TRUE;
    1162:     PetscCall(VecRestoreArray(Y, &y));
    1163:   } else {
    1164:     PetscCall(VecSet(Y, 0));
    1165:     *flag = PETSC_FALSE;
    1166:   }
    1167:   PetscFunctionReturn(PETSC_SUCCESS);
    1168: }

    1170: /* Solves the specified ODE and computes the error if exact solution is available */
    1171: PetscErrorCode SolveODE(char *ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
    1172: {
    1173:   TS        ts;          /* time-integrator                        */
    1174:   Vec       Y;           /* Solution vector                        */
    1175:   Vec       Yex;         /* Exact solution                         */
    1176:   PetscInt  N;           /* Size of the system of equations        */
    1177:   TSType    time_scheme; /* Type of time-integration scheme        */
    1178:   Mat       Jac = NULL;  /* Jacobian matrix                        */
    1179:   Vec       Yerr;        /* Auxiliary solution vector              */
    1180:   PetscReal err_norm;    /* Estimated error norm                   */
    1181:   PetscReal final_time;  /* Actual final time from the integrator  */

    1183:   PetscFunctionBeginUser;
    1184:   PetscCall(GetSize((const char *)&ptype[0], &N));
    1185:   PetscCheck(N >= 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Illegal problem specification.");
    1186:   PetscCall(VecCreate(PETSC_COMM_WORLD, &Y));
    1187:   PetscCall(VecSetSizes(Y, N, PETSC_DECIDE));
    1188:   PetscCall(VecSetUp(Y));
    1189:   PetscCall(VecSet(Y, 0));

    1191:   /* Initialize the problem */
    1192:   PetscCall(Initialize(Y, &ptype[0]));

    1194:   /* Create and initialize the time-integrator                            */
    1195:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    1196:   /* Default time integration options                                     */
    1197:   PetscCall(TSSetType(ts, TSRK));
    1198:   PetscCall(TSSetMaxSteps(ts, maxiter));
    1199:   PetscCall(TSSetMaxTime(ts, tfinal));
    1200:   PetscCall(TSSetTimeStep(ts, dt));
    1201:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    1202:   /* Read command line options for time integration                       */
    1203:   PetscCall(TSSetFromOptions(ts));
    1204:   /* Set solution vector                                                  */
    1205:   PetscCall(TSSetSolution(ts, Y));
    1206:   /* Specify left/right-hand side functions                               */
    1207:   PetscCall(TSGetType(ts, &time_scheme));

    1209:   if ((!strcmp(time_scheme, TSEULER)) || (!strcmp(time_scheme, TSRK)) || (!strcmp(time_scheme, TSSSP) || (!strcmp(time_scheme, TSGLEE)))) {
    1210:     /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
    1211:     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &ptype[0]));
    1212:     PetscCall(MatCreate(PETSC_COMM_WORLD, &Jac));
    1213:     PetscCall(MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N));
    1214:     PetscCall(MatSetFromOptions(Jac));
    1215:     PetscCall(MatSetUp(Jac));
    1216:     PetscCall(TSSetRHSJacobian(ts, Jac, Jac, RHSJacobian, &ptype[0]));
    1217:   } else if ((!strcmp(time_scheme, TSTHETA)) || (!strcmp(time_scheme, TSBEULER)) || (!strcmp(time_scheme, TSCN)) || (!strcmp(time_scheme, TSALPHA)) || (!strcmp(time_scheme, TSARKIMEX))) {
    1218:     /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
    1219:     /* and its Jacobian function                                                 */
    1220:     PetscCall(TSSetIFunction(ts, NULL, IFunction, &ptype[0]));
    1221:     PetscCall(MatCreate(PETSC_COMM_WORLD, &Jac));
    1222:     PetscCall(MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N));
    1223:     PetscCall(MatSetFromOptions(Jac));
    1224:     PetscCall(MatSetUp(Jac));
    1225:     PetscCall(TSSetIJacobian(ts, Jac, Jac, IJacobian, &ptype[0]));
    1226:   }

    1228:   /* Solve */
    1229:   PetscCall(TSSolve(ts, Y));
    1230:   PetscCall(TSGetTime(ts, &final_time));

    1232:   /* Get the estimated error, if available */
    1233:   PetscCall(VecDuplicate(Y, &Yerr));
    1234:   PetscCall(VecZeroEntries(Yerr));
    1235:   PetscCall(TSGetTimeError(ts, 0, &Yerr));
    1236:   PetscCall(VecNorm(Yerr, NORM_2, &err_norm));
    1237:   PetscCall(VecDestroy(&Yerr));
    1238:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Estimated Error = %e.\n", (double)err_norm));

    1240:   /* Exact solution */
    1241:   PetscCall(VecDuplicate(Y, &Yex));
    1242:   if (PetscAbsReal(final_time - tfinal) > 2. * PETSC_MACHINE_EPSILON) PetscCall(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));
    1243:   PetscCall(ExactSolution(Yex, &ptype[0], final_time, exact_flag));

    1245:   /* Calculate Error */
    1246:   PetscCall(VecAYPX(Yex, -1.0, Y));
    1247:   PetscCall(VecNorm(Yex, NORM_2, error));
    1248:   *error = PetscSqrtReal(((*error) * (*error)) / N);

    1250:   /* Clean up and finalize */
    1251:   PetscCall(MatDestroy(&Jac));
    1252:   PetscCall(TSDestroy(&ts));
    1253:   PetscCall(VecDestroy(&Yex));
    1254:   PetscCall(VecDestroy(&Y));
    1255:   PetscFunctionReturn(PETSC_SUCCESS);
    1256: }

    1258: int main(int argc, char **argv)
    1259: {
    1260:   char        ptype[256] = "hull1972a1"; /* Problem specification                                */
    1261:   PetscInt    n_refine   = 1;            /* Number of refinement levels for convergence analysis */
    1262:   PetscReal   refine_fac = 2.0;          /* Refinement factor for dt                             */
    1263:   PetscReal   dt_initial = 0.01;         /* Initial default value of dt                          */
    1264:   PetscReal   dt;
    1265:   PetscReal   tfinal  = 20.0;   /* Final time for the time-integration                  */
    1266:   PetscInt    maxiter = 100000; /* Maximum number of time-integration iterations        */
    1267:   PetscReal  *error;            /* Array to store the errors for convergence analysis   */
    1268:   PetscMPIInt size;             /* No of processors                                     */
    1269:   PetscBool   flag;             /* Flag denoting availability of exact solution         */
    1270:   PetscInt    r, N;

    1272:   /* Initialize program */
    1273:   PetscFunctionBeginUser;
    1274:   PetscCall(GetSize(&ptype[0], &N));
    1275:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

    1277:   /* Check if running with only 1 proc */
    1278:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    1279:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    1281:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL);
    1282:   PetscCall(PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL));
    1283:   PetscCall(PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL));
    1284:   PetscCall(PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL));
    1285:   PetscCall(PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL));
    1286:   PetscCall(PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL));
    1287:   PetscOptionsEnd();

    1289:   PetscCall(PetscMalloc1(n_refine, &error));
    1290:   for (r = 0, dt = dt_initial; r < n_refine; r++) {
    1291:     error[r] = 0;
    1292:     if (r > 0) dt /= refine_fac;

    1294:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving ODE \"%s\" with dt %f, final time %f and system size %" PetscInt_FMT ".\n", ptype, (double)dt, (double)tfinal, N));
    1295:     PetscCall(SolveODE(&ptype[0], dt, tfinal, maxiter, &error[r], &flag));
    1296:     if (flag) {
    1297:       /* If exact solution available for the specified ODE */
    1298:       if (r > 0) {
    1299:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r - 1])) / (-PetscLogReal(refine_fac));
    1300:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error           = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate));
    1301:       } else {
    1302:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error           = %E.\n", (double)error[r]));
    1303:       }
    1304:     }
    1305:   }
    1306:   PetscCall(PetscFree(error));
    1307:   PetscCall(PetscFinalize());
    1308:   return 0;
    1309: }

    1311: /*TEST

    1313:     test:
    1314:       suffix: 2
    1315:       args: -ts_type glee -final_time 5 -ts_adapt_type none
    1316:       timeoutfactor: 3
    1317:       requires: !single

    1319:     test:
    1320:       suffix: 3
    1321:       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
    1322:       timeoutfactor: 3
    1323:       requires: !single

    1325:     test:
    1326:       suffix: 4
    1327:       # the test is so sensitive that I could not even replace VecMAXPY with a sequence of VECAXPY, so I just disable this GEMV optimization
    1328:       args: -vec_maxpy_use_gemv 0
    1329:       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
    1330:       timeoutfactor: 3
    1331:       requires: !single !__float128

    1333: TEST*/