1 static char help[] = "Solves heat equation in 1d.\n";
2
3 /*
4 Solves the equation
5
6 u_t = kappa \Delta u
7 Periodic boundary conditions
8
9 Evolve the heat equation:
10 ---------------
11 ./heat -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor
12
13 Evolve the Allen-Cahn equation:
14 ---------------
15 ./heat -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -ts_type cn -da_refine 5 -allen-cahn -kappa .001 -ts_max_time 5 -mymonitor
16
17 Evolve the Allen-Cahn equation: zoom in on part of the domain
18 ---------------
19 ./heat -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -allen-cahn -kappa .001 -ts_max_time 5 -zoom .25,.45 -mymonitor
20
21 The option -square_initial indicates it should use a square wave initial condition otherwise it loads the file InitialSolution.heat as the initial solution. You should run with
22 ./heat -square_initial -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 9 -ts_max_time 1.e-4 -ts_time_step .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25 -ts_max_steps 15
23 to generate InitialSolution.heat
24
25 */
26 #include <petscdm.h>
27 #include <petscdmda.h>
28 #include <petscts.h>
29 #include <petscdraw.h>
30
31 /*
32 User-defined routines
33 */
34 extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, PetscCtx), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, PetscCtx);
35 extern PetscCtxDestroyFn MyDestroy;
36 typedef struct {
37 PetscReal kappa;
38 PetscBool allencahn;
39 PetscDrawViewPorts *ports;
40 } AppCtx;
41
main(int argc,char ** argv)42 int main(int argc, char **argv)
43 {
44 TS ts; /* time integrator */
45 Vec x, r; /* solution, residual vectors */
46 PetscInt steps, Mx;
47 DM da;
48 PetscReal dt;
49 AppCtx ctx;
50 PetscBool mymonitor;
51 PetscViewer viewer;
52 PetscBool flg;
53
54 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55 Initialize program
56 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57 PetscFunctionBeginUser;
58 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
59 ctx.kappa = 1.0;
60 PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL));
61 ctx.allencahn = PETSC_FALSE;
62 PetscCall(PetscOptionsHasName(NULL, NULL, "-allen-cahn", &ctx.allencahn));
63 PetscCall(PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor));
64
65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66 Create distributed array (DMDA) to manage parallel grid and vectors
67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da));
69 PetscCall(DMSetFromOptions(da));
70 PetscCall(DMSetUp(da));
71 PetscCall(DMDASetFieldName(da, 0, "Heat equation: u"));
72 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
73 dt = 1.0 / (ctx.kappa * Mx * Mx);
74
75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76 Extract global vectors from DMDA; then duplicate for remaining
77 vectors that are the same types
78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79 PetscCall(DMCreateGlobalVector(da, &x));
80 PetscCall(VecDuplicate(x, &r));
81
82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 Create timestepping solver context
84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
86 PetscCall(TSSetDM(ts, da));
87 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
88 PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &ctx));
89
90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 Customize nonlinear solver
92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 PetscCall(TSSetType(ts, TSCN));
94
95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96 Set initial conditions
97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98 PetscCall(FormInitialSolution(da, x));
99 PetscCall(TSSetTimeStep(ts, dt));
100 PetscCall(TSSetMaxTime(ts, .02));
101 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE));
102 PetscCall(TSSetSolution(ts, x));
103
104 if (mymonitor) {
105 ctx.ports = NULL;
106 PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
107 }
108
109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110 Set runtime options
111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112 PetscCall(TSSetFromOptions(ts));
113
114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115 Solve nonlinear system
116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117 PetscCall(TSSolve(ts, x));
118 PetscCall(TSGetStepNumber(ts, &steps));
119 PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg));
120 if (flg) {
121 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_WRITE, &viewer));
122 PetscCall(VecView(x, viewer));
123 PetscCall(PetscViewerDestroy(&viewer));
124 }
125
126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 Free work space. All PETSc objects should be destroyed when they
128 are no longer needed.
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130 PetscCall(VecDestroy(&x));
131 PetscCall(VecDestroy(&r));
132 PetscCall(TSDestroy(&ts));
133 PetscCall(DMDestroy(&da));
134
135 PetscCall(PetscFinalize());
136 return 0;
137 }
138 /* ------------------------------------------------------------------- */
139 /*
140 FormFunction - Evaluates nonlinear function, F(x).
141
142 Input Parameters:
143 . ts - the TS context
144 . X - input vector
145 . ptr - optional user-defined context, as set by SNESSetFunction()
146
147 Output Parameter:
148 . F - function vector
149 */
FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,PetscCtx Ctx)150 PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, PetscCtx Ctx)
151 {
152 DM da;
153 PetscInt i, Mx, xs, xm;
154 PetscReal hx, sx;
155 PetscScalar *x, *f;
156 Vec localX;
157 AppCtx *ctx = (AppCtx *)Ctx;
158
159 PetscFunctionBegin;
160 PetscCall(TSGetDM(ts, &da));
161 PetscCall(DMGetLocalVector(da, &localX));
162 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));
163
164 hx = 1.0 / (PetscReal)Mx;
165 sx = 1.0 / (hx * hx);
166
167 /*
168 Scatter ghost points to local vector,using the 2-step process
169 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
170 By placing code between these two statements, computations can be
171 done while messages are in transition.
172 */
173 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
174 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
175
176 /*
177 Get pointers to vector data
178 */
179 PetscCall(DMDAVecGetArrayRead(da, localX, &x));
180 PetscCall(DMDAVecGetArray(da, F, &f));
181
182 /*
183 Get local grid boundaries
184 */
185 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
186
187 /*
188 Compute function over the locally owned part of the grid
189 */
190 for (i = xs; i < xs + xm; i++) {
191 f[i] = ctx->kappa * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
192 if (ctx->allencahn) f[i] += (x[i] - x[i] * x[i] * x[i]);
193 }
194
195 /*
196 Restore vectors
197 */
198 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
199 PetscCall(DMDAVecRestoreArray(da, F, &f));
200 PetscCall(DMRestoreLocalVector(da, &localX));
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
204 /* ------------------------------------------------------------------- */
FormInitialSolution(DM da,Vec U)205 PetscErrorCode FormInitialSolution(DM da, Vec U)
206 {
207 PetscInt i, xs, xm, Mx, scale = 1, N;
208 PetscScalar *u;
209 const PetscScalar *f;
210 PetscReal hx, x, r;
211 Vec finesolution;
212 PetscViewer viewer;
213 PetscBool flg;
214
215 PetscFunctionBegin;
216 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));
217
218 hx = 1.0 / (PetscReal)Mx;
219
220 /*
221 Get pointers to vector data
222 */
223 PetscCall(DMDAVecGetArray(da, U, &u));
224
225 /*
226 Get local grid boundaries
227 */
228 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
229
230 /* InitialSolution is obtained with
231 ./heat -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 9 -ts_max_time 1.e-4 -ts_time_step .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25 -ts_max_steps 15
232 */
233 PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg));
234 if (!flg) {
235 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer));
236 PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution));
237 PetscCall(VecLoad(finesolution, viewer));
238 PetscCall(PetscViewerDestroy(&viewer));
239 PetscCall(VecGetSize(finesolution, &N));
240 scale = N / Mx;
241 PetscCall(VecGetArrayRead(finesolution, &f));
242 }
243
244 /*
245 Compute function over the locally owned part of the grid
246 */
247 for (i = xs; i < xs + xm; i++) {
248 x = i * hx;
249 r = PetscSqrtReal((x - .5) * (x - .5));
250 if (r < .125) u[i] = 1.0;
251 else u[i] = -.5;
252
253 /* With the initial condition above the method is first order in space */
254 /* this is a smooth initial condition so the method becomes second order in space */
255 /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
256 /* u[i] = f[scale*i];*/
257 if (!flg) u[i] = f[scale * i];
258 }
259 if (!flg) {
260 PetscCall(VecRestoreArrayRead(finesolution, &f));
261 PetscCall(VecDestroy(&finesolution));
262 }
263
264 /*
265 Restore vectors
266 */
267 PetscCall(DMDAVecRestoreArray(da, U, &u));
268 PetscFunctionReturn(PETSC_SUCCESS);
269 }
270
271 /*
272 This routine is not parallel
273 */
MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,PetscCtx Ctx)274 PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, PetscCtx Ctx)
275 {
276 AppCtx *ctx = (AppCtx *)Ctx;
277 PetscDrawLG lg;
278 PetscScalar *u;
279 PetscInt Mx, i, xs, xm, cnt;
280 PetscReal x, y, hx, pause, sx, len, max, xx[2], yy[2];
281 PetscDraw draw;
282 Vec localU;
283 DM da;
284 int colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE};
285 const char *const legend[] = {"-kappa (\\grad u,\\grad u)", "(1 - u^2)^2"};
286 PetscDrawAxis axis;
287 PetscDrawViewPorts *ports;
288 PetscReal vbounds[] = {-1.1, 1.1};
289
290 PetscFunctionBegin;
291 PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds));
292 PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1200, 800));
293 PetscCall(TSGetDM(ts, &da));
294 PetscCall(DMGetLocalVector(da, &localU));
295 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));
296 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
297 hx = 1.0 / (PetscReal)Mx;
298 sx = 1.0 / (hx * hx);
299 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
300 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
301 PetscCall(DMDAVecGetArrayRead(da, localU, &u));
302
303 PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg));
304 PetscCall(PetscDrawLGGetDraw(lg, &draw));
305 PetscCall(PetscDrawCheckResizedWindow(draw));
306 if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports));
307 ports = ctx->ports;
308 PetscCall(PetscDrawLGGetAxis(lg, &axis));
309 PetscCall(PetscDrawLGReset(lg));
310
311 xx[0] = 0.0;
312 xx[1] = 1.0;
313 cnt = 2;
314 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL));
315 xs = xx[0] / hx;
316 xm = (xx[1] - xx[0]) / hx;
317
318 /*
319 Plot the energies
320 */
321 PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->allencahn ? 1 : 0)));
322 PetscCall(PetscDrawLGSetColors(lg, colors + 1));
323 PetscCall(PetscDrawViewPortsSet(ports, 2));
324 x = hx * xs;
325 for (i = xs; i < xs + xm; i++) {
326 xx[0] = xx[1] = x;
327 yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
328 if (ctx->allencahn) yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
329 PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
330 x += hx;
331 }
332 PetscCall(PetscDrawGetPause(draw, &pause));
333 PetscCall(PetscDrawSetPause(draw, 0.0));
334 PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", ""));
335 PetscCall(PetscDrawLGSetLegend(lg, legend));
336 PetscCall(PetscDrawLGDraw(lg));
337
338 /*
339 Plot the forces
340 */
341 PetscCall(PetscDrawViewPortsSet(ports, 1));
342 PetscCall(PetscDrawLGReset(lg));
343 x = xs * hx;
344 max = 0.;
345 for (i = xs; i < xs + xm; i++) {
346 xx[0] = xx[1] = x;
347 yy[0] = PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
348 max = PetscMax(max, PetscAbs(yy[0]));
349 if (ctx->allencahn) {
350 yy[1] = PetscRealPart(u[i] - u[i] * u[i] * u[i]);
351 max = PetscMax(max, PetscAbs(yy[1]));
352 }
353 PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
354 x += hx;
355 }
356 PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", ""));
357 PetscCall(PetscDrawLGSetLegend(lg, NULL));
358 PetscCall(PetscDrawLGDraw(lg));
359
360 /*
361 Plot the solution
362 */
363 PetscCall(PetscDrawLGSetDimension(lg, 1));
364 PetscCall(PetscDrawViewPortsSet(ports, 0));
365 PetscCall(PetscDrawLGReset(lg));
366 x = hx * xs;
367 PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1));
368 PetscCall(PetscDrawLGSetColors(lg, colors));
369 for (i = xs; i < xs + xm; i++) {
370 xx[0] = x;
371 yy[0] = PetscRealPart(u[i]);
372 PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
373 x += hx;
374 }
375 PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", ""));
376 PetscCall(PetscDrawLGDraw(lg));
377
378 /*
379 Print the forces as arrows on the solution
380 */
381 x = hx * xs;
382 cnt = xm / 60;
383 cnt = (!cnt) ? 1 : cnt;
384
385 for (i = xs; i < xs + xm; i += cnt) {
386 y = PetscRealPart(u[i]);
387 len = .5 * PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
388 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED));
389 if (ctx->allencahn) {
390 len = .5 * PetscRealPart(u[i] - u[i] * u[i] * u[i]) / max;
391 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE));
392 }
393 x += cnt * hx;
394 }
395 PetscCall(DMDAVecRestoreArrayRead(da, localU, &x));
396 PetscCall(DMRestoreLocalVector(da, &localU));
397 PetscCall(PetscDrawStringSetSize(draw, .2, .2));
398 PetscCall(PetscDrawFlush(draw));
399 PetscCall(PetscDrawSetPause(draw, pause));
400 PetscCall(PetscDrawPause(draw));
401 PetscFunctionReturn(PETSC_SUCCESS);
402 }
403
MyDestroy(PetscCtxRt Ctx)404 PetscErrorCode MyDestroy(PetscCtxRt Ctx)
405 {
406 AppCtx *ctx = *(AppCtx **)Ctx;
407
408 PetscFunctionBegin;
409 PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
410 PetscFunctionReturn(PETSC_SUCCESS);
411 }
412
413 /*TEST
414
415 test:
416 args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial
417
418 test:
419 suffix: 2
420 args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001
421 requires: x
422
423 TEST*/
424