Actual source code: ex13.c
petsc-3.3-p7 2013-05-11
2: static char help[] = "Time-dependent PDE in 2d. Simplified from ex7.c for illustrating how to use TS on a structured domain. \n";
3: /*
4: u_t = uxx + uyy
5: 0 < x < 1, 0 < y < 1;
6: At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125
7: u(x,y) = 0.0 if r >= .125
9: Program usage:
10: mpiexec -n <procs> ./ex13 [-help] [all PETSc options]
11: e.g., mpiexec -n 2 ./ex13 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -use_coloring -snes_monitor -ksp_monitor
12: ./ex13 -use_coloring -drawcontours
13: ./ex13 -use_coloring -drawcontours -draw_pause -1
14: mpiexec -n 2 ./ex13 -drawcontours -ts_type sundials -ts_sundials_monitor_steps
15: */
17: #include <petscdmda.h>
18: #include <petscts.h>
20: /*
21: User-defined data structures and routines
22: */
23: typedef struct {
24: PetscBool drawcontours; /* flag - 1 indicates drawing contours */
25: } MonitorCtx;
26: typedef struct {
27: PetscReal c;
28: } AppCtx;
30: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
31: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
32: extern PetscErrorCode FormInitialSolution(DM,Vec,void*);
33: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
37: int main(int argc,char **argv)
38: {
39: TS ts; /* nonlinear solver */
40: Vec u,r; /* solution, residual vector */
41: Mat J; /* Jacobian matrix */
42: PetscInt steps,maxsteps = 1000; /* iterations for convergence */
44: DM da;
45: ISColoring iscoloring;
46: PetscReal ftime,dt;
47: MonitorCtx usermonitor; /* user-defined monitor context */
48: AppCtx user; /* user-defined work context */
49: PetscBool coloring=PETSC_FALSE;
50: MatFDColoring matfdcoloring;
52: PetscInitialize(&argc,&argv,(char *)0,help);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create distributed array (DMDA) to manage parallel grid and vectors
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
57: 1,1,PETSC_NULL,PETSC_NULL,&da);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Extract global vectors from DMDA;
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62: DMCreateGlobalVector(da,&u);
63: VecDuplicate(u,&r);
65: /* Initialize user application context */
66: user.c = -30.0;
68: usermonitor.drawcontours = PETSC_FALSE;
69: PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Create timestepping solver context
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: TSCreate(PETSC_COMM_WORLD,&ts);
75: TSSetDM(ts,da);
76: TSSetType(ts,TSBEULER);
77: TSSetRHSFunction(ts,r,RHSFunction,&user);
79: /* Set Jacobian */
80: DMCreateMatrix(da,MATAIJ,&J);
81: TSSetRHSJacobian(ts,J,J,RHSJacobian,PETSC_NULL);
83: /* Use coloring to compute Jacobian efficiently */
84: PetscOptionsGetBool(PETSC_NULL,"-use_coloring",&coloring,PETSC_NULL);
85: if (coloring){
86: SNES snes;
87: DMCreateColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
88: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
89: MatFDColoringSetFromOptions(matfdcoloring);
90: ISColoringDestroy(&iscoloring);
92: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
93: TSGetSNES(ts,&snes);
94: SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,SNESDefaultComputeJacobianColor,matfdcoloring);
95: }
97: ftime = 1.0;
98: TSSetDuration(ts,maxsteps,ftime);
99: TSMonitorSet(ts,MyTSMonitor,&usermonitor,PETSC_NULL);
100:
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Set initial conditions
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: FormInitialSolution(da,u,&user);
105: dt = .01;
106: TSSetInitialTimeStep(ts,0.0,dt);
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Set runtime options
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: TSSetFromOptions(ts);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Solve nonlinear system
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: TSSolve(ts,u,&ftime);
117: TSGetTimeStepNumber(ts,&steps);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Free work space.
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: MatDestroy(&J);
123: if (coloring){
124: MatFDColoringDestroy(&matfdcoloring);
125: }
126: VecDestroy(&u);
127: VecDestroy(&r);
128: TSDestroy(&ts);
129: DMDestroy(&da);
131: PetscFinalize();
132: return(0);
133: }
134: /* ------------------------------------------------------------------- */
137: /*
138: RHSFunction - Evaluates nonlinear function, F(u).
140: Input Parameters:
141: . ts - the TS context
142: . U - input vector
143: . ptr - optional user-defined context, as set by TSSetFunction()
145: Output Parameter:
146: . F - function vector
147: */
148: PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
149: {
150: /* PETSC_UNUSED AppCtx *user=(AppCtx*)ptr; */
151: DM da;
153: PetscInt i,j,Mx,My,xs,ys,xm,ym;
154: PetscReal two = 2.0,hx,hy,sx,sy;
155: PetscScalar u,uxx,uyy,**uarray,**f;
156: Vec localU;
159: TSGetDM(ts,&da);
160: DMGetLocalVector(da,&localU);
161: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
162: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
164: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
165: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
167: /*
168: Scatter ghost points to local vector,using the 2-step process
169: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
170: By placing code between these two statements, computations can be
171: done while messages are in transition.
172: */
173: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
174: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
176: /* Get pointers to vector data */
177: DMDAVecGetArray(da,localU,&uarray);
178: DMDAVecGetArray(da,F,&f);
180: /* Get local grid boundaries */
181: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
183: /* Compute function over the locally owned part of the grid */
184: for (j=ys; j<ys+ym; j++) {
185: for (i=xs; i<xs+xm; i++) {
186: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
187: f[j][i] = uarray[j][i];
188: continue;
189: }
190: u = uarray[j][i];
191: uxx = (-two*u + uarray[j][i-1] + uarray[j][i+1])*sx;
192: uyy = (-two*u + uarray[j-1][i] + uarray[j+1][i])*sy;
193: f[j][i] = uxx + uyy;
194: }
195: }
197: /* Restore vectors */
198: DMDAVecRestoreArray(da,localU,&uarray);
199: DMDAVecRestoreArray(da,F,&f);
200: DMRestoreLocalVector(da,&localU);
201: PetscLogFlops(11.0*ym*xm);
202: return(0);
203: }
205: /* --------------------------------------------------------------------- */
208: /*
209: RHSJacobian - User-provided routine to compute the Jacobian of
210: the nonlinear right-hand-side function of the ODE.
212: Input Parameters:
213: ts - the TS context
214: t - current time
215: U - global input vector
216: dummy - optional user-defined context, as set by TSetRHSJacobian()
218: Output Parameters:
219: J - Jacobian matrix
220: Jpre - optionally different preconditioning matrix
221: str - flag indicating matrix structure
222: */
223: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat *J,Mat *Jpre,MatStructure *str,void *ctx)
224: {
226: DM da;
227: DMDALocalInfo info;
228: PetscInt i,j;
229: PetscReal hx,hy,sx,sy;
232: TSGetDM(ts,&da);
233: DMDAGetLocalInfo(da,&info);
234: hx = 1.0/(PetscReal)(info.mx-1); sx = 1.0/(hx*hx);
235: hy = 1.0/(PetscReal)(info.my-1); sy = 1.0/(hy*hy);
236: for (j=info.ys; j<info.ys+info.ym; j++) {
237: for (i=info.xs; i<info.xs+info.xm; i++) {
238: PetscInt nc = 0;
239: MatStencil row,col[5];
240: PetscScalar val[5];
241: row.i = i; row.j = j;
242: if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
243: col[nc].i = i; col[nc].j = j; val[nc++] = 1.0;
244: } else {
245: col[nc].i = i-1; col[nc].j = j; val[nc++] = sx;
246: col[nc].i = i+1; col[nc].j = j; val[nc++] = sx;
247: col[nc].i = i; col[nc].j = j-1; val[nc++] = sy;
248: col[nc].i = i; col[nc].j = j+1; val[nc++] = sy;
249: col[nc].i = i; col[nc].j = j; val[nc++] = -2*sx - 2*sy;
250: }
251: MatSetValuesStencil(*Jpre,1,&row,nc,col,val,INSERT_VALUES);
252: }
253: }
254: MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);
255: MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);
256: if (*J != *Jpre) {
257: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
258: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
259: }
260: return(0);
261: }
263: /* ------------------------------------------------------------------- */
266: PetscErrorCode FormInitialSolution(DM da,Vec U,void* ptr)
267: {
268: AppCtx *user=(AppCtx*)ptr;
269: PetscReal c=user->c;
271: PetscInt i,j,xs,ys,xm,ym,Mx,My;
272: PetscScalar **u;
273: PetscReal hx,hy,x,y,r;
276: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
277: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
279: hx = 1.0/(PetscReal)(Mx-1);
280: hy = 1.0/(PetscReal)(My-1);
282: /* Get pointers to vector data */
283: DMDAVecGetArray(da,U,&u);
285: /* Get local grid boundaries */
286: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
288: /* Compute function over the locally owned part of the grid */
289: for (j=ys; j<ys+ym; j++) {
290: y = j*hy;
291: for (i=xs; i<xs+xm; i++) {
292: x = i*hx;
293: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
294: if (r < .125) {
295: u[j][i] = PetscExpScalar(c*r*r*r);
296: } else {
297: u[j][i] = 0.0;
298: }
299: }
300: }
302: /* Restore vectors */
303: DMDAVecRestoreArray(da,U,&u);
304: return(0);
305: }
309: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ptr)
310: {
312: PetscReal norm;
313: MPI_Comm comm;
314: MonitorCtx *user = (MonitorCtx*)ptr;
317: VecNorm(v,NORM_2,&norm);
318: PetscObjectGetComm((PetscObject)ts,&comm);
319: PetscPrintf(comm,"timestep %D: time %G, solution norm %G\n",step,ptime,norm);
320: if (user->drawcontours){
321: VecView(v,PETSC_VIEWER_DRAW_WORLD);
322: }
323: return(0);
324: }