Actual source code: biharmonic.c

      1: static char help[] = "Solves biharmonic equation in 1d.\n";

      3: /*
      4:   Solves the equation

      6:     u_t = - kappa  \Delta \Delta u
      7:     Periodic boundary conditions

      9: Evolve the biharmonic heat equation:
     10: ---------------
     11: ./biharmonic -ts_monitor -snes_monitor   -pc_type lu  -draw_pause .1 -snes_converged_reason  -draw_pause -2   -ts_type cn  -da_refine 5 -mymonitor

     13: Evolve with the restriction that -1 <= u <= 1; i.e. as a variational inequality
     14: ---------------
     15: ./biharmonic -ts_monitor -snes_monitor   -pc_type lu  -draw_pause .1 -snes_converged_reason  -draw_pause -2   -ts_type cn   -da_refine 5  -mymonitor

     17:    u_t =  kappa \Delta \Delta u +   6.*u*(u_x)^2 + (3*u^2 - 12) \Delta u
     18:     -1 <= u <= 1
     19:     Periodic boundary conditions

     21: Evolve the Cahn-Hillard equations: double well Initial hump shrinks then grows
     22: ---------------
     23: ./biharmonic -ts_monitor -snes_monitor   -pc_type lu  -draw_pause .1 -snes_converged_reason   -draw_pause -2   -ts_type cn    -da_refine 6   -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -ts_monitor_draw_solution --mymonitor

     25: Initial hump neither shrinks nor grows when degenerate (otherwise similar solution)

     27: ./biharmonic -ts_monitor -snes_monitor   -pc_type lu  -draw_pause .1 -snes_converged_reason   -draw_pause -2   -ts_type cn    -da_refine 6   -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -degenerate -ts_monitor_draw_solution --mymonitor

     29: ./biharmonic -ts_monitor -snes_monitor   -pc_type lu  -draw_pause .1 -snes_converged_reason   -draw_pause -2   -ts_type cn    -da_refine 6   -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -snes_vi_ignore_function_sign -ts_monitor_draw_solution --mymonitor

     31: Evolve the Cahn-Hillard equations: double obstacle
     32: ---------------
     33: ./biharmonic -ts_monitor -snes_monitor  -pc_type lu  -draw_pause .1 -snes_converged_reason   -draw_pause -2   -ts_type cn    -da_refine 5   -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 2 -snes_linesearch_monitor    -ts_monitor_draw_solution --mymonitor

     35: Evolve the Cahn-Hillard equations: logarithmic + double well (never shrinks and then grows)
     36: ---------------
     37: ./biharmonic -ts_monitor -snes_monitor  -pc_type lu  --snes_converged_reason  -draw_pause -2   -ts_type cn    -da_refine 5   -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001    -ts_monitor_draw_solution --ts_max_time 1. -mymonitor

     39: ./biharmonic -ts_monitor -snes_monitor  -pc_type lu  --snes_converged_reason  -draw_pause -2   -ts_type cn    -da_refine 5   -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001    -ts_monitor_draw_solution --ts_max_time 1. -degenerate -mymonitor

     41: Evolve the Cahn-Hillard equations: logarithmic +  double obstacle (never shrinks, never grows)
     42: ---------------
     43: ./biharmonic -ts_monitor -snes_monitor  -pc_type lu  --snes_converged_reason  -draw_pause -2   -ts_type cn    -da_refine 5   -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 4 -snes_linesearch_monitor -theta .00000001   -ts_monitor_draw_solution --mymonitor

     45: */
     46: #include <petscdm.h>
     47: #include <petscdmda.h>
     48: #include <petscts.h>
     49: #include <petscdraw.h>

     51: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), MyDestroy(void **), FormJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
     52: typedef struct {
     53:   PetscBool           cahnhillard;
     54:   PetscBool           degenerate;
     55:   PetscReal           kappa;
     56:   PetscInt            energy;
     57:   PetscReal           tol;
     58:   PetscReal           theta, theta_c;
     59:   PetscInt            truncation;
     60:   PetscBool           netforce;
     61:   PetscDrawViewPorts *ports;
     62: } UserCtx;

     64: int main(int argc, char **argv)
     65: {
     66:   TS        ts;   /* nonlinear solver */
     67:   Vec       x, r; /* solution, residual vectors */
     68:   Mat       J;    /* Jacobian matrix */
     69:   PetscInt  steps, Mx;
     70:   DM        da;
     71:   PetscReal dt;
     72:   PetscBool mymonitor;
     73:   UserCtx   ctx;

     75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:      Initialize program
     77:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     78:   PetscFunctionBeginUser;
     79:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     80:   ctx.kappa = 1.0;
     81:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL));
     82:   ctx.degenerate = PETSC_FALSE;
     83:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-degenerate", &ctx.degenerate, NULL));
     84:   ctx.cahnhillard = PETSC_FALSE;
     85:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-cahn-hillard", &ctx.cahnhillard, NULL));
     86:   ctx.netforce = PETSC_FALSE;
     87:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-netforce", &ctx.netforce, NULL));
     88:   ctx.energy = 1;
     89:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-energy", &ctx.energy, NULL));
     90:   ctx.tol = 1.0e-8;
     91:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &ctx.tol, NULL));
     92:   ctx.theta   = .001;
     93:   ctx.theta_c = 1.0;
     94:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-theta", &ctx.theta, NULL));
     95:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-theta_c", &ctx.theta_c, NULL));
     96:   ctx.truncation = 1;
     97:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-truncation", &ctx.truncation, NULL));
     98:   PetscCall(PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor));

    100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101:      Create distributed array (DMDA) to manage parallel grid and vectors
    102:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    103:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da));
    104:   PetscCall(DMSetFromOptions(da));
    105:   PetscCall(DMSetUp(da));
    106:   PetscCall(DMDASetFieldName(da, 0, "Biharmonic heat equation: u"));
    107:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    108:   dt = 1.0 / (10. * ctx.kappa * Mx * Mx * Mx * Mx);

    110:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    111:      Extract global vectors from DMDA; then duplicate for remaining
    112:      vectors that are the same types
    113:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    114:   PetscCall(DMCreateGlobalVector(da, &x));
    115:   PetscCall(VecDuplicate(x, &r));

    117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    118:      Create timestepping solver context
    119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    120:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    121:   PetscCall(TSSetDM(ts, da));
    122:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    123:   PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &ctx));
    124:   PetscCall(DMSetMatType(da, MATAIJ));
    125:   PetscCall(DMCreateMatrix(da, &J));
    126:   PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &ctx));
    127:   PetscCall(TSSetMaxTime(ts, .02));
    128:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE));

    130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    131:      Create matrix data structure; set Jacobian evaluation routine

    133:      Set Jacobian matrix data structure and default Jacobian evaluation
    134:      routine. User can override with:
    135:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
    136:                 (unless user explicitly sets preconditioner)
    137:      -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
    138:                          but use matrix-free approx for Jacobian-vector
    139:                          products within Newton-Krylov method

    141:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    143:      Customize nonlinear solver
    144:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    145:   PetscCall(TSSetType(ts, TSCN));

    147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    148:      Set initial conditions
    149:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    150:   PetscCall(FormInitialSolution(da, x));
    151:   PetscCall(TSSetTimeStep(ts, dt));
    152:   PetscCall(TSSetSolution(ts, x));

    154:   if (mymonitor) {
    155:     ctx.ports = NULL;
    156:     PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
    157:   }

    159:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    160:      Set runtime options
    161:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    162:   PetscCall(TSSetFromOptions(ts));

    164:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    165:      Solve nonlinear system
    166:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    167:   PetscCall(TSSolve(ts, x));
    168:   PetscCall(TSGetStepNumber(ts, &steps));
    169:   PetscCall(VecView(x, PETSC_VIEWER_BINARY_WORLD));

    171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172:      Free work space.  All PETSc objects should be destroyed when they
    173:      are no longer needed.
    174:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    175:   PetscCall(MatDestroy(&J));
    176:   PetscCall(VecDestroy(&x));
    177:   PetscCall(VecDestroy(&r));
    178:   PetscCall(TSDestroy(&ts));
    179:   PetscCall(DMDestroy(&da));

    181:   PetscCall(PetscFinalize());
    182:   return 0;
    183: }
    184: /* ------------------------------------------------------------------- */
    185: /*
    186:    FormFunction - Evaluates nonlinear function, F(x).

    188:    Input Parameters:
    189: .  ts - the TS context
    190: .  X - input vector
    191: .  ptr - optional user-defined context, as set by SNESSetFunction()

    193:    Output Parameter:
    194: .  F - function vector
    195:  */
    196: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
    197: {
    198:   DM           da;
    199:   PetscInt     i, Mx, xs, xm;
    200:   PetscReal    hx, sx;
    201:   PetscScalar *x, *f, c, r, l;
    202:   Vec          localX;
    203:   UserCtx     *ctx = (UserCtx *)ptr;
    204:   PetscReal    tol = ctx->tol, theta = ctx->theta, theta_c = ctx->theta_c, a, b; /* a and b are used in the cubic truncation of the log function */

    206:   PetscFunctionBegin;
    207:   PetscCall(TSGetDM(ts, &da));
    208:   PetscCall(DMGetLocalVector(da, &localX));
    209:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

    211:   hx = 1.0 / (PetscReal)Mx;
    212:   sx = 1.0 / (hx * hx);

    214:   /*
    215:      Scatter ghost points to local vector,using the 2-step process
    216:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    217:      By placing code between these two statements, computations can be
    218:      done while messages are in transition.
    219:   */
    220:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    221:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

    223:   /*
    224:      Get pointers to vector data
    225:   */
    226:   PetscCall(DMDAVecGetArrayRead(da, localX, &x));
    227:   PetscCall(DMDAVecGetArray(da, F, &f));

    229:   /*
    230:      Get local grid boundaries
    231:   */
    232:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    234:   /*
    235:      Compute function over the locally owned part of the grid
    236:   */
    237:   for (i = xs; i < xs + xm; i++) {
    238:     if (ctx->degenerate) {
    239:       c = (1. - x[i] * x[i]) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    240:       r = (1. - x[i + 1] * x[i + 1]) * (x[i] + x[i + 2] - 2.0 * x[i + 1]) * sx;
    241:       l = (1. - x[i - 1] * x[i - 1]) * (x[i - 2] + x[i] - 2.0 * x[i - 1]) * sx;
    242:     } else {
    243:       c = (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    244:       r = (x[i] + x[i + 2] - 2.0 * x[i + 1]) * sx;
    245:       l = (x[i - 2] + x[i] - 2.0 * x[i - 1]) * sx;
    246:     }
    247:     f[i] = -ctx->kappa * (l + r - 2.0 * c) * sx;
    248:     if (ctx->cahnhillard) {
    249:       switch (ctx->energy) {
    250:       case 1: /*  double well */
    251:         f[i] += 6. * .25 * x[i] * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (3. * x[i] * x[i] - 1.) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    252:         break;
    253:       case 2: /* double obstacle */
    254:         f[i] += -(x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    255:         break;
    256:       case 3: /* logarithmic + double well */
    257:         f[i] += 6. * .25 * x[i] * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (3. * x[i] * x[i] - 1.) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    258:         if (ctx->truncation == 2) { /* log function with approximated with a quadratic polynomial outside -1.0+2*tol, 1.0-2*tol */
    259:           if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    260:           else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    261:           else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    262:         } else { /* log function is approximated with a cubic polynomial outside -1.0+2*tol, 1.0-2*tol */
    263:           a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
    264:           b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
    265:           if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += -1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (-1.0 * a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    266:           else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += 1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    267:           else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    268:         }
    269:         break;
    270:       case 4: /* logarithmic + double obstacle */
    271:         f[i] += -theta_c * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    272:         if (ctx->truncation == 2) { /* quadratic */
    273:           if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    274:           else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    275:           else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    276:         } else { /* cubic */
    277:           a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
    278:           b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
    279:           if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += -1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (-1.0 * a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    280:           else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += 1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    281:           else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    282:         }
    283:         break;
    284:       }
    285:     }
    286:   }

    288:   /*
    289:      Restore vectors
    290:   */
    291:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
    292:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    293:   PetscCall(DMRestoreLocalVector(da, &localX));
    294:   PetscFunctionReturn(PETSC_SUCCESS);
    295: }

    297: /* ------------------------------------------------------------------- */
    298: /*
    299:    FormJacobian - Evaluates nonlinear function's Jacobian

    301: */
    302: PetscErrorCode FormJacobian(TS ts, PetscReal ftime, Vec X, Mat A, Mat B, void *ptr)
    303: {
    304:   DM           da;
    305:   PetscInt     i, Mx, xs, xm;
    306:   MatStencil   row, cols[5];
    307:   PetscReal    hx, sx;
    308:   PetscScalar *x, vals[5];
    309:   Vec          localX;
    310:   UserCtx     *ctx = (UserCtx *)ptr;

    312:   PetscFunctionBegin;
    313:   PetscCall(TSGetDM(ts, &da));
    314:   PetscCall(DMGetLocalVector(da, &localX));
    315:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

    317:   hx = 1.0 / (PetscReal)Mx;
    318:   sx = 1.0 / (hx * hx);

    320:   /*
    321:      Scatter ghost points to local vector,using the 2-step process
    322:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    323:      By placing code between these two statements, computations can be
    324:      done while messages are in transition.
    325:   */
    326:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    327:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

    329:   /*
    330:      Get pointers to vector data
    331:   */
    332:   PetscCall(DMDAVecGetArrayRead(da, localX, &x));

    334:   /*
    335:      Get local grid boundaries
    336:   */
    337:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    339:   /*
    340:      Compute function over the locally owned part of the grid
    341:   */
    342:   for (i = xs; i < xs + xm; i++) {
    343:     row.i = i;
    344:     if (ctx->degenerate) {
    345:       /*PetscScalar c,r,l;
    346:       c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
    347:       r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx;
    348:       l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx; */
    349:     } else {
    350:       cols[0].i = i - 2;
    351:       vals[0]   = -ctx->kappa * sx * sx;
    352:       cols[1].i = i - 1;
    353:       vals[1]   = 4.0 * ctx->kappa * sx * sx;
    354:       cols[2].i = i;
    355:       vals[2]   = -6.0 * ctx->kappa * sx * sx;
    356:       cols[3].i = i + 1;
    357:       vals[3]   = 4.0 * ctx->kappa * sx * sx;
    358:       cols[4].i = i + 2;
    359:       vals[4]   = -ctx->kappa * sx * sx;
    360:     }
    361:     PetscCall(MatSetValuesStencil(B, 1, &row, 5, cols, vals, INSERT_VALUES));

    363:     if (ctx->cahnhillard) {
    364:       switch (ctx->energy) {
    365:       case 1: /* double well */
    366:         /*  f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; */
    367:         break;
    368:       case 2: /* double obstacle */
    369:         /*        f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx; */
    370:         break;
    371:       case 3: /* logarithmic + double well */
    372:         break;
    373:       case 4: /* logarithmic + double obstacle */
    374:         break;
    375:       }
    376:     }
    377:   }

    379:   /*
    380:      Restore vectors
    381:   */
    382:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
    383:   PetscCall(DMRestoreLocalVector(da, &localX));
    384:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    385:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    386:   if (A != B) {
    387:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    388:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    389:   }
    390:   PetscFunctionReturn(PETSC_SUCCESS);
    391: }
    392: /* ------------------------------------------------------------------- */
    393: PetscErrorCode FormInitialSolution(DM da, Vec U)
    394: {
    395:   PetscInt           i, xs, xm, Mx, N, scale;
    396:   PetscScalar       *u;
    397:   PetscReal          r, hx, x;
    398:   const PetscScalar *f;
    399:   Vec                finesolution;
    400:   PetscViewer        viewer;

    402:   PetscFunctionBegin;
    403:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

    405:   hx = 1.0 / (PetscReal)Mx;

    407:   /*
    408:      Get pointers to vector data
    409:   */
    410:   PetscCall(DMDAVecGetArray(da, U, &u));

    412:   /*
    413:      Get local grid boundaries
    414:   */
    415:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    417:   /*
    418:       Seee heat.c for how to generate InitialSolution.heat
    419:   */
    420:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer));
    421:   PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution));
    422:   PetscCall(VecLoad(finesolution, viewer));
    423:   PetscCall(PetscViewerDestroy(&viewer));
    424:   PetscCall(VecGetSize(finesolution, &N));
    425:   scale = N / Mx;
    426:   PetscCall(VecGetArrayRead(finesolution, &f));

    428:   /*
    429:      Compute function over the locally owned part of the grid
    430:   */
    431:   for (i = xs; i < xs + xm; i++) {
    432:     x = i * hx;
    433:     r = PetscSqrtReal((x - .5) * (x - .5));
    434:     if (r < .125) u[i] = 1.0;
    435:     else u[i] = -.5;

    437:     /* With the initial condition above the method is first order in space */
    438:     /* this is a smooth initial condition so the method becomes second order in space */
    439:     /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
    440:     u[i] = f[scale * i];
    441:   }
    442:   PetscCall(VecRestoreArrayRead(finesolution, &f));
    443:   PetscCall(VecDestroy(&finesolution));

    445:   /*
    446:      Restore vectors
    447:   */
    448:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    449:   PetscFunctionReturn(PETSC_SUCCESS);
    450: }

    452: /*
    453:     This routine is not parallel
    454: */
    455: PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr)
    456: {
    457:   UserCtx     *ctx = (UserCtx *)ptr;
    458:   PetscDrawLG  lg;
    459:   PetscScalar *u, l, r, c;
    460:   PetscInt     Mx, i, xs, xm, cnt;
    461:   PetscReal    x, y, hx, pause, sx, len, max, xx[4], yy[4], xx_netforce, yy_netforce, yup, ydown, y2, len2;
    462:   PetscDraw    draw;
    463:   Vec          localU;
    464:   DM           da;
    465:   int          colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE, PETSC_DRAW_PLUM, PETSC_DRAW_BLACK};
    466:   /*
    467:   const char *const  legend[3][3] = {{"-kappa (\\grad u,\\grad u)","(1 - u^2)^2"},{"-kappa (\\grad u,\\grad u)","(1 - u^2)"},{"-kappa (\\grad u,\\grad u)","logarithmic"}};
    468:    */
    469:   PetscDrawAxis       axis;
    470:   PetscDrawViewPorts *ports;
    471:   PetscReal           tol = ctx->tol, theta = ctx->theta, theta_c = ctx->theta_c, a, b; /* a and b are used in the cubic truncation of the log function */
    472:   PetscReal           vbounds[] = {-1.1, 1.1};

    474:   PetscFunctionBegin;
    475:   PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds));
    476:   PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 800, 600));
    477:   PetscCall(TSGetDM(ts, &da));
    478:   PetscCall(DMGetLocalVector(da, &localU));
    479:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
    480:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
    481:   hx = 1.0 / (PetscReal)Mx;
    482:   sx = 1.0 / (hx * hx);
    483:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
    484:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
    485:   PetscCall(DMDAVecGetArrayRead(da, localU, &u));

    487:   PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg));
    488:   PetscCall(PetscDrawLGGetDraw(lg, &draw));
    489:   PetscCall(PetscDrawCheckResizedWindow(draw));
    490:   if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports));
    491:   ports = ctx->ports;
    492:   PetscCall(PetscDrawLGGetAxis(lg, &axis));
    493:   PetscCall(PetscDrawLGReset(lg));

    495:   xx[0] = 0.0;
    496:   xx[1] = 1.0;
    497:   cnt   = 2;
    498:   PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL));
    499:   xs = xx[0] / hx;
    500:   xm = (xx[1] - xx[0]) / hx;

    502:   /*
    503:       Plot the  energies
    504:   */
    505:   PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->cahnhillard ? 1 : 0) + (ctx->energy == 3)));
    506:   PetscCall(PetscDrawLGSetColors(lg, colors + 1));
    507:   PetscCall(PetscDrawViewPortsSet(ports, 2));
    508:   x = hx * xs;
    509:   for (i = xs; i < xs + xm; i++) {
    510:     xx[0] = xx[1] = xx[2] = x;
    511:     if (ctx->degenerate) yy[0] = PetscRealPart(.25 * (1. - u[i] * u[i]) * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
    512:     else yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);

    514:     if (ctx->cahnhillard) {
    515:       switch (ctx->energy) {
    516:       case 1: /* double well */
    517:         yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
    518:         break;
    519:       case 2: /* double obstacle */
    520:         yy[1] = .5 * PetscRealPart(1. - u[i] * u[i]);
    521:         break;
    522:       case 3: /* logarithm + double well */
    523:         yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
    524:         if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = .5 * theta * (2.0 * tol * PetscLogReal(tol) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1. - u[i]) / 2.0));
    525:         else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + 2.0 * tol * PetscLogReal(tol));
    526:         else yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1.0 - u[i]) / 2.0));
    527:         break;
    528:       case 4: /* logarithm + double obstacle */
    529:         yy[1] = .5 * theta_c * PetscRealPart(1.0 - u[i] * u[i]);
    530:         if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = .5 * theta * (2.0 * tol * PetscLogReal(tol) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1. - u[i]) / 2.0));
    531:         else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + 2.0 * tol * PetscLogReal(tol));
    532:         else yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1.0 - u[i]) / 2.0));
    533:         break;
    534:       default:
    535:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "It will always be one of the values");
    536:       }
    537:     }
    538:     PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    539:     x += hx;
    540:   }
    541:   PetscCall(PetscDrawGetPause(draw, &pause));
    542:   PetscCall(PetscDrawSetPause(draw, 0.0));
    543:   PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", ""));
    544:   /*  PetscCall(PetscDrawLGSetLegend(lg,legend[ctx->energy-1])); */
    545:   PetscCall(PetscDrawLGDraw(lg));

    547:   /*
    548:       Plot the  forces
    549:   */
    550:   PetscCall(PetscDrawLGSetDimension(lg, 0 + (ctx->cahnhillard ? 2 : 0) + (ctx->energy == 3)));
    551:   PetscCall(PetscDrawLGSetColors(lg, colors + 1));
    552:   PetscCall(PetscDrawViewPortsSet(ports, 1));
    553:   PetscCall(PetscDrawLGReset(lg));
    554:   x   = xs * hx;
    555:   max = 0.;
    556:   for (i = xs; i < xs + xm; i++) {
    557:     xx[0] = xx[1] = xx[2] = xx[3] = x;
    558:     xx_netforce                   = x;
    559:     if (ctx->degenerate) {
    560:       c = (1. - u[i] * u[i]) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    561:       r = (1. - u[i + 1] * u[i + 1]) * (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx;
    562:       l = (1. - u[i - 1] * u[i - 1]) * (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx;
    563:     } else {
    564:       c = (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    565:       r = (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx;
    566:       l = (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx;
    567:     }
    568:     yy[0]       = PetscRealPart(-ctx->kappa * (l + r - 2.0 * c) * sx);
    569:     yy_netforce = yy[0];
    570:     max         = PetscMax(max, PetscAbs(yy[0]));
    571:     if (ctx->cahnhillard) {
    572:       switch (ctx->energy) {
    573:       case 1: /* double well */
    574:         yy[1] = PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    575:         break;
    576:       case 2: /* double obstacle */
    577:         yy[1] = -PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    578:         break;
    579:       case 3: /* logarithmic + double well */
    580:         yy[1] = PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    581:         if (ctx->truncation == 2) { /* quadratic */
    582:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    583:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    584:           else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    585:         } else { /* cubic */
    586:           a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
    587:           b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
    588:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    589:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = PetscRealPart(1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    590:           else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    591:         }
    592:         break;
    593:       case 4: /* logarithmic + double obstacle */
    594:         yy[1] = theta_c * PetscRealPart(-(u[i - 1] + u[i + 1] - 2.0 * u[i])) * sx;
    595:         if (ctx->truncation == 2) {
    596:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    597:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    598:           else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    599:         } else {
    600:           a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
    601:           b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
    602:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    603:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = PetscRealPart(1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    604:           else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    605:         }
    606:         break;
    607:       default:
    608:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "It will always be one of the values");
    609:       }
    610:       if (ctx->energy < 3) {
    611:         max         = PetscMax(max, PetscAbs(yy[1]));
    612:         yy[2]       = yy[0] + yy[1];
    613:         yy_netforce = yy[2];
    614:       } else {
    615:         max         = PetscMax(max, PetscAbs(yy[1] + yy[2]));
    616:         yy[3]       = yy[0] + yy[1] + yy[2];
    617:         yy_netforce = yy[3];
    618:       }
    619:     }
    620:     if (ctx->netforce) {
    621:       PetscCall(PetscDrawLGAddPoint(lg, &xx_netforce, &yy_netforce));
    622:     } else {
    623:       PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    624:     }
    625:     x += hx;
    626:     /*if (max > 7200150000.0) */
    627:     /* printf("max very big when i = %d\n",i); */
    628:   }
    629:   PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", ""));
    630:   PetscCall(PetscDrawLGSetLegend(lg, NULL));
    631:   PetscCall(PetscDrawLGDraw(lg));

    633:   /*
    634:         Plot the solution
    635:   */
    636:   PetscCall(PetscDrawLGSetDimension(lg, 1));
    637:   PetscCall(PetscDrawViewPortsSet(ports, 0));
    638:   PetscCall(PetscDrawLGReset(lg));
    639:   x = hx * xs;
    640:   PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1));
    641:   PetscCall(PetscDrawLGSetColors(lg, colors));
    642:   for (i = xs; i < xs + xm; i++) {
    643:     xx[0] = x;
    644:     yy[0] = PetscRealPart(u[i]);
    645:     PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    646:     x += hx;
    647:   }
    648:   PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", ""));
    649:   PetscCall(PetscDrawLGDraw(lg));

    651:   /*
    652:       Print the  forces as arrows on the solution
    653:   */
    654:   x   = hx * xs;
    655:   cnt = xm / 60;
    656:   cnt = (!cnt) ? 1 : cnt;

    658:   for (i = xs; i < xs + xm; i += cnt) {
    659:     y = yup = ydown = PetscRealPart(u[i]);
    660:     c               = (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
    661:     r               = (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx;
    662:     l               = (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx;
    663:     len             = -.5 * PetscRealPart(ctx->kappa * (l + r - 2.0 * c) * sx) / max;
    664:     PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED));
    665:     if (ctx->cahnhillard) {
    666:       if (len < 0.) ydown += len;
    667:       else yup += len;

    669:       switch (ctx->energy) {
    670:       case 1: /* double well */
    671:         len = .5 * PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    672:         break;
    673:       case 2: /* double obstacle */
    674:         len = -.5 * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
    675:         break;
    676:       case 3: /* logarithmic + double well */
    677:         len = .5 * PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    678:         if (len < 0.) ydown += len;
    679:         else yup += len;

    681:         if (ctx->truncation == 2) { /* quadratic */
    682:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
    683:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
    684:           else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
    685:         } else { /* cubic */
    686:           a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
    687:           b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
    688:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = PetscRealPart(.5 * (-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
    689:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = PetscRealPart(.5 * (a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
    690:           else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
    691:         }
    692:         y2 = len < 0 ? ydown : yup;
    693:         PetscCall(PetscDrawArrow(draw, x, y2, x, y2 + len2, PETSC_DRAW_PLUM));
    694:         break;
    695:       case 4: /* logarithmic + double obstacle */
    696:         len = -.5 * theta_c * PetscRealPart(-(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max);
    697:         if (len < 0.) ydown += len;
    698:         else yup += len;

    700:         if (ctx->truncation == 2) { /* quadratic */
    701:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
    702:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
    703:           else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
    704:         } else { /* cubic */
    705:           a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
    706:           b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
    707:           if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = .5 * PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    708:           else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = .5 * PetscRealPart(a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    709:           else len2 = .5 * PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    710:         }
    711:         y2 = len < 0 ? ydown : yup;
    712:         PetscCall(PetscDrawArrow(draw, x, y2, x, y2 + len2, PETSC_DRAW_PLUM));
    713:         break;
    714:       }
    715:       PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE));
    716:     }
    717:     x += cnt * hx;
    718:   }
    719:   PetscCall(DMDAVecRestoreArrayRead(da, localU, &x));
    720:   PetscCall(DMRestoreLocalVector(da, &localU));
    721:   PetscCall(PetscDrawStringSetSize(draw, .2, .2));
    722:   PetscCall(PetscDrawFlush(draw));
    723:   PetscCall(PetscDrawSetPause(draw, pause));
    724:   PetscCall(PetscDrawPause(draw));
    725:   PetscFunctionReturn(PETSC_SUCCESS);
    726: }

    728: PetscErrorCode MyDestroy(void **ptr)
    729: {
    730:   UserCtx *ctx = *(UserCtx **)ptr;

    732:   PetscFunctionBegin;
    733:   PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
    734:   PetscFunctionReturn(PETSC_SUCCESS);
    735: }

    737: /*TEST

    739:    test:
    740:      TODO: currently requires initial condition file generated by heat

    742: TEST*/