Actual source code: ex7.c

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

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


  7: /* 
  8:    Include "petscda.h" so that we can use distributed arrays (DAs).
  9:    Include "petscts.h" so that we can use SNES solvers.  Note that this
 10:    file automatically includes:
 11:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 12:      petscmat.h - matrices
 13:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 14:      petscviewer.h - viewers               petscpc.h  - preconditioners
 15:      petscksp.h   - linear solvers
 16: */
 17:  #include petscda.h
 18:  #include petscts.h


 21: /* 
 22:    User-defined routines
 23: */

 30: int main(int argc,char **argv)
 31: {
 32:   TS                     ts;                 /* nonlinear solver */
 33:   Vec                    x,r;                  /* solution, residual vectors */
 34:   Mat                    J;                    /* Jacobian matrix */
 35:   PetscInt               steps,maxsteps = 100;     /* iterations for convergence */
 36:   PetscErrorCode         ierr;
 37:   DA                     da;
 38:   MatFDColoring          matfdcoloring;
 39:   ISColoring             iscoloring;
 40:   PetscReal              ftime;
 41:   SNES                   ts_snes;

 43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 44:      Initialize program
 45:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 46:   PetscInitialize(&argc,&argv,(char *)0,help);
 47:   PetscOptionsGetInt(PETSC_NULL,"-max_steps",&maxsteps,PETSC_NULL);

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Create distributed array (DA) to manage parallel grid and vectors
 51:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 52:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,
 53:                     1,1,PETSC_NULL,PETSC_NULL,&da);

 55:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:      Extract global vectors from DA; then duplicate for remaining
 57:      vectors that are the same types
 58:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 59:   DACreateGlobalVector(da,&x);
 60:   VecDuplicate(x,&r);

 62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63:      Create timestepping solver context
 64:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 65:   TSCreate(PETSC_COMM_WORLD,&ts);
 66:   TSSetProblemType(ts,TS_NONLINEAR);
 67:   TSSetRHSFunction(ts,FormFunction,da);

 69:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 70:      Create matrix data structure; set Jacobian evaluation routine

 72:      Set Jacobian matrix data structure and default Jacobian evaluation
 73:      routine. User can override with:
 74:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 75:                 (unless user explicitly sets preconditioner) 
 76:      -snes_mf_operator : form preconditioning matrix as set by the user,
 77:                          but use matrix-free approx for Jacobian-vector
 78:                          products within Newton-Krylov method

 80:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 81:   DAGetColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
 82:   DAGetMatrix(da,MATAIJ,&J);
 83:   MatFDColoringCreate(J,iscoloring,&matfdcoloring);
 84:   ISColoringDestroy(iscoloring);
 85:   MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,da);
 86:   MatFDColoringSetFromOptions(matfdcoloring);
 87:   TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);

 89:   TSSetDuration(ts,maxsteps,1.0);
 90:   TSMonitorSet(ts,MyTSMonitor,0,0);
 91: 
 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Customize nonlinear solver
 94:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   TSSetType(ts,TSBEULER);
 96:   TSGetSNES(ts,&ts_snes);
 97:   SNESMonitorSet(ts_snes,MySNESMonitor,PETSC_NULL,PETSC_NULL);
 98: 
 99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:      Set initial conditions
101:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102:   FormInitialSolution(da,x);
103:   TSSetInitialTimeStep(ts,0.0,.0001);
104:   TSSetSolution(ts,x);

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Set runtime options
108:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109:   TSSetFromOptions(ts);

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:      Solve nonlinear system
113:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114:   TSStep(ts,&steps,&ftime);

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Free work space.  All PETSc objects should be destroyed when they
118:      are no longer needed.
119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   MatDestroy(J);
121:   MatFDColoringDestroy(matfdcoloring);
122:   VecDestroy(x);
123:   VecDestroy(r);
124:   TSDestroy(ts);
125:   DADestroy(da);

