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 {PetscReal kappa;PetscBool allencahn;PetscDrawViewPorts *ports;} UserCtx;

 38: int main(int argc,char **argv)
 39: {
 40:   TS             ts;                           /* time integrator */
 41:   Vec            x,r;                          /* solution, residual vectors */
 42:   PetscInt       steps,Mx;
 44:   DM             da;
 45:   PetscReal      dt;
 46:   UserCtx        ctx;
 47:   PetscBool      mymonitor;
 48:   PetscViewer    viewer;
 49:   PetscBool      flg;

 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52:      Initialize program
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 54:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 55:   ctx.kappa     = 1.0;
 56:   PetscOptionsGetReal(NULL,NULL,"-kappa",&ctx.kappa,NULL);
 57:   ctx.allencahn = PETSC_FALSE;
 58:   PetscOptionsHasName(NULL,NULL,"-allen-cahn",&ctx.allencahn);
 59:   PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor);

 61:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:      Create distributed array (DMDA) to manage parallel grid and vectors
 63:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 64:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da);
 65:   DMSetFromOptions(da);
 66:   DMSetUp(da);
 67:   DMDASetFieldName(da,0,"Heat equation: u");
 68:   DMDAGetInfo(da,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
 69:   dt   = 1.0/(ctx.kappa*Mx*Mx);

 71:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 72:      Extract global vectors from DMDA; then duplicate for remaining
 73:      vectors that are the same types
 74:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 75:   DMCreateGlobalVector(da,&x);
 76:   VecDuplicate(x,&r);

 78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:      Create timestepping solver context
 80:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 81:   TSCreate(PETSC_COMM_WORLD,&ts);
 82:   TSSetDM(ts,da);
 83:   TSSetProblemType(ts,TS_NONLINEAR);
 84:   TSSetRHSFunction(ts,NULL,FormFunction,&ctx);

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:      Customize nonlinear solver
 88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 89:   TSSetType(ts,TSCN);

 91:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:      Set initial conditions
 93:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 94:   FormInitialSolution(da,x);
 95:   TSSetTimeStep(ts,dt);
 96:   TSSetMaxTime(ts,.02);
 97:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_INTERPOLATE);
 98:   TSSetSolution(ts,x);

100:   if (mymonitor) {
101:     ctx.ports = NULL;
102:     TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy);
103:   }

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Set runtime options
107:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108:   TSSetFromOptions(ts);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Solve nonlinear system
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   TSSolve(ts,x);
114:   TSGetStepNumber(ts,&steps);
115:   PetscOptionsHasName(NULL,NULL,"-square_initial",&flg);
116:   if (flg) {
117:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_WRITE,&viewer);
118:     VecView(x,viewer);
119:     PetscViewerDestroy(&viewer);
120:   }

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:      Free work space.  All PETSc objects should be destroyed when they
124:      are no longer needed.
125:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126:   VecDestroy(&x);
127:   VecDestroy(&r);
128:   TSDestroy(&ts);
129:   DMDestroy(&da);

131:   PetscFinalize();
132:   return ierr;
133: }
134: /* ------------------------------------------------------------------- */
135: /*
136:    FormFunction - Evaluates nonlinear function, F(x).

138:    Input Parameters:
139: .  ts - the TS context
140: .  X - input vector
141: .  ptr - optional user-defined context, as set by SNESSetFunction()

143:    Output Parameter:
144: .  F - function vector
145:  */
146: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
147: {
148:   DM             da;
150:   PetscInt       i,Mx,xs,xm;
151:   PetscReal      hx,sx;
152:   PetscScalar    *x,*f;
153:   Vec            localX;
154:   UserCtx        *ctx = (UserCtx*)ptr;

157:   TSGetDM(ts,&da);
158:   DMGetLocalVector(da,&localX);
159:   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);

161:   hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);

163:   /*
164:      Scatter ghost points to local vector,using the 2-step process
165:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
166:      By placing code between these two statements, computations can be
167:      done while messages are in transition.
168:   */
169:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
170:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

