Actual source code: heat.c
2: static char help[] = "Solves heat equation in 1d.\n";
4: /*
5: Solves the equation
7: u_t = kappa \Delta u
8: Periodic boundary conditions
10: Evolve the heat equation:
11: ---------------
12: ./heat -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor
14: Evolve the Allen-Cahn equation:
15: ---------------
16: ./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
18: Evolve the Allen-Cahn equation: zoom in on part of the domain
19: ---------------
20: ./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
22: 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
23: ./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_dt .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25 -ts_max_steps 15
24: to generate InitialSolution.heat
26: */
27: #include <petscdm.h>
28: #include <petscdmda.h>
29: #include <petscts.h>
30: #include <petscdraw.h>
32: /*
33: User-defined routines
34: */
35: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), MyDestroy(void **);
36: typedef struct {
37: PetscReal kappa;
38: PetscBool allencahn;
39: PetscDrawViewPorts *ports;
40: } UserCtx;
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: UserCtx ctx;
50: PetscBool mymonitor;
51: PetscViewer viewer;
52: PetscBool flg;
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Initialize program
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: PetscFunctionBeginUser;
58: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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));
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);
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));
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));
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Customize nonlinear solver
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: PetscCall(TSSetType(ts, TSCN));
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));
104: if (mymonitor) {
105: ctx.ports = NULL;
106: PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
107: }
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Set runtime options
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: PetscCall(TSSetFromOptions(ts));
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: }
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));
135: PetscCall(PetscFinalize());
136: return 0;
137: }
138: /* ------------------------------------------------------------------- */
139: /*
140: FormFunction - Evaluates nonlinear function, F(x).
142: Input Parameters:
143: . ts - the TS context
144: . X - input vector
145: . ptr - optional user-defined context, as set by SNESSetFunction()
147: Output Parameter:
148: . F - function vector
149: */
150: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
151: {
152: DM da;
153: PetscInt i, Mx, xs, xm;
154: PetscReal hx, sx;
155: PetscScalar *x, *f;
156: Vec localX;
157: UserCtx *ctx = (UserCtx *)ptr;
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));
164: hx = 1.0 / (PetscReal)Mx;
165: sx = 1.0 / (hx * hx);
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));
176: /*
177: Get pointers to vector data
178: */
179: PetscCall(DMDAVecGetArrayRead(da, localX, &x));
180: PetscCall(DMDAVecGetArray(da, F, &f));
182: /*
183: Get local grid boundaries
184: */
185: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
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: }
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: }
204: /* ------------------------------------------------------------------- */
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;
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));
218: hx = 1.0 / (PetscReal)Mx;
220: /*
221: Get pointers to vector data
222: */
223: PetscCall(DMDAVecGetArray(da, U, &u));
225: /*
226: Get local grid boundaries
227: */
228: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
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_dt .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: }
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;
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: }
264: /*
265: Restore vectors
266: */
267: PetscCall(DMDAVecRestoreArray(da, U, &u));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*
272: This routine is not parallel
273: */
274: PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr)
275: {
276: UserCtx *ctx = (UserCtx *)ptr;
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};
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));
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));
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;
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));
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));
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));
378: /*
379: Print the forces as arrows on the solution
380: */
381: x = hx * xs;
382: cnt = xm / 60;
383: cnt = (!cnt) ? 1 : cnt;
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: }
404: PetscErrorCode MyDestroy(void **ptr)
405: {
406: UserCtx *ctx = *(UserCtx **)ptr;
408: PetscFunctionBegin;
409: PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*TEST
415: test:
416: args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial
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
423: TEST*/