127:   PetscFinalize();
128:   return(0);
129: }
130: /* ------------------------------------------------------------------- */
133: /* 
134:    FormFunction - Evaluates nonlinear function, F(x).

136:    Input Parameters:
137: .  ts - the TS context
138: .  X - input vector
139: .  ptr - optional user-defined context, as set by SNESSetFunction()

141:    Output Parameter:
142: .  F - function vector
143:  */
144: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
145: {
146:   DA             da = (DA)ptr;
148:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
149:   PetscReal      two = 2.0,hx,hy,hxdhy,hydhx,sx,sy;
150:   PetscScalar    u,uxx,uyy,**x,**f;
151:   Vec            localX;

154:   DAGetLocalVector(da,&localX);
155:   DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
156:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

158:   hx     = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
159:   hy     = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
160:   hxdhy  = hx/hy;
161:   hydhx  = hy/hx;

163:   /*
164:      Scatter ghost points to local vector,using the 2-step process
165:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
166:      By placing code between these two statements, computations can be
167:      done while messages are in transition.
168:   */
169:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
170:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

172:   /*
173:      Get pointers to vector data
174:   */
175:   DAVecGetArray(da,localX,&x);
176:   DAVecGetArray(da,F,&f);

178:   /*
179:      Get local grid boundaries
180:   */
181:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

183:   /*
184:      Compute function over the locally owned part of the grid
185:   */
186:   for (j=ys; j<ys+ym; j++) {
187:     for (i=xs; i<xs+xm; i++) {
188:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
189:         f[j][i] = x[j][i];
190:         continue;
191:       }
192:       u       = x[j][i];
193:       uxx     = (two*u - x[j][i-1] - x[j][i+1])*sx;
194:       uyy     = (two*u - x[j-1][i] - x[j+1][i])*sy;
195:       /*      f[j][i] = -(uxx + uyy); */
196:       f[j][i] = -u*(uxx + uyy) - (4.0 - 1.0)*((x[j][i+1] - x[j][i-1])*(x[j][i+1] - x[j][i-1])*.25*sx +
197:                                             (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
198:     }
199:   }

201:   /*
202:      Restore vectors
203:   */
204:   DAVecRestoreArray(da,localX,&x);
205:   DAVecRestoreArray(da,F,&f);
206:   DARestoreLocalVector(da,&localX);
207:   PetscLogFlops(11.0*ym*xm);
208:   return(0);
209: }

211: /* ------------------------------------------------------------------- */
214: PetscErrorCode FormInitialSolution(DA da,Vec U)
215: {
217:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
218:   PetscScalar    **u;
219:   PetscReal      hx,hy,x,y,r;

222:   DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
223:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

225:   hx     = 1.0/(PetscReal)(Mx-1);
226:   hy     = 1.0/(PetscReal)(My-1);

228:   /*
229:      Get pointers to vector data
230:   */
231:   DAVecGetArray(da,U,&u);

233:   /*
234:      Get local grid boundaries
235:   */
236:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

238:   /*
239:      Compute function over the locally owned part of the grid
240:   */
241:   for (j=ys; j<ys+ym; j++) {
242:     y = j*hy;
243:     for (i=xs; i<xs+xm; i++) {
244:       x = i*hx;
245:       r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
246:       if (r < .125) {
247:         u[j][i] = PetscExpScalar(-30.0*r*r*r);
248:       } else {
249:         u[j][i] = 0.0;
250:       }
251:     }
252:   }

254:   /*
255:      Restore vectors
256:   */
257:   DAVecRestoreArray(da,U,&u);
258:   return(0);
259: }

263: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
264: {
266:   PetscReal      norm;
267:   MPI_Comm       comm;

270:   VecNorm(v,NORM_2,&norm);
271:   PetscObjectGetComm((PetscObject)ts,&comm);
272:   PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
273:   return(0);
274: }

278: /*
279:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
280:    Input Parameters:
281:      snes - the SNES context
282:      its - iteration number
283:      fnorm - 2-norm function value (may be estimated)
284:      ctx - optional user-defined context for private data for the 
285:          monitor routine, as set by SNESMonitorSet()
286:  */
287: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
288: {
290: 
292:   SNESMonitorDefaultShort(snes,its,fnorm,ctx);
293:   return(0);
294: }