172:   /*
173:      Get pointers to vector data
174:   */
175:   DMDAVecGetArrayRead(da,localX,&x);
176:   DMDAVecGetArray(da,F,&f);

178:   /*
179:      Get local grid boundaries
180:   */
181:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

183:   /*
184:      Compute function over the locally owned part of the grid
185:   */
186:   for (i=xs; i<xs+xm; i++) {
187:     f[i] = ctx->kappa*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
188:     if (ctx->allencahn) f[i] += (x[i] - x[i]*x[i]*x[i]);
189:   }

191:   /*
192:      Restore vectors
193:   */
194:   DMDAVecRestoreArrayRead(da,localX,&x);
195:   DMDAVecRestoreArray(da,F,&f);
196:   DMRestoreLocalVector(da,&localX);
197:   return(0);
198: }

200: /* ------------------------------------------------------------------- */
201: PetscErrorCode FormInitialSolution(DM da,Vec U)
202: {
203:   PetscErrorCode    ierr;
204:   PetscInt          i,xs,xm,Mx,scale=1,N;
205:   PetscScalar       *u;
206:   const PetscScalar *f;
207:   PetscReal         hx,x,r;
208:   Vec               finesolution;
209:   PetscViewer       viewer;
210:   PetscBool         flg;

213:   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);

215:   hx = 1.0/(PetscReal)Mx;

217:   /*
218:      Get pointers to vector data
219:   */
220:   DMDAVecGetArray(da,U,&u);

222:   /*
223:      Get local grid boundaries
224:   */
225:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

227:   /*  InitialSolution is obtained with
228:       ./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
229:   */
230:   PetscOptionsHasName(NULL,NULL,"-square_initial",&flg);
231:   if (!flg) {
232:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer);
233:     VecCreate(PETSC_COMM_WORLD,&finesolution);
234:     VecLoad(finesolution,viewer);
235:     PetscViewerDestroy(&viewer);
236:     VecGetSize(finesolution,&N);
237:     scale = N/Mx;
238:     VecGetArrayRead(finesolution,&f);
239:   }

241:   /*
242:      Compute function over the locally owned part of the grid
243:   */
244:   for (i=xs; i<xs+xm; i++) {
245:     x = i*hx;
246:     r = PetscSqrtScalar((x-.5)*(x-.5));
247:     if (r < .125) u[i] = 1.0;
248:     else u[i] = -.5;

250:     /* With the initial condition above the method is first order in space */
251:     /* this is a smooth initial condition so the method becomes second order in space */
252:     /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
253:     /*  u[i] = f[scale*i];*/
254:     if (!flg) u[i] = f[scale*i];
255:   }
256:   if (!flg) {
257:     VecRestoreArrayRead(finesolution,&f);
258:     VecDestroy(&finesolution);
259:   }

261:   /*
262:      Restore vectors
263:   */
264:   DMDAVecRestoreArray(da,U,&u);
265:   return(0);
266: }

268: /*
269:     This routine is not parallel
270: */
271: PetscErrorCode  MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,void *ptr)
272: {
273:   UserCtx            *ctx = (UserCtx*)ptr;
274:   PetscDrawLG        lg;
275:   PetscErrorCode     ierr;
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};

289:   PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds);
290:   PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1200,800);
291:   TSGetDM(ts,&da);
292:   DMGetLocalVector(da,&localU);
293:   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);
294:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
295:   hx   = 1.0/(PetscReal)Mx; 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) {
304:     PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports);
305:   }
306:   ports = ctx->ports;
307:   PetscDrawLGGetAxis(lg,&axis);
308:   PetscDrawLGReset(lg);

310:   xx[0] = 0.0; xx[1] = 1.0; cnt = 2;
311:   PetscOptionsGetRealArray(NULL,NULL,"-zoom",xx,&cnt,NULL);
312:   xs    = xx[0]/hx; xm = (xx[1] - xx[0])/hx;

