Actual source code: heat.c

petsc-3.4.2 2013-07-02
  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  -wait   -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  -wait   -ts_type cn  -da_refine 5   -allen-cahn -kappa .001 -ts_final_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_final_time 5  -zoom .25,.45 -wait  -mymonitor


 23: 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
 24: ./heat -square_initial -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason    -ts_type cn  -da_refine 9 -ts_final_time 1.e-4 -ts_dt .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25  -ts_max_steps 15
 25: to generate binaryoutput then do mv binaryoutput InitialSolution.heat to obtain the initial solution file

 27: */
 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;

 40: int main(int argc,char **argv)
 41: {
 42:   TS             ts;                           /* time integrator */
 43:   Vec            x,r;                          /* solution, residual vectors */
 44:   Mat            J;                            /* Jacobian matrix */
 45:   PetscInt       steps,Mx, maxsteps = 1000000;
 47:   DM             da;
 48:   PetscReal      dt;
 49:   PetscReal      vbounds[] = {-1.1,1.1};
 50:   PetscBool      wait;
 51:   UserCtx        ctx;
 52:   PetscBool      mymonitor;

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

 64:   PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds);
 65:   PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1200,800);

 67:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68:      Create distributed array (DMDA) to manage parallel grid and vectors
 69:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 70:   DMDACreate1d(PETSC_COMM_WORLD, DMDA_BOUNDARY_PERIODIC, -10,1,2,NULL,&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:   TSSetInitialTimeStep(ts,0.0,dt);
100:   TSSetDuration(ts,maxsteps,.02);
101:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_INTERPOLATE);
102:   TSSetSolution(ts,x);


105:   if (mymonitor) {
106:     ctx.ports = NULL;
107:     TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy);
108:   }

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Set runtime options
112:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   TSSetFromOptions(ts);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Solve nonlinear system
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118:   TSSolve(ts,x);
119:   wait = PETSC_FALSE;
120:   PetscOptionsGetBool(NULL,"-wait",&wait,NULL);
121:   if (wait) {
122:     PetscSleep(-1);
123:   }
124:   TSGetTimeStepNumber(ts,&steps);
125:   VecView(x,PETSC_VIEWER_BINARY_WORLD);

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Free work space.  All PETSc objects should be destroyed when they
129:      are no longer needed.
130:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131:   MatDestroy(&J);
132:   VecDestroy(&x);
133:   VecDestroy(&r);
134:   TSDestroy(&ts);
135:   DMDestroy(&da);

137:   PetscFinalize();
138:   return(0);
139: }
140: /* ------------------------------------------------------------------- */
143: /*
144:    FormFunction - Evaluates nonlinear function, F(x).

146:    Input Parameters:
147: .  ts - the TS context
148: .  X - input vector
149: .  ptr - optional user-defined context, as set by SNESSetFunction()

151:    Output Parameter:
152: .  F - function vector
153:  */
154: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
155: {
156:   DM             da;
158:   PetscInt       i,Mx,xs,xm;
159:   PetscReal      hx,sx;
160:   PetscScalar    *x,*f;
161:   Vec            localX;
162:   UserCtx        *ctx = (UserCtx*)ptr;

165:   TSGetDM(ts,&da);
166:   DMGetLocalVector(da,&localX);
167:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
168:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

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

172:   /*
173:      Scatter ghost points to local vector,using the 2-step process
174:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
175:      By placing code between these two statements, computations can be
176:      done while messages are in transition.
177:   */
178:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
179:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

181:   /*
182:      Get pointers to vector data
183:   */
184:   DMDAVecGetArray(da,localX,&x);
185:   DMDAVecGetArray(da,F,&f);

187:   /*
188:      Get local grid boundaries
189:   */
190:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

192:   /*
193:      Compute function over the locally owned part of the grid
194:   */
195:   for (i=xs; i<xs+xm; i++) {
196:     f[i] = ctx->kappa*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
197:     if (ctx->allencahn) f[i] += (x[i] - x[i]*x[i]*x[i]);
198:   }

200:   /*
201:      Restore vectors
202:   */
203:   DMDAVecRestoreArray(da,localX,&x);
204:   DMDAVecRestoreArray(da,F,&f);
205:   DMRestoreLocalVector(da,&localX);
206:   return(0);
207: }

