Actual source code: biharmonic.c


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

  4: /*
  5:   Solves the equation

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

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

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

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

 22: Evolve the Cahn-Hillard equations: double well Initial hump shrinks then grows
 23: ---------------
 24: ./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

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

 28: ./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

 30: ./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

 32: Evolve the Cahn-Hillard equations: double obstacle
 33: ---------------
 34: ./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

 36: Evolve the Cahn-Hillard equations: logarithmic + double well (never shrinks and then grows)
 37: ---------------
 38: ./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

 40: ./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

 42: Evolve the Cahn-Hillard equations: logarithmic +  double obstacle (never shrinks, never grows)
 43: ---------------
 44: ./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

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

 52: 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 *);
 53: typedef struct {
 54:   PetscBool           cahnhillard;
 55:   PetscBool           degenerate;
 56:   PetscReal           kappa;
 57:   PetscInt            energy;
 58:   PetscReal           tol;
 59:   PetscReal           theta, theta_c;
 60:   PetscInt            truncation;
 61:   PetscBool           netforce;
 62:   PetscDrawViewPorts *ports;
 63: } UserCtx;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

194:    Output Parameter:
195: .  F - function vector
196:  */
197: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
198: {
199:   DM           da;
200:   PetscInt     i, Mx, xs, xm;
201:   PetscReal    hx, sx;
202:   PetscScalar *x, *f, c, r, l;
203:   Vec          localX;
204:   UserCtx     *ctx = (UserCtx *)ptr;
205:   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 */

207:   PetscFunctionBegin;
208:   PetscCall(TSGetDM(ts, &da));
209:   PetscCall(DMGetLocalVector(da, &localX));
210:   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));

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

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

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

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

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

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

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

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

313:   PetscFunctionBegin;
314:   PetscCall(TSGetDM(ts, &da));
315:   PetscCall(DMGetLocalVector(da, &localX));
316:   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));

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

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

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

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

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

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

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

403:   PetscFunctionBegin;
404:   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));

406:   hx = 1.0 / (PetscReal)Mx;

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

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

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

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

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

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

453: /*
454:     This routine is not parallel
455: */
456: PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr)
457: {
458:   UserCtx     *ctx = (UserCtx *)ptr;
459:   PetscDrawLG  lg;
460:   PetscScalar *u, l, r, c;
461:   PetscInt     Mx, i, xs, xm, cnt;
462:   PetscReal    x, y, hx, pause, sx, len, max, xx[4], yy[4], xx_netforce, yy_netforce, yup, ydown, y2, len2;
463:   PetscDraw    draw;
464:   Vec          localU;
465:   DM           da;
466:   int          colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE, PETSC_DRAW_PLUM, PETSC_DRAW_BLACK};
467:   /*
468:   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"}};
469:    */
470:   PetscDrawAxis       axis;
471:   PetscDrawViewPorts *ports;
472:   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 */
473:   PetscReal           vbounds[] = {-1.1, 1.1};

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

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

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

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

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

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

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

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

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

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

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

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

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

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

738: /*TEST

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

743: TEST*/