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*/