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