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