314:   /*
315:       Plot the  energies
316:   */
317:   PetscDrawLGSetDimension(lg,1 + (ctx->allencahn ? 1 : 0));
318:   PetscDrawLGSetColors(lg,colors+1);
319:   PetscDrawViewPortsSet(ports,2);
320:   x    = hx*xs;
321:   for (i=xs; i<xs+xm; i++) {
322:     xx[0] = xx[1] = x;
323:     yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx);
324:     if (ctx->allencahn) yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i]));
325:     PetscDrawLGAddPoint(lg,xx,yy);
326:     x   += hx;
327:   }
328:   PetscDrawGetPause(draw,&pause);
329:   PetscDrawSetPause(draw,0.0);
330:   PetscDrawAxisSetLabels(axis,"Energy","","");
331:   PetscDrawLGSetLegend(lg,legend);
332:   PetscDrawLGDraw(lg);

334:   /*
335:       Plot the  forces
336:   */
337:   PetscDrawViewPortsSet(ports,1);
338:   PetscDrawLGReset(lg);
339:   x    = xs*hx;
340:   max  = 0.;
341:   for (i=xs; i<xs+xm; i++) {
342:     xx[0] = xx[1] = x;
343:     yy[0] = PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
344:     max   = PetscMax(max,PetscAbs(yy[0]));
345:     if (ctx->allencahn) {
346:       yy[1] = PetscRealPart(u[i] - u[i]*u[i]*u[i]);
347:       max   = PetscMax(max,PetscAbs(yy[1]));
348:     }
349:     PetscDrawLGAddPoint(lg,xx,yy);
350:     x   += hx;
351:   }
352:   PetscDrawAxisSetLabels(axis,"Right hand side","","");
353:   PetscDrawLGSetLegend(lg,NULL);
354:   PetscDrawLGDraw(lg);

356:   /*
357:         Plot the solution
358:   */
359:   PetscDrawLGSetDimension(lg,1);
360:   PetscDrawViewPortsSet(ports,0);
361:   PetscDrawLGReset(lg);
362:   x    = hx*xs;
363:   PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1);
364:   PetscDrawLGSetColors(lg,colors);
365:   for (i=xs; i<xs+xm; i++) {
366:     xx[0] = x;
367:     yy[0] = PetscRealPart(u[i]);
368:     PetscDrawLGAddPoint(lg,xx,yy);
369:     x    += hx;
370:   }
371:   PetscDrawAxisSetLabels(axis,"Solution","","");
372:   PetscDrawLGDraw(lg);

374:   /*
375:       Print the  forces as arrows on the solution
376:   */
377:   x   = hx*xs;
378:   cnt = xm/60;
379:   cnt = (!cnt) ? 1 : cnt;

381:   for (i=xs; i<xs+xm; i += cnt) {
382:     y    = PetscRealPart(u[i]);
383:     len  = .5*PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
384:     PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED);
385:     if (ctx->allencahn) {
386:       len  = .5*PetscRealPart(u[i] - u[i]*u[i]*u[i])/max;
387:       PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE);
388:     }
389:     x += cnt*hx;
390:   }
391:   DMDAVecRestoreArrayRead(da,localU,&x);
392:   DMRestoreLocalVector(da,&localU);
393:   PetscDrawStringSetSize(draw,.2,.2);
394:   PetscDrawFlush(draw);
395:   PetscDrawSetPause(draw,pause);
396:   PetscDrawPause(draw);
397:   return(0);
398: }

400: PetscErrorCode  MyDestroy(void **ptr)
401: {
402:   UserCtx        *ctx = *(UserCtx**)ptr;

406:   PetscDrawViewPortsDestroy(ctx->ports);
407:   return(0);
408: }

410: /*TEST

412:    test:
413:      args: -ts_monitor -snes_monitor  -pc_type lu  -snes_converged_reason   -ts_type cn  -da_refine 2 -square_initial

415:    test:
416:      suffix: 2
417:      args: -ts_monitor -snes_monitor  -pc_type lu  -snes_converged_reason   -ts_type cn  -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001
418:      requires: x

420: TEST*/