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;
43: DM da;
44: PetscReal dt;
45: UserCtx ctx;
46: PetscBool mymonitor;
47: PetscViewer viewer;
48: PetscBool flg;
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Initialize program
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: PetscInitialize(&argc,&argv,(char*)0,help);
54: ctx.kappa = 1.0;
55: PetscOptionsGetReal(NULL,NULL,"-kappa",&ctx.kappa,NULL);
56: ctx.allencahn = PETSC_FALSE;
57: PetscOptionsHasName(NULL,NULL,"-allen-cahn",&ctx.allencahn);
58: PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create distributed array (DMDA) to manage parallel grid and vectors
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da);
64: DMSetFromOptions(da);
65: DMSetUp(da);
66: DMDASetFieldName(da,0,"Heat equation: u");
67: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
68: dt = 1.0/(ctx.kappa*Mx*Mx);
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Extract global vectors from DMDA; then duplicate for remaining
72: vectors that are the same types
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: DMCreateGlobalVector(da,&x);
75: VecDuplicate(x,&r);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create timestepping solver context
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: TSCreate(PETSC_COMM_WORLD,&ts);
81: TSSetDM(ts,da);
82: TSSetProblemType(ts,TS_NONLINEAR);
83: TSSetRHSFunction(ts,NULL,FormFunction,&ctx);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Customize nonlinear solver
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: TSSetType(ts,TSCN);
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Set initial conditions
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: FormInitialSolution(da,x);
94: TSSetTimeStep(ts,dt);
95: TSSetMaxTime(ts,.02);
96: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_INTERPOLATE);
97: TSSetSolution(ts,x);
99: if (mymonitor) {
100: ctx.ports = NULL;
101: TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy);
102: }
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Set runtime options
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: TSSetFromOptions(ts);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Solve nonlinear system
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: TSSolve(ts,x);
113: TSGetStepNumber(ts,&steps);
114: PetscOptionsHasName(NULL,NULL,"-square_initial",&flg);
115: if (flg) {
116: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_WRITE,&viewer);
117: VecView(x,viewer);
118: PetscViewerDestroy(&viewer);
119: }
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Free work space. All PETSc objects should be destroyed when they
123: are no longer needed.
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: VecDestroy(&x);
126: VecDestroy(&r);
127: TSDestroy(&ts);
128: DMDestroy(&da);
130: PetscFinalize();
131: return 0;
132: }
133: /* ------------------------------------------------------------------- */
134: /*
135: FormFunction - Evaluates nonlinear function, F(x).
137: Input Parameters:
138: . ts - the TS context
139: . X - input vector
140: . ptr - optional user-defined context, as set by SNESSetFunction()
142: Output Parameter:
143: . F - function vector
144: */
145: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
146: {
147: DM da;
148: PetscInt i,Mx,xs,xm;
149: PetscReal hx,sx;
150: PetscScalar *x,*f;
151: Vec localX;
152: UserCtx *ctx = (UserCtx*)ptr;
154: TSGetDM(ts,&da);
155: DMGetLocalVector(da,&localX);
156: 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);
158: hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);
160: /*
161: Scatter ghost points to local vector,using the 2-step process
162: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
163: By placing code between these two statements, computations can be
164: done while messages are in transition.
165: */
166: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
167: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
169: /*
170: Get pointers to vector data
171: */
172: DMDAVecGetArrayRead(da,localX,&x);
173: DMDAVecGetArray(da,F,&f);
175: /*
176: Get local grid boundaries
177: */
178: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
180: /*
181: Compute function over the locally owned part of the grid
182: */
183: for (i=xs; i<xs+xm; i++) {
184: f[i] = ctx->kappa*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
185: if (ctx->allencahn) f[i] += (x[i] - x[i]*x[i]*x[i]);
186: }
188: /*
189: Restore vectors
190: */
191: DMDAVecRestoreArrayRead(da,localX,&x);
192: DMDAVecRestoreArray(da,F,&f);
193: DMRestoreLocalVector(da,&localX);
194: return 0;
195: }
197: /* ------------------------------------------------------------------- */
198: PetscErrorCode FormInitialSolution(DM da,Vec U)
199: {
200: PetscInt i,xs,xm,Mx,scale=1,N;
201: PetscScalar *u;
202: const PetscScalar *f;
203: PetscReal hx,x,r;
204: Vec finesolution;
205: PetscViewer viewer;
206: PetscBool flg;
208: 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);
210: hx = 1.0/(PetscReal)Mx;
212: /*
213: Get pointers to vector data
214: */
215: DMDAVecGetArray(da,U,&u);
217: /*
218: Get local grid boundaries
219: */
220: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
222: /* InitialSolution is obtained with
223: ./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
224: */
225: PetscOptionsHasName(NULL,NULL,"-square_initial",&flg);
226: if (!flg) {
227: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer);
228: VecCreate(PETSC_COMM_WORLD,&finesolution);
229: VecLoad(finesolution,viewer);
230: PetscViewerDestroy(&viewer);
231: VecGetSize(finesolution,&N);
232: scale = N/Mx;
233: VecGetArrayRead(finesolution,&f);
234: }
236: /*
237: Compute function over the locally owned part of the grid
238: */
239: for (i=xs; i<xs+xm; i++) {
240: x = i*hx;
241: r = PetscSqrtScalar((x-.5)*(x-.5));
242: if (r < .125) u[i] = 1.0;
243: else u[i] = -.5;
245: /* With the initial condition above the method is first order in space */
246: /* this is a smooth initial condition so the method becomes second order in space */
247: /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
248: /* u[i] = f[scale*i];*/
249: if (!flg) u[i] = f[scale*i];
250: }
251: if (!flg) {
252: VecRestoreArrayRead(finesolution,&f);
253: VecDestroy(&finesolution);
254: }
256: /*
257: Restore vectors
258: */
259: DMDAVecRestoreArray(da,U,&u);
260: return 0;
261: }
263: /*
264: This routine is not parallel
265: */
266: PetscErrorCode MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,void *ptr)
267: {
268: UserCtx *ctx = (UserCtx*)ptr;
269: PetscDrawLG lg;
270: PetscScalar *u;
271: PetscInt Mx,i,xs,xm,cnt;
272: PetscReal x,y,hx,pause,sx,len,max,xx[2],yy[2];
273: PetscDraw draw;
274: Vec localU;
275: DM da;
276: int colors[] = {PETSC_DRAW_YELLOW,PETSC_DRAW_RED,PETSC_DRAW_BLUE};
277: const char*const legend[] = {"-kappa (\\grad u,\\grad u)","(1 - u^2)^2"};
278: PetscDrawAxis axis;
279: PetscDrawViewPorts *ports;
280: PetscReal vbounds[] = {-1.1,1.1};
282: PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds);
283: PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1200,800);
284: TSGetDM(ts,&da);
285: DMGetLocalVector(da,&localU);
286: 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);
287: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
288: hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);
289: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
290: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
291: DMDAVecGetArrayRead(da,localU,&u);
293: PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,&lg);
294: PetscDrawLGGetDraw(lg,&draw);
295: PetscDrawCheckResizedWindow(draw);
296: if (!ctx->ports) {
297: PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports);
298: }
299: ports = ctx->ports;
300: PetscDrawLGGetAxis(lg,&axis);
301: PetscDrawLGReset(lg);
303: xx[0] = 0.0; xx[1] = 1.0; cnt = 2;
304: PetscOptionsGetRealArray(NULL,NULL,"-zoom",xx,&cnt,NULL);
305: xs = xx[0]/hx; xm = (xx[1] - xx[0])/hx;
307: /*
308: Plot the energies
309: */
310: PetscDrawLGSetDimension(lg,1 + (ctx->allencahn ? 1 : 0));
311: PetscDrawLGSetColors(lg,colors+1);
312: PetscDrawViewPortsSet(ports,2);
313: x = hx*xs;
314: for (i=xs; i<xs+xm; i++) {
315: xx[0] = xx[1] = x;
316: yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx);
317: if (ctx->allencahn) yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i]));
318: PetscDrawLGAddPoint(lg,xx,yy);
319: x += hx;
320: }
321: PetscDrawGetPause(draw,&pause);
322: PetscDrawSetPause(draw,0.0);
323: PetscDrawAxisSetLabels(axis,"Energy","","");
324: PetscDrawLGSetLegend(lg,legend);
325: PetscDrawLGDraw(lg);
327: /*
328: Plot the forces
329: */
330: PetscDrawViewPortsSet(ports,1);
331: PetscDrawLGReset(lg);
332: x = xs*hx;
333: max = 0.;
334: for (i=xs; i<xs+xm; i++) {
335: xx[0] = xx[1] = x;
336: yy[0] = PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
337: max = PetscMax(max,PetscAbs(yy[0]));
338: if (ctx->allencahn) {
339: yy[1] = PetscRealPart(u[i] - u[i]*u[i]*u[i]);
340: max = PetscMax(max,PetscAbs(yy[1]));
341: }
342: PetscDrawLGAddPoint(lg,xx,yy);
343: x += hx;
344: }
345: PetscDrawAxisSetLabels(axis,"Right hand side","","");
346: PetscDrawLGSetLegend(lg,NULL);
347: PetscDrawLGDraw(lg);
349: /*
350: Plot the solution
351: */
352: PetscDrawLGSetDimension(lg,1);
353: PetscDrawViewPortsSet(ports,0);
354: PetscDrawLGReset(lg);
355: x = hx*xs;
356: PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1);
357: PetscDrawLGSetColors(lg,colors);
358: for (i=xs; i<xs+xm; i++) {
359: xx[0] = x;
360: yy[0] = PetscRealPart(u[i]);
361: PetscDrawLGAddPoint(lg,xx,yy);
362: x += hx;
363: }
364: PetscDrawAxisSetLabels(axis,"Solution","","");
365: PetscDrawLGDraw(lg);
367: /*
368: Print the forces as arrows on the solution
369: */
370: x = hx*xs;
371: cnt = xm/60;
372: cnt = (!cnt) ? 1 : cnt;
374: for (i=xs; i<xs+xm; i += cnt) {
375: y = PetscRealPart(u[i]);
376: len = .5*PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
377: PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED);
378: if (ctx->allencahn) {
379: len = .5*PetscRealPart(u[i] - u[i]*u[i]*u[i])/max;
380: PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE);
381: }
382: x += cnt*hx;
383: }
384: DMDAVecRestoreArrayRead(da,localU,&x);
385: DMRestoreLocalVector(da,&localU);
386: PetscDrawStringSetSize(draw,.2,.2);
387: PetscDrawFlush(draw);
388: PetscDrawSetPause(draw,pause);
389: PetscDrawPause(draw);
390: return 0;
391: }
393: PetscErrorCode MyDestroy(void **ptr)
394: {
395: UserCtx *ctx = *(UserCtx**)ptr;
397: PetscDrawViewPortsDestroy(ctx->ports);
398: return 0;
399: }
401: /*TEST
403: test:
404: args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial
406: test:
407: suffix: 2
408: args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001
409: requires: x
411: TEST*/