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*/