1 static char help[] = "Solves biharmonic equation in 1d.\n";
2
3 /*
4 Solves the equation
5
6 u_t = - kappa \Delta \Delta u
7 Periodic boundary conditions
8
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
12
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
16
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
20
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_time_step 5.96046e-06 -cahn-hillard -ts_monitor_draw_solution --mymonitor
24
25 Initial hump neither shrinks nor grows when degenerate (otherwise similar solution)
26
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_time_step 5.96046e-06 -cahn-hillard -degenerate -ts_monitor_draw_solution --mymonitor
28
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_time_step 5.96046e-06 -cahn-hillard -snes_vi_ignore_function_sign -ts_monitor_draw_solution --mymonitor
30
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_time_step 5.96046e-06 -cahn-hillard -energy 2 -snes_linesearch_monitor -ts_monitor_draw_solution --mymonitor
34
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_time_step 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -mymonitor
38
39 ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_time_step 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -degenerate -mymonitor
40
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_time_step 5.96046e-06 -cahn-hillard -energy 4 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --mymonitor
44
45 */
46 #include <petscdm.h>
47 #include <petscdmda.h>
48 #include <petscts.h>
49 #include <petscdraw.h>
50
51 extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), FormJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
52 extern PetscCtxDestroyFn MyDestroy;
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;
64
main(int argc,char ** argv)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;
75
76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77 Initialize program
78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79 PetscFunctionBeginUser;
80 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
100
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);
110
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));
117
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));
130
131 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 Create matrix data structure; set Jacobian evaluation routine
133
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 matrix used to construct the preconditioner as set by the user,
139 but use matrix-free approx for Jacobian-vector
140 products within Newton-Krylov method
141
142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 Customize nonlinear solver
145 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146 PetscCall(TSSetType(ts, TSCN));
147
148 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149 Set initial conditions
150 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151 PetscCall(FormInitialSolution(da, x));
152 PetscCall(TSSetTimeStep(ts, dt));
153 PetscCall(TSSetSolution(ts, x));
154
155 if (mymonitor) {
156 ctx.ports = NULL;
157 PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
158 }
159
160 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161 Set runtime options
162 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163 PetscCall(TSSetFromOptions(ts));
164
165 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166 Solve nonlinear system
167 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168 PetscCall(TSSolve(ts, x));
169 PetscCall(TSGetStepNumber(ts, &steps));
170 PetscCall(VecView(x, PETSC_VIEWER_BINARY_WORLD));
171
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));
181
182 PetscCall(PetscFinalize());
183 return 0;
184 }
185 /* ------------------------------------------------------------------- */
186 /*
187 FormFunction - Evaluates nonlinear function, F(x).
188
189 Input Parameters:
190 . ts - the TS context
191 . X - input vector
192 . ptr - optional user-defined context, as set by SNESSetFunction()
193
194 Output Parameter:
195 . F - function vector
196 */
FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void * ptr)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 */
206
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));
211
212 hx = 1.0 / (PetscReal)Mx;
213 sx = 1.0 / (hx * hx);
214
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));
223
224 /*
225 Get pointers to vector data
226 */
227 PetscCall(DMDAVecGetArrayRead(da, localX, &x));
228 PetscCall(DMDAVecGetArray(da, F, &f));
229
230 /*
231 Get local grid boundaries
232 */
233 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
234
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 }
288
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 }
297
298 /* ------------------------------------------------------------------- */
299 /*
300 FormJacobian - Evaluates nonlinear function's Jacobian
301
302 */
FormJacobian(TS ts,PetscReal ftime,Vec X,Mat A,Mat B,void * ptr)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;
312
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));
317
318 hx = 1.0 / (PetscReal)Mx;
319 sx = 1.0 / (hx * hx);
320
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));
329
330 /*
331 Get pointers to vector data
332 */
333 PetscCall(DMDAVecGetArrayRead(da, localX, &x));
334
335 /*
336 Get local grid boundaries
337 */
338 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
339
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));
363
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 }
379
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 /* ------------------------------------------------------------------- */
FormInitialSolution(DM da,Vec U)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;
402
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));
405
406 hx = 1.0 / (PetscReal)Mx;
407
408 /*
409 Get pointers to vector data
410 */
411 PetscCall(DMDAVecGetArray(da, U, &u));
412
413 /*
414 Get local grid boundaries
415 */
416 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
417
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));
428
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;
437
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));
445
446 /*
447 Restore vectors
448 */
449 PetscCall(DMDAVecRestoreArray(da, U, &u));
450 PetscFunctionReturn(PETSC_SUCCESS);
451 }
452
453 /*
454 This routine is not parallel
455 */
MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,void * ptr)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};
474
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));
487
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));
495
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;
502
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);
514
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));
547
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));
633
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));
651
652 /*
653 Print the forces as arrows on the solution
654 */
655 x = hx * xs;
656 cnt = xm / 60;
657 cnt = (!cnt) ? 1 : cnt;
658
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;
669
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;
681
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;
700
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 }
728
MyDestroy(PetscCtxRt Ctx)729 PetscErrorCode MyDestroy(PetscCtxRt Ctx)
730 {
731 UserCtx *ctx = *(UserCtx **)Ctx;
732
733 PetscFunctionBegin;
734 PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
735 PetscFunctionReturn(PETSC_SUCCESS);
736 }
737
738 /*TEST
739
740 test:
741 TODO: currently requires initial condition file generated by heat
742
743 TEST*/
744