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: }