Actual source code: ex12.c

  2: /* Program usage:  mpiexec -n <procs> ex5 [-help] [all PETSc options] */

  4: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";

  6: /*
  7:   Solves the equation

  9:     u_tt - \Delta u = 0

 11:   which we split into two first-order systems

 13:     u_t -     v    = 0    so that     F(u,v).u = -v
 14:     v_t - \Delta u = 0    so that     F(u,v).v = -Delta u

 16:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 17:    Include "petscts.h" so that we can use SNES solvers.  Note that this
 18:    file automatically includes:
 19:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 20:      petscmat.h - matrices
 21:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 22:      petscviewer.h - viewers               petscpc.h  - preconditioners
 23:      petscksp.h   - linear solvers
 24: */
 25: #include <petscdmda.h>
 26: #include <petscts.h>


 29: /* 
 30:    User-defined routines
 31: */

 38: int main(int argc,char **argv)
 39: {
 40:   TS                     ts;                 /* nonlinear solver */
 41:   Vec                    x,r;                  /* solution, residual vectors */
 42:   Mat                    J;                    /* Jacobian matrix */
 43:   PetscInt               steps,maxsteps = 100;     /* iterations for convergence */
 44:   PetscErrorCode         ierr;
 45:   DM                     da;
 46:   MatFDColoring          matfdcoloring;
 47:   ISColoring             iscoloring;
 48:   PetscReal              ftime;
 49:   SNES                   ts_snes;

 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52:      Initialize program
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 54:   PetscInitialize(&argc,&argv,(char *)0,help);
 55:   PetscOptionsGetInt(PETSC_NULL,"-max_steps",&maxsteps,PETSC_NULL);

 57:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 58:      Create distributed array (DMDA) to manage parallel grid and vectors
 59:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 60:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
 61:                     2,1,PETSC_NULL,PETSC_NULL,&da);

 63:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 64:      Extract global vectors from DMDA; then duplicate for remaining
 65:      vectors that are the same types
 66:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 67:   DMCreateGlobalVector(da,&x);
 68:   VecDuplicate(x,&r);

 70:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71:      Create timestepping solver context
 72:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 73:   TSCreate(PETSC_COMM_WORLD,&ts);
 74:   TSSetProblemType(ts,TS_NONLINEAR);
 75:   TSSetRHSFunction(ts,PETSC_NULL,FormFunction,da);

 77:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 78:      Create matrix data structure; set Jacobian evaluation routine

 80:      Set Jacobian matrix data structure and default Jacobian evaluation
 81:      routine. User can override with:
 82:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 83:                 (unless user explicitly sets preconditioner) 
 84:      -snes_mf_operator : form preconditioning matrix as set by the user,
 85:                          but use matrix-free approx for Jacobian-vector
 86:                          products within Newton-Krylov method

 88:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 89:   DMGetColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
 90:   DMGetMatrix(da,MATAIJ,&J);
 91:   MatFDColoringCreate(J,iscoloring,&matfdcoloring);
 92:   ISColoringDestroy(&iscoloring);
 93:   MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,da);
 94:   MatFDColoringSetFromOptions(matfdcoloring);
 95:   TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);

 97:   TSSetDuration(ts,maxsteps,1.0);
 98:   TSMonitorSet(ts,MyTSMonitor,0,0);
 99: 
100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Customize nonlinear solver
102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   TSSetType(ts,TSBEULER);
104:   TSGetSNES(ts,&ts_snes);
105:   SNESMonitorSet(ts_snes,MySNESMonitor,PETSC_NULL,PETSC_NULL);
106: 
107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Set initial conditions
109:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110:   FormInitialSolution(da,x);
111:   TSSetInitialTimeStep(ts,0.0,.0001);
112:   TSSetSolution(ts,x);

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Set runtime options
116:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   TSSetFromOptions(ts);

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Solve nonlinear system
121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122:   TSSolve(ts,x,&ftime);
123:   TSGetTimeStepNumber(ts,&steps);

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:      Free work space.  All PETSc objects should be destroyed when they
127:      are no longer needed.
128:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129:   MatDestroy(&J);
130:   MatFDColoringDestroy(&matfdcoloring);
131:   VecDestroy(&x);
132:   VecDestroy(&r);
133:   TSDestroy(&ts);
134:   DMDestroy(&da);

