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