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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: PetscInitialize(&argc, &argv, (char *)0, help);
59: ctx.kappa = 1.0;
60: PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL);
61: ctx.allencahn = PETSC_FALSE;
62: PetscOptionsHasName(NULL, NULL, "-allen-cahn", &ctx.allencahn);
63: PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor);
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Create distributed array (DMDA) to manage parallel grid and vectors
67: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da);
69: DMSetFromOptions(da);
70: DMSetUp(da);
71: DMDASetFieldName(da, 0, "Heat equation: u");
72: 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: DMCreateGlobalVector(da, &x);
80: VecDuplicate(x, &r);
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Create timestepping solver context
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: TSCreate(PETSC_COMM_WORLD, &ts);
86: TSSetDM(ts, da);
87: TSSetProblemType(ts, TS_NONLINEAR);
88: TSSetRHSFunction(ts, NULL, FormFunction, &ctx);
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Customize nonlinear solver
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: TSSetType(ts, TSCN);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Set initial conditions
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: FormInitialSolution(da, x);
99: TSSetTimeStep(ts, dt);
100: TSSetMaxTime(ts, .02);
101: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE);
102: TSSetSolution(ts, x);
104: if (mymonitor) {
105: ctx.ports = NULL;
106: TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy);
107: }
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Set runtime options
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: TSSetFromOptions(ts);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Solve nonlinear system
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: TSSolve(ts, x);
118: TSGetStepNumber(ts, &steps);
119: PetscOptionsHasName(NULL, NULL, "-square_initial", &flg);
120: if (flg) {
121: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_WRITE, &viewer);
122: VecView(x, viewer);
123: PetscViewerDestroy(&viewer);
124: }
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Free work space. All PETSc objects should be destroyed when they
128: are no longer needed.
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: VecDestroy(&x);
131: VecDestroy(&r);
132: TSDestroy(&ts);
133: DMDestroy(&da);
135: 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: TSGetDM(ts, &da);
160: DMGetLocalVector(da, &localX);
161: 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: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
173: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
175: /*
176: Get pointers to vector data
177: */
178: DMDAVecGetArrayRead(da, localX, &x);
179: DMDAVecGetArray(da, F, &f);
181: /*
182: Get local grid boundaries
183: */
184: 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: DMDAVecRestoreArrayRead(da, localX, &x);
198: DMDAVecRestoreArray(da, F, &f);
199: DMRestoreLocalVector(da, &localX);
200: return 0;
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: 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);
216: hx = 1.0 / (PetscReal)Mx;
218: /*
219: Get pointers to vector data
220: */
221: DMDAVecGetArray(da, U, &u);
223: /*
224: Get local grid boundaries
225: */
226: DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
228: /* InitialSolution is obtained with
229: ./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
230: */
231: PetscOptionsHasName(NULL, NULL, "-square_initial", &flg);
232: if (!flg) {
233: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer);
234: VecCreate(PETSC_COMM_WORLD, &finesolution);
235: VecLoad(finesolution, viewer);
236: PetscViewerDestroy(&viewer);
237: VecGetSize(finesolution, &N);
238: scale = N / Mx;
239: VecGetArrayRead(finesolution, &f);
240: }
242: /*
243: Compute function over the locally owned part of the grid
244: */
245: for (i = xs; i < xs + xm; i++) {
246: x = i * hx;
247: r = PetscSqrtScalar((x - .5) * (x - .5));
248: if (r < .125) u[i] = 1.0;
249: else u[i] = -.5;
251: /* With the initial condition above the method is first order in space */
252: /* this is a smooth initial condition so the method becomes second order in space */
253: /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
254: /* u[i] = f[scale*i];*/
255: if (!flg) u[i] = f[scale * i];
256: }
257: if (!flg) {
258: VecRestoreArrayRead(finesolution, &f);
259: VecDestroy(&finesolution);
260: }
262: /*
263: Restore vectors
264: */
265: DMDAVecRestoreArray(da, U, &u);
266: return 0;
267: }
269: /*
270: This routine is not parallel
271: */
272: PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr)
273: {
274: UserCtx *ctx = (UserCtx *)ptr;
275: PetscDrawLG lg;
276: PetscScalar *u;
277: PetscInt Mx, i, xs, xm, cnt;
278: PetscReal x, y, hx, pause, sx, len, max, xx[2], yy[2];
279: PetscDraw draw;
280: Vec localU;
281: DM da;
282: int colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE};
283: const char *const legend[] = {"-kappa (\\grad u,\\grad u)", "(1 - u^2)^2"};
284: PetscDrawAxis axis;
285: PetscDrawViewPorts *ports;
286: PetscReal vbounds[] = {-1.1, 1.1};
288: PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds);
289: PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1200, 800);
290: TSGetDM(ts, &da);
291: DMGetLocalVector(da, &localU);
292: 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);
293: DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
294: hx = 1.0 / (PetscReal)Mx;
295: sx = 1.0 / (hx * hx);
296: DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
297: DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
298: DMDAVecGetArrayRead(da, localU, &u);
300: PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg);
301: PetscDrawLGGetDraw(lg, &draw);
302: PetscDrawCheckResizedWindow(draw);
303: if (!ctx->ports) PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports);
304: ports = ctx->ports;
305: PetscDrawLGGetAxis(lg, &axis);
306: PetscDrawLGReset(lg);
308: xx[0] = 0.0;
309: xx[1] = 1.0;
310: cnt = 2;
311: PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL);
312: xs = xx[0] / hx;
313: xm = (xx[1] - xx[0]) / hx;
315: /*
316: Plot the energies
317: */
318: PetscDrawLGSetDimension(lg, 1 + (ctx->allencahn ? 1 : 0));
319: PetscDrawLGSetColors(lg, colors + 1);
320: PetscDrawViewPortsSet(ports, 2);
321: x = hx * xs;
322: for (i = xs; i < xs + xm; i++) {
323: xx[0] = xx[1] = x;
324: yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
325: if (ctx->allencahn) yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
326: PetscDrawLGAddPoint(lg, xx, yy);
327: x += hx;
328: }
329: PetscDrawGetPause(draw, &pause);
330: PetscDrawSetPause(draw, 0.0);
331: PetscDrawAxisSetLabels(axis, "Energy", "", "");
332: PetscDrawLGSetLegend(lg, legend);
333: PetscDrawLGDraw(lg);
335: /*
336: Plot the forces
337: */
338: PetscDrawViewPortsSet(ports, 1);
339: PetscDrawLGReset(lg);
340: x = xs * hx;
341: max = 0.;
342: for (i = xs; i < xs + xm; i++) {
343: xx[0] = xx[1] = x;
344: yy[0] = PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
345: max = PetscMax(max, PetscAbs(yy[0]));
346: if (ctx->allencahn) {
347: yy[1] = PetscRealPart(u[i] - u[i] * u[i] * u[i]);
348: max = PetscMax(max, PetscAbs(yy[1]));
349: }
350: PetscDrawLGAddPoint(lg, xx, yy);
351: x += hx;
352: }
353: PetscDrawAxisSetLabels(axis, "Right hand side", "", "");
354: PetscDrawLGSetLegend(lg, NULL);
355: PetscDrawLGDraw(lg);
357: /*
358: Plot the solution
359: */
360: PetscDrawLGSetDimension(lg, 1);
361: PetscDrawViewPortsSet(ports, 0);
362: PetscDrawLGReset(lg);
363: x = hx * xs;
364: PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1);
365: PetscDrawLGSetColors(lg, colors);
366: for (i = xs; i < xs + xm; i++) {
367: xx[0] = x;
368: yy[0] = PetscRealPart(u[i]);
369: PetscDrawLGAddPoint(lg, xx, yy);
370: x += hx;
371: }
372: PetscDrawAxisSetLabels(axis, "Solution", "", "");
373: PetscDrawLGDraw(lg);
375: /*
376: Print the forces as arrows on the solution
377: */
378: x = hx * xs;
379: cnt = xm / 60;
380: cnt = (!cnt) ? 1 : cnt;
382: for (i = xs; i < xs + xm; i += cnt) {
383: y = PetscRealPart(u[i]);
384: len = .5 * PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
385: PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED);
386: if (ctx->allencahn) {
387: len = .5 * PetscRealPart(u[i] - u[i] * u[i] * u[i]) / max;
388: PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE);
389: }
390: x += cnt * hx;
391: }
392: DMDAVecRestoreArrayRead(da, localU, &x);
393: DMRestoreLocalVector(da, &localU);
394: PetscDrawStringSetSize(draw, .2, .2);
395: PetscDrawFlush(draw);
396: PetscDrawSetPause(draw, pause);
397: PetscDrawPause(draw);
398: return 0;
399: }
401: PetscErrorCode MyDestroy(void **ptr)
402: {
403: UserCtx *ctx = *(UserCtx **)ptr;
405: PetscDrawViewPortsDestroy(ctx->ports);
406: return 0;
407: }
409: /*TEST
411: test:
412: args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial
414: test:
415: suffix: 2
416: args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001
417: requires: x
419: TEST*/