Actual source code: ex7.c


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

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

 18: /*
 19:    User-defined routines
 20: */
 21: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
 22: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
 23: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,PetscViewerAndFormat*);

 25: int main(int argc,char **argv)
 26: {
 27:   TS                   ts;                         /* time integrator */
 28:   SNES                 snes;
 29:   Vec                  x,r;                        /* solution, residual vectors */
 30:   PetscErrorCode       ierr;
 31:   DM                   da;
 32:   PetscViewerAndFormat *vf;

 34:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 35:      Initialize program
 36:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 37:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 38:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 39:      Create distributed array (DMDA) to manage parallel grid and vectors
 40:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 41:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
 42:   DMSetFromOptions(da);
 43:   DMSetUp(da);

 45:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 46:      Extract global vectors from DMDA; then duplicate for remaining
 47:      vectors that are the same types
 48:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 49:   DMCreateGlobalVector(da,&x);
 50:   VecDuplicate(x,&r);

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:      Create timestepping solver context
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 55:   TSCreate(PETSC_COMM_WORLD,&ts);
 56:   TSSetProblemType(ts,TS_NONLINEAR);
 57:   TSSetRHSFunction(ts,NULL,FormFunction,da);

 59:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 60:      Create matrix data structure; set Jacobian evaluation routine

 62:      Set Jacobian matrix data structure and default Jacobian evaluation
 63:      routine. User can override with:
 64:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 65:                 (unless user explicitly sets preconditioner)
 66:      -snes_mf_operator : form preconditioning matrix as set by the user,
 67:                          but use matrix-free approx for Jacobian-vector
 68:                          products within Newton-Krylov method

 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 72:   TSSetMaxTime(ts,1.0);
 73:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
 74:   TSMonitorSet(ts,MyTSMonitor,PETSC_VIEWER_STDOUT_WORLD,NULL);
 75:   TSSetDM(ts,da);
 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:      Customize nonlinear solver
 78:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 79:   TSSetType(ts,TSBEULER);
 80:   TSGetSNES(ts,&snes);
 81:   PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
 82:   SNESMonitorSet(snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void*))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Set initial conditions
 86:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 87:   FormInitialSolution(da,x);
 88:   TSSetTimeStep(ts,.0001);
 89:   TSSetSolution(ts,x);

 91:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:      Set runtime options
 93:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 94:   TSSetFromOptions(ts);

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:      Solve nonlinear system
 98:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 99:   TSSolve(ts,x);

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Free work space.  All PETSc objects should be destroyed when they
103:      are no longer needed.
104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105:   VecDestroy(&x);
106:   VecDestroy(&r);
107:   TSDestroy(&ts);
108:   DMDestroy(&da);

110:   PetscFinalize();
111:   return ierr;
112: }
113: /* ------------------------------------------------------------------- */
114: /*
115:    FormFunction - Evaluates nonlinear function, F(x).

117:    Input Parameters:
118: .  ts - the TS context
119: .  X - input vector
120: .  ptr - optional user-defined context, as set by SNESSetFunction()

122:    Output Parameter:
123: .  F - function vector
124:  */
125: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
126: {
127:   DM             da;
129:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
130:   PetscReal      two = 2.0,hx,hy,sx,sy;
131:   PetscScalar    u,uxx,uyy,**x,**f;
132:   Vec            localX;

135:   TSGetDM(ts,&da);
136:   DMGetLocalVector(da,&localX);
137:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

139:   hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
140:   hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);

142:   /*
143:      Scatter ghost points to local vector,using the 2-step process
144:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
145:      By placing code between these two statements, computations can be
146:      done while messages are in transition.
147:   */
148:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
149:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

151:   /*
152:      Get pointers to vector data
153:   */
154:   DMDAVecGetArrayRead(da,localX,&x);
155:   DMDAVecGetArray(da,F,&f);

157:   /*
158:      Get local grid boundaries
159:   */
160:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

162:   /*
163:      Compute function over the locally owned part of the grid
164:   */
165:   for (j=ys; j<ys+ym; j++) {
166:     for (i=xs; i<xs+xm; i++) {
167:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
168:         f[j][i] = x[j][i];
169:         continue;
170:       }
171:       u   = x[j][i];
172:       uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
173:       uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
174:       /*      f[j][i] = -(uxx + uyy); */
175:       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 +
176:                                               (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
177:     }
178:   }

180:   /*
181:      Restore vectors
182:   */
183:   DMDAVecRestoreArrayRead(da,localX,&x);
184:   DMDAVecRestoreArray(da,F,&f);
185:   DMRestoreLocalVector(da,&localX);
186:   PetscLogFlops(11.0*ym*xm);
187:   return(0);
188: }

190: /* ------------------------------------------------------------------- */
191: PetscErrorCode FormInitialSolution(DM da,Vec U)
192: {
194:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
195:   PetscScalar    **u;
196:   PetscReal      hx,hy,x,y,r;

199:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

201:   hx = 1.0/(PetscReal)(Mx-1);
202:   hy = 1.0/(PetscReal)(My-1);

204:   /*
205:      Get pointers to vector data
206:   */
207:   DMDAVecGetArray(da,U,&u);

209:   /*
210:      Get local grid boundaries
211:   */
212:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

214:   /*
215:      Compute function over the locally owned part of the grid
216:   */
217:   for (j=ys; j<ys+ym; j++) {
218:     y = j*hy;
219:     for (i=xs; i<xs+xm; i++) {
220:       x = i*hx;
221:       r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
222:       if (r < .125) u[j][i] = PetscExpReal(-30.0*r*r*r);
223:       else          u[j][i] = 0.0;
224:     }
225:   }

227:   /*
228:      Restore vectors
229:   */
230:   DMDAVecRestoreArray(da,U,&u);
231:   return(0);
232: }

234: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
235: {
237:   PetscReal      norm;
238:   MPI_Comm       comm;

241:   if (step < 0) return(0); /* step of -1 indicates an interpolated solution */
242:   VecNorm(v,NORM_2,&norm);
243:   PetscObjectGetComm((PetscObject)ts,&comm);
244:   PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
245:   return(0);
246: }

248: /*
249:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
250:    Input Parameters:
251:      snes - the SNES context
252:      its - iteration number
253:      fnorm - 2-norm function value (may be estimated)
254:      ctx - optional user-defined context for private data for the
255:          monitor routine, as set by SNESMonitorSet()
256:  */
257: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
258: {

262:   SNESMonitorDefaultShort(snes,its,fnorm,vf);
263:   return(0);
264: }

266: /*TEST

268:     test:
269:       args: -ts_max_steps 5

271:     test:
272:       suffix: 2
273:       args: -ts_max_steps 5  -snes_mf_operator

275:     test:
276:       suffix: 3
277:       args: -ts_max_steps 5  -snes_mf -pc_type none

279: TEST*/