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