Actual source code: biharmonic.c
1: static char help[] = "Solves biharmonic equation in 1d.\n";
3: /*
4: Solves the equation
6: u_t = - kappa \Delta \Delta u
7: Periodic boundary conditions
9: Evolve the biharmonic heat equation:
10: ---------------
11: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor
13: Evolve with the restriction that -1 <= u <= 1; i.e. as a variational inequality
14: ---------------
15: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor
17: u_t = kappa \Delta \Delta u + 6.*u*(u_x)^2 + (3*u^2 - 12) \Delta u
18: -1 <= u <= 1
19: Periodic boundary conditions
21: Evolve the Cahn-Hillard equations: double well Initial hump shrinks then grows
22: ---------------
23: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -ts_monitor_draw_solution --mymonitor
25: Initial hump neither shrinks nor grows when degenerate (otherwise similar solution)
27: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -degenerate -ts_monitor_draw_solution --mymonitor
29: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -snes_vi_ignore_function_sign -ts_monitor_draw_solution --mymonitor
31: Evolve the Cahn-Hillard equations: double obstacle
32: ---------------
33: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 2 -snes_linesearch_monitor -ts_monitor_draw_solution --mymonitor
35: Evolve the Cahn-Hillard equations: logarithmic + double well (never shrinks and then grows)
36: ---------------
37: ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -mymonitor
39: ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -degenerate -mymonitor
41: Evolve the Cahn-Hillard equations: logarithmic + double obstacle (never shrinks, never grows)
42: ---------------
43: ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 4 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --mymonitor
45: */
46: #include <petscdm.h>
47: #include <petscdmda.h>
48: #include <petscts.h>
49: #include <petscdraw.h>
51: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), MyDestroy(void **), FormJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
52: typedef struct {
53: PetscBool cahnhillard;
54: PetscBool degenerate;
55: PetscReal kappa;
56: PetscInt energy;
57: PetscReal tol;
58: PetscReal theta, theta_c;
59: PetscInt truncation;
60: PetscBool netforce;
61: PetscDrawViewPorts *ports;
62: } UserCtx;
64: int main(int argc, char **argv)
65: {
66: TS ts; /* nonlinear solver */
67: Vec x, r; /* solution, residual vectors */
68: Mat J; /* Jacobian matrix */
69: PetscInt steps, Mx;
70: DM da;
71: PetscReal dt;
72: PetscBool mymonitor;
73: UserCtx ctx;
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Initialize program
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: PetscFunctionBeginUser;
79: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
80: ctx.kappa = 1.0;
81: PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL));
82: ctx.degenerate = PETSC_FALSE;
83: PetscCall(PetscOptionsGetBool(NULL, NULL, "-degenerate", &ctx.degenerate, NULL));
84: ctx.cahnhillard = PETSC_FALSE;
85: PetscCall(PetscOptionsGetBool(NULL, NULL, "-cahn-hillard", &ctx.cahnhillard, NULL));
86: ctx.netforce = PETSC_FALSE;
87: PetscCall(PetscOptionsGetBool(NULL, NULL, "-netforce", &ctx.netforce, NULL));
88: ctx.energy = 1;
89: PetscCall(PetscOptionsGetInt(NULL, NULL, "-energy", &ctx.energy, NULL));
90: ctx.tol = 1.0e-8;
91: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &ctx.tol, NULL));
92: ctx.theta = .001;
93: ctx.theta_c = 1.0;
94: PetscCall(PetscOptionsGetReal(NULL, NULL, "-theta", &ctx.theta, NULL));
95: PetscCall(PetscOptionsGetReal(NULL, NULL, "-theta_c", &ctx.theta_c, NULL));
96: ctx.truncation = 1;
97: PetscCall(PetscOptionsGetInt(NULL, NULL, "-truncation", &ctx.truncation, NULL));
98: PetscCall(PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor));
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create distributed array (DMDA) to manage parallel grid and vectors
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da));
104: PetscCall(DMSetFromOptions(da));
105: PetscCall(DMSetUp(da));
106: PetscCall(DMDASetFieldName(da, 0, "Biharmonic heat equation: u"));
107: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
108: dt = 1.0 / (10. * ctx.kappa * Mx * Mx * Mx * Mx);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Extract global vectors from DMDA; then duplicate for remaining
112: vectors that are the same types
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: PetscCall(DMCreateGlobalVector(da, &x));
115: PetscCall(VecDuplicate(x, &r));
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Create timestepping solver context
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
121: PetscCall(TSSetDM(ts, da));
122: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
123: PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &ctx));
124: PetscCall(DMSetMatType(da, MATAIJ));
125: PetscCall(DMCreateMatrix(da, &J));
126: PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &ctx));
127: PetscCall(TSSetMaxTime(ts, .02));
128: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE));
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Create matrix data structure; set Jacobian evaluation routine
133: Set Jacobian matrix data structure and default Jacobian evaluation
134: routine. User can override with:
135: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
136: (unless user explicitly sets preconditioner)
137: -snes_mf_operator : form preconditioning matrix as set by the user,
138: but use matrix-free approx for Jacobian-vector
139: products within Newton-Krylov method
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Customize nonlinear solver
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: PetscCall(TSSetType(ts, TSCN));
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Set initial conditions
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: PetscCall(FormInitialSolution(da, x));
151: PetscCall(TSSetTimeStep(ts, dt));
152: PetscCall(TSSetSolution(ts, x));
154: if (mymonitor) {
155: ctx.ports = NULL;
156: PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
157: }
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Set runtime options
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: PetscCall(TSSetFromOptions(ts));
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Solve nonlinear system
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: PetscCall(TSSolve(ts, x));
168: PetscCall(TSGetStepNumber(ts, &steps));
169: PetscCall(VecView(x, PETSC_VIEWER_BINARY_WORLD));
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Free work space. All PETSc objects should be destroyed when they
173: are no longer needed.
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: PetscCall(MatDestroy(&J));
176: PetscCall(VecDestroy(&x));
177: PetscCall(VecDestroy(&r));
178: PetscCall(TSDestroy(&ts));
179: PetscCall(DMDestroy(&da));
181: PetscCall(PetscFinalize());
182: return 0;
183: }
184: /* ------------------------------------------------------------------- */
185: /*
186: FormFunction - Evaluates nonlinear function, F(x).
188: Input Parameters:
189: . ts - the TS context
190: . X - input vector
191: . ptr - optional user-defined context, as set by SNESSetFunction()
193: Output Parameter:
194: . F - function vector
195: */
196: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
197: {
198: DM da;
199: PetscInt i, Mx, xs, xm;
200: PetscReal hx, sx;
201: PetscScalar *x, *f, c, r, l;
202: Vec localX;
203: UserCtx *ctx = (UserCtx *)ptr;
204: PetscReal tol = ctx->tol, theta = ctx->theta, theta_c = ctx->theta_c, a, b; /* a and b are used in the cubic truncation of the log function */
206: PetscFunctionBegin;
207: PetscCall(TSGetDM(ts, &da));
208: PetscCall(DMGetLocalVector(da, &localX));
209: PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
211: hx = 1.0 / (PetscReal)Mx;
212: sx = 1.0 / (hx * hx);
214: /*
215: Scatter ghost points to local vector,using the 2-step process
216: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
217: By placing code between these two statements, computations can be
218: done while messages are in transition.
219: */
220: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
221: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
223: /*
224: Get pointers to vector data
225: */
226: PetscCall(DMDAVecGetArrayRead(da, localX, &x));
227: PetscCall(DMDAVecGetArray(da, F, &f));
229: /*
230: Get local grid boundaries
231: */
232: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
234: /*
235: Compute function over the locally owned part of the grid
236: */
237: for (i = xs; i < xs + xm; i++) {
238: if (ctx->degenerate) {
239: c = (1. - x[i] * x[i]) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
240: r = (1. - x[i + 1] * x[i + 1]) * (x[i] + x[i + 2] - 2.0 * x[i + 1]) * sx;
241: l = (1. - x[i - 1] * x[i - 1]) * (x[i - 2] + x[i] - 2.0 * x[i - 1]) * sx;
242: } else {
243: c = (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
244: r = (x[i] + x[i + 2] - 2.0 * x[i + 1]) * sx;
245: l = (x[i - 2] + x[i] - 2.0 * x[i - 1]) * sx;
246: }
247: f[i] = -ctx->kappa * (l + r - 2.0 * c) * sx;
248: if (ctx->cahnhillard) {
249: switch (ctx->energy) {
250: case 1: /* double well */
251: f[i] += 6. * .25 * x[i] * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (3. * x[i] * x[i] - 1.) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
252: break;
253: case 2: /* double obstacle */
254: f[i] += -(x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
255: break;
256: case 3: /* logarithmic + double well */
257: f[i] += 6. * .25 * x[i] * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (3. * x[i] * x[i] - 1.) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
258: if (ctx->truncation == 2) { /* log function with approximated with a quadratic polynomial outside -1.0+2*tol, 1.0-2*tol */
259: if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
260: else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
261: else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
262: } else { /* log function is approximated with a cubic polynomial outside -1.0+2*tol, 1.0-2*tol */
263: a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
264: b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
265: if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += -1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (-1.0 * a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
266: else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += 1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
267: else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
268: }
269: break;
270: case 4: /* logarithmic + double obstacle */
271: f[i] += -theta_c * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
272: if (ctx->truncation == 2) { /* quadratic */
273: if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
274: else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += (.25 * theta / (tol - tol * tol)) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
275: else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
276: } else { /* cubic */
277: a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
278: b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
279: if (PetscRealPart(x[i]) < -1.0 + 2.0 * tol) f[i] += -1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (-1.0 * a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
280: else if (PetscRealPart(x[i]) > 1.0 - 2.0 * tol) f[i] += 1.0 * a * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (a * x[i] + b) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
281: else f[i] += 2.0 * theta * x[i] / ((1.0 - x[i] * x[i]) * (1.0 - x[i] * x[i])) * .25 * (x[i + 1] - x[i - 1]) * (x[i + 1] - x[i - 1]) * sx + (theta / (1.0 - x[i] * x[i])) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
282: }
283: break;
284: }
285: }
286: }
288: /*
289: Restore vectors
290: */
291: PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
292: PetscCall(DMDAVecRestoreArray(da, F, &f));
293: PetscCall(DMRestoreLocalVector(da, &localX));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /* ------------------------------------------------------------------- */
298: /*
299: FormJacobian - Evaluates nonlinear function's Jacobian
301: */
302: PetscErrorCode FormJacobian(TS ts, PetscReal ftime, Vec X, Mat A, Mat B, void *ptr)
303: {
304: DM da;
305: PetscInt i, Mx, xs, xm;
306: MatStencil row, cols[5];
307: PetscReal hx, sx;
308: PetscScalar *x, vals[5];
309: Vec localX;
310: UserCtx *ctx = (UserCtx *)ptr;
312: PetscFunctionBegin;
313: PetscCall(TSGetDM(ts, &da));
314: PetscCall(DMGetLocalVector(da, &localX));
315: PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
317: hx = 1.0 / (PetscReal)Mx;
318: sx = 1.0 / (hx * hx);
320: /*
321: Scatter ghost points to local vector,using the 2-step process
322: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
323: By placing code between these two statements, computations can be
324: done while messages are in transition.
325: */
326: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
327: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
329: /*
330: Get pointers to vector data
331: */
332: PetscCall(DMDAVecGetArrayRead(da, localX, &x));
334: /*
335: Get local grid boundaries
336: */
337: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
339: /*
340: Compute function over the locally owned part of the grid
341: */
342: for (i = xs; i < xs + xm; i++) {
343: row.i = i;
344: if (ctx->degenerate) {
345: /*PetscScalar c,r,l;
346: c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
347: r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx;
348: l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx; */
349: } else {
350: cols[0].i = i - 2;
351: vals[0] = -ctx->kappa * sx * sx;
352: cols[1].i = i - 1;
353: vals[1] = 4.0 * ctx->kappa * sx * sx;
354: cols[2].i = i;
355: vals[2] = -6.0 * ctx->kappa * sx * sx;
356: cols[3].i = i + 1;
357: vals[3] = 4.0 * ctx->kappa * sx * sx;
358: cols[4].i = i + 2;
359: vals[4] = -ctx->kappa * sx * sx;
360: }
361: PetscCall(MatSetValuesStencil(B, 1, &row, 5, cols, vals, INSERT_VALUES));
363: if (ctx->cahnhillard) {
364: switch (ctx->energy) {
365: case 1: /* double well */
366: /* f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; */
367: break;
368: case 2: /* double obstacle */
369: /* f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx; */
370: break;
371: case 3: /* logarithmic + double well */
372: break;
373: case 4: /* logarithmic + double obstacle */
374: break;
375: }
376: }
377: }
379: /*
380: Restore vectors
381: */
382: PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
383: PetscCall(DMRestoreLocalVector(da, &localX));
384: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
385: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
386: if (A != B) {
387: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
388: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
389: }
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
392: /* ------------------------------------------------------------------- */
393: PetscErrorCode FormInitialSolution(DM da, Vec U)
394: {
395: PetscInt i, xs, xm, Mx, N, scale;
396: PetscScalar *u;
397: PetscReal r, hx, x;
398: const PetscScalar *f;
399: Vec finesolution;
400: PetscViewer viewer;
402: PetscFunctionBegin;
403: PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
405: hx = 1.0 / (PetscReal)Mx;
407: /*
408: Get pointers to vector data
409: */
410: PetscCall(DMDAVecGetArray(da, U, &u));
412: /*
413: Get local grid boundaries
414: */
415: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
417: /*
418: Seee heat.c for how to generate InitialSolution.heat
419: */
420: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer));
421: PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution));
422: PetscCall(VecLoad(finesolution, viewer));
423: PetscCall(PetscViewerDestroy(&viewer));
424: PetscCall(VecGetSize(finesolution, &N));
425: scale = N / Mx;
426: PetscCall(VecGetArrayRead(finesolution, &f));
428: /*
429: Compute function over the locally owned part of the grid
430: */
431: for (i = xs; i < xs + xm; i++) {
432: x = i * hx;
433: r = PetscSqrtReal((x - .5) * (x - .5));
434: if (r < .125) u[i] = 1.0;
435: else u[i] = -.5;
437: /* With the initial condition above the method is first order in space */
438: /* this is a smooth initial condition so the method becomes second order in space */
439: /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
440: u[i] = f[scale * i];
441: }
442: PetscCall(VecRestoreArrayRead(finesolution, &f));
443: PetscCall(VecDestroy(&finesolution));
445: /*
446: Restore vectors
447: */
448: PetscCall(DMDAVecRestoreArray(da, U, &u));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*
453: This routine is not parallel
454: */
455: PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr)
456: {
457: UserCtx *ctx = (UserCtx *)ptr;
458: PetscDrawLG lg;
459: PetscScalar *u, l, r, c;
460: PetscInt Mx, i, xs, xm, cnt;
461: PetscReal x, y, hx, pause, sx, len, max, xx[4], yy[4], xx_netforce, yy_netforce, yup, ydown, y2, len2;
462: PetscDraw draw;
463: Vec localU;
464: DM da;
465: int colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE, PETSC_DRAW_PLUM, PETSC_DRAW_BLACK};
466: /*
467: const char *const legend[3][3] = {{"-kappa (\\grad u,\\grad u)","(1 - u^2)^2"},{"-kappa (\\grad u,\\grad u)","(1 - u^2)"},{"-kappa (\\grad u,\\grad u)","logarithmic"}};
468: */
469: PetscDrawAxis axis;
470: PetscDrawViewPorts *ports;
471: PetscReal tol = ctx->tol, theta = ctx->theta, theta_c = ctx->theta_c, a, b; /* a and b are used in the cubic truncation of the log function */
472: PetscReal vbounds[] = {-1.1, 1.1};
474: PetscFunctionBegin;
475: PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds));
476: PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 800, 600));
477: PetscCall(TSGetDM(ts, &da));
478: PetscCall(DMGetLocalVector(da, &localU));
479: PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
480: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
481: hx = 1.0 / (PetscReal)Mx;
482: sx = 1.0 / (hx * hx);
483: PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
484: PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
485: PetscCall(DMDAVecGetArrayRead(da, localU, &u));
487: PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg));
488: PetscCall(PetscDrawLGGetDraw(lg, &draw));
489: PetscCall(PetscDrawCheckResizedWindow(draw));
490: if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports));
491: ports = ctx->ports;
492: PetscCall(PetscDrawLGGetAxis(lg, &axis));
493: PetscCall(PetscDrawLGReset(lg));
495: xx[0] = 0.0;
496: xx[1] = 1.0;
497: cnt = 2;
498: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL));
499: xs = xx[0] / hx;
500: xm = (xx[1] - xx[0]) / hx;
502: /*
503: Plot the energies
504: */
505: PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->cahnhillard ? 1 : 0) + (ctx->energy == 3)));
506: PetscCall(PetscDrawLGSetColors(lg, colors + 1));
507: PetscCall(PetscDrawViewPortsSet(ports, 2));
508: x = hx * xs;
509: for (i = xs; i < xs + xm; i++) {
510: xx[0] = xx[1] = xx[2] = x;
511: if (ctx->degenerate) yy[0] = PetscRealPart(.25 * (1. - u[i] * u[i]) * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
512: else yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
514: if (ctx->cahnhillard) {
515: switch (ctx->energy) {
516: case 1: /* double well */
517: yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
518: break;
519: case 2: /* double obstacle */
520: yy[1] = .5 * PetscRealPart(1. - u[i] * u[i]);
521: break;
522: case 3: /* logarithm + double well */
523: yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
524: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = .5 * theta * (2.0 * tol * PetscLogReal(tol) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1. - u[i]) / 2.0));
525: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + 2.0 * tol * PetscLogReal(tol));
526: else yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1.0 - u[i]) / 2.0));
527: break;
528: case 4: /* logarithm + double obstacle */
529: yy[1] = .5 * theta_c * PetscRealPart(1.0 - u[i] * u[i]);
530: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = .5 * theta * (2.0 * tol * PetscLogReal(tol) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1. - u[i]) / 2.0));
531: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + 2.0 * tol * PetscLogReal(tol));
532: else yy[2] = .5 * theta * (PetscRealPart(1.0 + u[i]) * PetscLogReal(PetscRealPart(1.0 + u[i]) / 2.0) + PetscRealPart(1.0 - u[i]) * PetscLogReal(PetscRealPart(1.0 - u[i]) / 2.0));
533: break;
534: default:
535: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "It will always be one of the values");
536: }
537: }
538: PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
539: x += hx;
540: }
541: PetscCall(PetscDrawGetPause(draw, &pause));
542: PetscCall(PetscDrawSetPause(draw, 0.0));
543: PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", ""));
544: /* PetscCall(PetscDrawLGSetLegend(lg,legend[ctx->energy-1])); */
545: PetscCall(PetscDrawLGDraw(lg));
547: /*
548: Plot the forces
549: */
550: PetscCall(PetscDrawLGSetDimension(lg, 0 + (ctx->cahnhillard ? 2 : 0) + (ctx->energy == 3)));
551: PetscCall(PetscDrawLGSetColors(lg, colors + 1));
552: PetscCall(PetscDrawViewPortsSet(ports, 1));
553: PetscCall(PetscDrawLGReset(lg));
554: x = xs * hx;
555: max = 0.;
556: for (i = xs; i < xs + xm; i++) {
557: xx[0] = xx[1] = xx[2] = xx[3] = x;
558: xx_netforce = x;
559: if (ctx->degenerate) {
560: c = (1. - u[i] * u[i]) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
561: r = (1. - u[i + 1] * u[i + 1]) * (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx;
562: l = (1. - u[i - 1] * u[i - 1]) * (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx;
563: } else {
564: c = (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
565: r = (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx;
566: l = (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx;
567: }
568: yy[0] = PetscRealPart(-ctx->kappa * (l + r - 2.0 * c) * sx);
569: yy_netforce = yy[0];
570: max = PetscMax(max, PetscAbs(yy[0]));
571: if (ctx->cahnhillard) {
572: switch (ctx->energy) {
573: case 1: /* double well */
574: yy[1] = PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
575: break;
576: case 2: /* double obstacle */
577: yy[1] = -PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
578: break;
579: case 3: /* logarithmic + double well */
580: yy[1] = PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
581: if (ctx->truncation == 2) { /* quadratic */
582: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
583: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
584: else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
585: } else { /* cubic */
586: a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
587: b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
588: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
589: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = PetscRealPart(1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
590: else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
591: }
592: break;
593: case 4: /* logarithmic + double obstacle */
594: yy[1] = theta_c * PetscRealPart(-(u[i - 1] + u[i + 1] - 2.0 * u[i])) * sx;
595: if (ctx->truncation == 2) {
596: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
597: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
598: else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
599: } else {
600: a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
601: b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
602: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
603: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = PetscRealPart(1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
604: else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
605: }
606: break;
607: default:
608: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "It will always be one of the values");
609: }
610: if (ctx->energy < 3) {
611: max = PetscMax(max, PetscAbs(yy[1]));
612: yy[2] = yy[0] + yy[1];
613: yy_netforce = yy[2];
614: } else {
615: max = PetscMax(max, PetscAbs(yy[1] + yy[2]));
616: yy[3] = yy[0] + yy[1] + yy[2];
617: yy_netforce = yy[3];
618: }
619: }
620: if (ctx->netforce) {
621: PetscCall(PetscDrawLGAddPoint(lg, &xx_netforce, &yy_netforce));
622: } else {
623: PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
624: }
625: x += hx;
626: /*if (max > 7200150000.0) */
627: /* printf("max very big when i = %d\n",i); */
628: }
629: PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", ""));
630: PetscCall(PetscDrawLGSetLegend(lg, NULL));
631: PetscCall(PetscDrawLGDraw(lg));
633: /*
634: Plot the solution
635: */
636: PetscCall(PetscDrawLGSetDimension(lg, 1));
637: PetscCall(PetscDrawViewPortsSet(ports, 0));
638: PetscCall(PetscDrawLGReset(lg));
639: x = hx * xs;
640: PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1));
641: PetscCall(PetscDrawLGSetColors(lg, colors));
642: for (i = xs; i < xs + xm; i++) {
643: xx[0] = x;
644: yy[0] = PetscRealPart(u[i]);
645: PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
646: x += hx;
647: }
648: PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", ""));
649: PetscCall(PetscDrawLGDraw(lg));
651: /*
652: Print the forces as arrows on the solution
653: */
654: x = hx * xs;
655: cnt = xm / 60;
656: cnt = (!cnt) ? 1 : cnt;
658: for (i = xs; i < xs + xm; i += cnt) {
659: y = yup = ydown = PetscRealPart(u[i]);
660: c = (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx;
661: r = (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx;
662: l = (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx;
663: len = -.5 * PetscRealPart(ctx->kappa * (l + r - 2.0 * c) * sx) / max;
664: PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED));
665: if (ctx->cahnhillard) {
666: if (len < 0.) ydown += len;
667: else yup += len;
669: switch (ctx->energy) {
670: case 1: /* double well */
671: len = .5 * PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
672: break;
673: case 2: /* double obstacle */
674: len = -.5 * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
675: break;
676: case 3: /* logarithmic + double well */
677: len = .5 * PetscRealPart(6. * .25 * u[i] * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (3. * u[i] * u[i] - 1.) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
678: if (len < 0.) ydown += len;
679: else yup += len;
681: if (ctx->truncation == 2) { /* quadratic */
682: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
683: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
684: else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
685: } else { /* cubic */
686: a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
687: b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
688: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = PetscRealPart(.5 * (-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
689: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = PetscRealPart(.5 * (a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
690: else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
691: }
692: y2 = len < 0 ? ydown : yup;
693: PetscCall(PetscDrawArrow(draw, x, y2, x, y2 + len2, PETSC_DRAW_PLUM));
694: break;
695: case 4: /* logarithmic + double obstacle */
696: len = -.5 * theta_c * PetscRealPart(-(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max);
697: if (len < 0.) ydown += len;
698: else yup += len;
700: if (ctx->truncation == 2) { /* quadratic */
701: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
702: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = .5 * (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max;
703: else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max);
704: } else { /* cubic */
705: a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol));
706: b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol);
707: if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) len2 = .5 * PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
708: else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) len2 = .5 * PetscRealPart(a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
709: else len2 = .5 * PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
710: }
711: y2 = len < 0 ? ydown : yup;
712: PetscCall(PetscDrawArrow(draw, x, y2, x, y2 + len2, PETSC_DRAW_PLUM));
713: break;
714: }
715: PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE));
716: }
717: x += cnt * hx;
718: }
719: PetscCall(DMDAVecRestoreArrayRead(da, localU, &x));
720: PetscCall(DMRestoreLocalVector(da, &localU));
721: PetscCall(PetscDrawStringSetSize(draw, .2, .2));
722: PetscCall(PetscDrawFlush(draw));
723: PetscCall(PetscDrawSetPause(draw, pause));
724: PetscCall(PetscDrawPause(draw));
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: PetscErrorCode MyDestroy(void **ptr)
729: {
730: UserCtx *ctx = *(UserCtx **)ptr;
732: PetscFunctionBegin;
733: PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: /*TEST
739: test:
740: TODO: currently requires initial condition file generated by heat
742: TEST*/