136:   PetscFinalize();
137:   return(0);
138: }
139: /* ------------------------------------------------------------------- */
142: /* 
143:    FormFunction - Evaluates nonlinear function, F(x).

145:    Input Parameters:
146: .  ts - the TS context
147: .  X - input vector
148: .  ptr - optional user-defined context, as set by SNESSetFunction()

150:    Output Parameter:
151: .  F - function vector
152:  */
153: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
154: {
155:   DM             da = (DM)ptr;
157:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
158:   PetscReal      two = 2.0,hx,hy,/*hxdhy,hydhx,*/sx,sy;
159:   PetscScalar    u,uxx,uyy,v,***x,***f;
160:   Vec            localX;

163:   DMGetLocalVector(da,&localX);
164:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
165:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

167:   hx     = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
168:   hy     = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
169:   /*hxdhy  = hx/hy;*/
170:   /*hydhx  = hy/hx;*/

172:   /*
173:      Scatter ghost points to local vector,using the 2-step process
174:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
175:      By placing code between these two statements, computations can be
176:      done while messages are in transition.
177:   */
178:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
179:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

181:   /*
182:      Get pointers to vector data
183:   */
184:   DMDAVecGetArrayDOF(da,localX,&x);
185:   DMDAVecGetArrayDOF(da,F,&f);

187:   /*
188:      Get local grid boundaries
189:   */
190:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

192:   /*
193:      Compute function over the locally owned part of the grid
194:   */
195:   for (j=ys; j<ys+ym; j++) {
196:     for (i=xs; i<xs+xm; i++) {
197:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
198:         f[j][i][0] = x[j][i][0];
199:         f[j][i][1] = x[j][i][1];
200:         continue;
201:       }
202:       u          = x[j][i][0];
203:       v          = x[j][i][1];
204:       uxx        = (two*u - x[j][i-1][0] - x[j][i+1][0])*sx;
205:       uyy        = (two*u - x[j-1][i][0] - x[j+1][i][0])*sy;
206:       f[j][i][0] = -v;
207:       f[j][i][1] = -(uxx + uyy);
208:     }
209:   }

211:   /*
212:      Restore vectors
213:   */
214:   DMDAVecRestoreArrayDOF(da,localX,&x);
215:   DMDAVecRestoreArrayDOF(da,F,&f);
216:   DMRestoreLocalVector(da,&localX);
217:   PetscLogFlops(11.0*ym*xm);
218:   return(0);
219: }

221: /* ------------------------------------------------------------------- */
224: PetscErrorCode FormInitialSolution(DM da,Vec U)
225: {
227:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
228:   PetscScalar    ***u;
229:   PetscReal      hx,hy,x,y,r;

232:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
233:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

235:   hx     = 1.0/(PetscReal)(Mx-1);
236:   hy     = 1.0/(PetscReal)(My-1);

238:   /*
239:      Get pointers to vector data
240:   */
241:   DMDAVecGetArrayDOF(da,U,&u);

243:   /*
244:      Get local grid boundaries
245:   */
246:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

248:   /*
249:      Compute function over the locally owned part of the grid
250:   */
251:   for (j=ys; j<ys+ym; j++) {
252:     y = j*hy;
253:     for (i=xs; i<xs+xm; i++) {
254:       x = i*hx;
255:       r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
256:       if (r < .125) {
257:         u[j][i][0] = PetscExpScalar(-30.0*r*r*r);
258:         u[j][i][1] = 0.0;
259:       } else {
260:         u[j][i][0] = 0.0;
261:         u[j][i][1] = 0.0;
262:       }
263:     }
264:   }

266:   /*
267:      Restore vectors
268:   */
269:   DMDAVecRestoreArrayDOF(da,U,&u);
270:   return(0);
271: }

275: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
276: {
278:   PetscReal      norm;
279:   MPI_Comm       comm;

282:   VecNorm(v,NORM_2,&norm);
283:   PetscObjectGetComm((PetscObject)ts,&comm);
284:   PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
285:   return(0);
286: }

290: /*
291:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
292:    Input Parameters:
293:      snes - the SNES context
294:      its - iteration number
295:      fnorm - 2-norm function value (may be estimated)
296:      ctx - optional user-defined context for private data for the 
297:          monitor routine, as set by SNESMonitorSet()
298:  */
299: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
300: {
302: 
304:   SNESMonitorDefaultShort(snes,its,fnorm,ctx);
305:   return(0);
306: }