209: /* ------------------------------------------------------------------- */
212: PetscErrorCode FormInitialSolution(DM da,Vec U)
213: {
214:   PetscErrorCode    ierr;
215:   PetscInt          i,xs,xm,Mx,scale,N;
216:   PetscScalar       *u;
217:   const PetscScalar *f;
218:   PetscReal         hx,x,r;
219:   Vec               finesolution;
220:   PetscViewer       viewer;
221:   PetscBool         flg;

224:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
225:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

227:   hx = 1.0/(PetscReal)Mx;

229:   /*
230:      Get pointers to vector data
231:   */
232:   DMDAVecGetArray(da,U,&u);

234:   /*
235:      Get local grid boundaries
236:   */
237:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

239:   /*  InitialSolution is obtained with
240:       ./heat -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason    -ts_type cn  -da_refine 9 -ts_final_time 1.e-4 -ts_dt .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25  -ts_max_steps 15
241:    */

243:   /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution",FILE_MODE_READ,&viewer); */
244:   /* VecCreate(PETSC_COMM_WORLD,&finesolution); */
245:   /* VecLoad(finesolution,viewer); */
246:   /* PetscViewerDestroy(&viewer); */
247:   /* VecGetSize(finesolution,&N); */
248:   /* scale = N/Mx; */
249:   /* VecGetArrayRead(finesolution,&f); */

251:   PetscOptionsHasName(NULL,"-square_initial",&flg);
252:   if (!flg) {
253:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer);
254:     VecCreate(PETSC_COMM_WORLD,&finesolution);
255:     VecLoad(finesolution,viewer);
256:     PetscViewerDestroy(&viewer);
257:     VecGetSize(finesolution,&N);
258:     scale = N/Mx;
259:     VecGetArrayRead(finesolution,&f);
260:   }

262:   /*
263:      Compute function over the locally owned part of the grid
264:   */
265:   for (i=xs; i<xs+xm; i++) {
266:     x = i*hx;
267:     r = PetscSqrtScalar((x-.5)*(x-.5));
268:     if (r < .125) u[i] = 1.0;
269:     else u[i] = -.5;

271:     /* With the initial condition above the method is first order in space */
272:     /* this is a smooth initial condition so the method becomes second order in space */
273:     /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
274:     /*  u[i] = f[scale*i];*/
275:     if (!flg) u[i] = f[scale*i];
276:   }
277:   if (!flg) {
278:     VecRestoreArrayRead(finesolution,&f);
279:     VecDestroy(&finesolution);
280:   }
281:   /* VecRestoreArrayRead(finesolution,&f);
282:    VecDestroy(&finesolution);*/

284:   /*
285:      Restore vectors
286:   */
287:   DMDAVecRestoreArray(da,U,&u);
288:   return(0);
289: }

293: /*
294:     This routine is not parallel
295: */
296: PetscErrorCode  MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,void *ptr)
297: {
298:   UserCtx            *ctx = (UserCtx*)ptr;
299:   PetscDrawLG        lg;
300:   PetscErrorCode     ierr;
301:   PetscScalar        *u;
302:   PetscInt           Mx,i,xs,xm,cnt;
303:   PetscReal          x,y,hx,pause,sx,len,max,xx[2],yy[2];
304:   PetscDraw          draw;
305:   Vec                localU;
306:   DM                 da;
307:   int                colors[] = {PETSC_DRAW_YELLOW,PETSC_DRAW_RED,PETSC_DRAW_BLUE};
308:   const char*const   legend[] = {"-kappa (\\grad u,\\grad u)","(1 - u^2)^2"};
309:   PetscDrawAxis      axis;
310:   PetscDrawViewPorts *ports;


314:   TSGetDM(ts,&da);
315:   DMGetLocalVector(da,&localU);
316:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
317:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
318:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
319:   hx   = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);
320:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
321:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
322:   DMDAVecGetArray(da,localU,&u);

