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

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Create distributed array (DMDA) to manage parallel grid and vectors
103:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da);
105:   DMSetFromOptions(da);
106:   DMSetUp(da);
107:   DMDASetFieldName(da, 0, "Biharmonic heat equation: u");
108:   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:   DMCreateGlobalVector(da, &x);
116:   VecDuplicate(x, &r);

118:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119:      Create timestepping solver context
120:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121:   TSCreate(PETSC_COMM_WORLD, &ts);
122:   TSSetDM(ts, da);
123:   TSSetProblemType(ts, TS_NONLINEAR);
124:   TSSetRHSFunction(ts, NULL, FormFunction, &ctx);
125:   DMSetMatType(da, MATAIJ);
126:   DMCreateMatrix(da, &J);
127:   TSSetRHSJacobian(ts, J, J, FormJacobian, &ctx);
128:   TSSetMaxTime(ts, .02);
129:   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:   TSSetType(ts, TSCN);

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

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

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

165:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166:      Solve nonlinear system
167:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168:   TSSolve(ts, x);
169:   TSGetStepNumber(ts, &steps);
170:   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:   MatDestroy(&J);
177:   VecDestroy(&x);
178:   VecDestroy(&r);
179:   TSDestroy(&ts);
180:   DMDestroy(&da);

182:   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:   TSGetDM(ts, &da);
208:   DMGetLocalVector(da, &localX);
209:   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:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
221:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);

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

229:   /*
230:      Get local grid boundaries
231:   */
232:   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:   DMDAVecRestoreArrayRead(da, localX, &x);
292:   DMDAVecRestoreArray(da, F, &f);
293:   DMRestoreLocalVector(da, &localX);
294:   return 0;
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:   TSGetDM(ts, &da);
313:   DMGetLocalVector(da, &localX);
314:   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);

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

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

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

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

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

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

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

401:   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);

403:   hx = 1.0 / (PetscReal)Mx;

405:   /*
406:      Get pointers to vector data
407:   */
408:   DMDAVecGetArray(da, U, &u);

410:   /*
411:      Get local grid boundaries
412:   */
413:   DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);

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

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

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

443:   /*
444:      Restore vectors
445:   */
446:   DMDAVecRestoreArray(da, U, &u);
447:   return 0;
448: }

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

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

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

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

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

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

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

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

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

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

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

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

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

725: PetscErrorCode MyDestroy(void **ptr)
726: {
727:   UserCtx *ctx = *(UserCtx **)ptr;

729:   PetscDrawViewPortsDestroy(ctx->ports);
730:   return 0;
731: }

733: /*TEST

735:    test:
736:      TODO: currently requires initial condition file generated by heat

738: TEST*/