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