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