Actual source code: ex13.c
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=sqrt((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: DM da;
28: PetscReal c;
29: } AppCtx;
38: int main(int argc,char **argv)
39: {
40: TS ts; /* nonlinear solver */
41: Vec u,r; /* solution, residual vector */
42: Mat J; /* Jacobian matrix */
43: PetscInt steps,maxsteps = 1000; /* iterations for convergence */
45: DM da;
46: ISColoring iscoloring;
47: PetscReal ftime,dt;
48: MonitorCtx usermonitor; /* user-defined monitor context */
49: AppCtx user; /* user-defined work context */
50: PetscBool coloring=PETSC_FALSE;
51: MatFDColoring matfdcoloring;
53: PetscInitialize(&argc,&argv,(char *)0,help);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Create distributed array (DMDA) to manage parallel grid and vectors
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
58: 1,1,PETSC_NULL,PETSC_NULL,&da);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Extract global vectors from DMDA;
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: DMCreateGlobalVector(da,&u);
64: VecDuplicate(u,&r);
66: /* Initialize user application context */
67: user.da = da;
68: user.c = -30.0;
69:
70: usermonitor.drawcontours = PETSC_FALSE;
71: PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create timestepping solver context
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: TSCreate(PETSC_COMM_WORLD,&ts);
77: TSSetType(ts,TSBEULER);
78: TSSetRHSFunction(ts,r,RHSFunction,&user);
80: /* Set Jacobian */
81: DMGetMatrix(da,MATAIJ,&J);
82: TSSetRHSJacobian(ts,J,J,RHSJacobian,PETSC_NULL);
84: /* Use coloring to compute rhs Jacobian efficiently */
85: PetscOptionsGetBool(PETSC_NULL,"-use_coloring",&coloring,PETSC_NULL);
86: if (coloring){
87: DMGetColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
88: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
89: MatFDColoringSetFromOptions(matfdcoloring);
90: ISColoringDestroy(&iscoloring);
92: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))RHSFunction,&user);
93: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
94: }
96: ftime = 1.0;
97: TSSetDuration(ts,maxsteps,ftime);
98: TSMonitorSet(ts,MyTSMonitor,&usermonitor,PETSC_NULL);
99:
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Set initial conditions
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: FormInitialSolution(u,&user);
104: dt = .01;
105: TSSetInitialTimeStep(ts,0.0,dt);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Set runtime options
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: TSSetFromOptions(ts);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Solve nonlinear system
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: TSSolve(ts,u,&ftime);
116: TSGetTimeStepNumber(ts,&steps);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Free work space.
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: MatDestroy(&J);
122: if (coloring){
123: MatFDColoringDestroy(&matfdcoloring);
124: }
125: VecDestroy(&u);
126: VecDestroy(&r);
127: TSDestroy(&ts);
128: DMDestroy(&da);
130: PetscFinalize();
131: return(0);
132: }
133: /* ------------------------------------------------------------------- */
136: /*
137: RHSFunction - Evaluates nonlinear function, F(u).
139: Input Parameters:
140: . ts - the TS context
141: . U - input vector
142: . ptr - optional user-defined context, as set by TSSetFunction()
144: Output Parameter:
145: . F - function vector
146: */
147: PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
148: {
149: AppCtx *user=(AppCtx*)ptr;
150: DM da = (DM)user->da;
152: PetscInt i,j,Mx,My,xs,ys,xm,ym;
153: PetscReal two = 2.0,hx,hy,sx,sy;
154: PetscScalar u,uxx,uyy,**uarray,**f;
155: Vec localU;
158: DMGetLocalVector(da,&localU);
159: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
160: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
162: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
163: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
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,U,INSERT_VALUES,localU);
172: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
174: /* Get pointers to vector data */
175: DMDAVecGetArray(da,localU,&uarray);
176: DMDAVecGetArray(da,F,&f);
178: /* Get local grid boundaries */
179: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
181: /* Compute function over the locally owned part of the grid */
182: for (j=ys; j<ys+ym; j++) {
183: for (i=xs; i<xs+xm; i++) {
184: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
185: f[j][i] = uarray[j][i];
186: continue;
187: }
188: u = uarray[j][i];
189: uxx = (-two*u + uarray[j][i-1] + uarray[j][i+1])*sx;
190: uyy = (-two*u + uarray[j-1][i] + uarray[j+1][i])*sy;
191: f[j][i] = uxx + uyy;
192: }
193: }
195: /* Restore vectors */
196: DMDAVecRestoreArray(da,localU,&uarray);
197: DMDAVecRestoreArray(da,F,&f);
198: DMRestoreLocalVector(da,&localU);
199: PetscLogFlops(11.0*ym*xm);
200: return(0);
201: }
203: /* --------------------------------------------------------------------- */
206: /*
207: RHSJacobian - User-provided routine to compute the Jacobian of
208: the nonlinear right-hand-side function of the ODE.
210: Input Parameters:
211: ts - the TS context
212: t - current time
213: U - global input vector
214: dummy - optional user-defined context, as set by TSetRHSJacobian()
216: Output Parameters:
217: J - Jacobian matrix
218: Jpre - optionally different preconditioning matrix
219: str - flag indicating matrix structure
220: */
221: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat *J,Mat *Jpre,MatStructure *str,void *ctx)
222: {
226: TSDefaultComputeJacobian(ts,t,U,J,Jpre,str,ctx);
227: return(0);
228: }
230: /* ------------------------------------------------------------------- */
233: PetscErrorCode FormInitialSolution(Vec U,void* ptr)
234: {
235: AppCtx *user=(AppCtx*)ptr;
236: DM da=user->da;
237: PetscReal c=user->c;
239: PetscInt i,j,xs,ys,xm,ym,Mx,My;
240: PetscScalar **u;
241: PetscReal hx,hy,x,y,r;
244: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
245: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
247: hx = 1.0/(PetscReal)(Mx-1);
248: hy = 1.0/(PetscReal)(My-1);
250: /* Get pointers to vector data */
251: DMDAVecGetArray(da,U,&u);
253: /* Get local grid boundaries */
254: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
256: /* Compute function over the locally owned part of the grid */
257: for (j=ys; j<ys+ym; j++) {
258: y = j*hy;
259: for (i=xs; i<xs+xm; i++) {
260: x = i*hx;
261: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
262: if (r < .125) {
263: u[j][i] = PetscExpScalar(c*r*r*r);
264: } else {
265: u[j][i] = 0.0;
266: }
267: }
268: }
270: /* Restore vectors */
271: DMDAVecRestoreArray(da,U,&u);
272: return(0);
273: }
277: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ptr)
278: {
280: PetscReal norm;
281: MPI_Comm comm;
282: MonitorCtx *user = (MonitorCtx*)ptr;
285: VecNorm(v,NORM_2,&norm);
286: PetscObjectGetComm((PetscObject)ts,&comm);
287: PetscPrintf(comm,"timestep %D: time %G, solution norm %G\n",step,ptime,norm);
288: if (user->drawcontours){
289: VecView(v,PETSC_VIEWER_DRAW_WORLD);
290: }
291: return(0);
292: }