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: }