324:   PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,&lg);
325:   PetscDrawLGGetDraw(lg,&draw);
326:   PetscDrawCheckResizedWindow(draw);
327:   if (!ctx->ports) {
328:     PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports);
329:   }
330:   ports = ctx->ports;
331:   PetscDrawLGGetAxis(lg,&axis);
332:   PetscDrawLGReset(lg);

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

338:   /*
339:       Plot the  energies
340:   */
341:   PetscDrawLGSetDimension(lg,1 + (ctx->allencahn ? 1 : 0));
342:   PetscDrawLGSetColors(lg,colors+1);
343:   PetscDrawViewPortsSet(ports,2);
344:   x    = hx*xs;
345:   for (i=xs; i<xs+xm; i++) {
346:     xx[0] = xx[1] = x;
347:     yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx);
348:     if (ctx->allencahn) yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i]));
349:     PetscDrawLGAddPoint(lg,xx,yy);
350:     x   += hx;
351:   }
352:   PetscDrawGetPause(draw,&pause);
353:   PetscDrawSetPause(draw,0.0);
354:   PetscDrawAxisSetLabels(axis,"Energy","","");
355:   PetscDrawLGSetLegend(lg,legend);
356:   PetscDrawLGDraw(lg);

358:   /*
359:       Plot the  forces
360:   */
361:   PetscDrawViewPortsSet(ports,1);
362:   PetscDrawLGReset(lg);
363:   x    = xs*hx;;
364:   max  = 0.;
365:   for (i=xs; i<xs+xm; i++) {
366:     xx[0] = xx[1] = x;
367:     yy[0] = PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
368:     max   = PetscMax(max,PetscAbs(yy[0]));
369:     if (ctx->allencahn) {
370:       yy[1] = PetscRealPart(u[i] - u[i]*u[i]*u[i]);
371:       max   = PetscMax(max,PetscAbs(yy[1]));
372:     }
373:     PetscDrawLGAddPoint(lg,xx,yy);
374:     x   += hx;
375:   }
376:   PetscDrawAxisSetLabels(axis,"Right hand side","","");
377:   PetscDrawLGSetLegend(lg,NULL);
378:   PetscDrawLGDraw(lg);

380:   /*
381:         Plot the solution
382:   */
383:   PetscDrawLGSetDimension(lg,1);
384:   PetscDrawViewPortsSet(ports,0);
385:   PetscDrawLGReset(lg);
386:   x    = hx*xs;
387:   PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1);
388:   PetscDrawLGSetColors(lg,colors);
389:   for (i=xs; i<xs+xm; i++) {
390:     xx[0] = x;
391:     yy[0] = PetscRealPart(u[i]);
392:     PetscDrawLGAddPoint(lg,xx,yy);
393:     x    += hx;
394:   }
395:   PetscDrawAxisSetLabels(axis,"Solution","","");
396:   PetscDrawLGDraw(lg);

398:   /*
399:       Print the  forces as arrows on the solution
400:   */
401:   x   = hx*xs;
402:   cnt = xm/60;
403:   cnt = (!cnt) ? 1 : cnt;

405:   for (i=xs; i<xs+xm; i += cnt) {
406:     y    = PetscRealPart(u[i]);
407:     len  = .5*PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
408:     PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED);
409:     if (ctx->allencahn) {
410:       len  = .5*PetscRealPart(u[i] - u[i]*u[i]*u[i])/max;
411:       PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE);
412:     }
413:     x += cnt*hx;
414:   }
415:   DMDAVecRestoreArray(da,localU,&x);
416:   DMRestoreLocalVector(da,&localU);
417:   PetscDrawStringSetSize(draw,.2,.2);
418:   PetscDrawFlush(draw);
419:   PetscDrawSetPause(draw,pause);
420:   PetscDrawPause(draw);
421:   return(0);
422: }

426: PetscErrorCode  MyDestroy(void **ptr)
427: {
428:   UserCtx        *ctx = *(UserCtx**)ptr;

432:   PetscDrawViewPortsDestroy(ctx->ports);
433:   return(0);
434: }