Actual source code: ex12.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  2: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";

  4: /*
  5:   Solves the equation

  7:     u_tt - \Delta u = 0

  9:   which we split into two first-order systems

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

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


 28: /*
 29:    User-defined routines
 30: */
 31: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
 32: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
 33: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,PetscViewerAndFormat*);

 37: int main(int argc,char **argv)
 38: {
 39:   TS                   ts;                         /* nonlinear solver */
 40:   Vec                  x,r;                        /* solution, residual vectors */
 41:   PetscInt             steps,maxsteps = 100;       /* iterations for convergence */
 42:   PetscErrorCode       ierr;
 43:   DM                   da;
 44:   PetscReal            ftime;
 45:   SNES                 ts_snes;
 46:   PetscViewerAndFormat *vf;

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:      Initialize program
 50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 51:   PetscInitialize(&argc,&argv,(char*)0,help);

 53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:      Create distributed array (DMDA) to manage parallel grid and vectors
 55:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 56:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
 57:                       2,1,NULL,NULL,&da);
 58:   DMDASetFieldName(da,0,"u");
 59:   DMDASetFieldName(da,1,"v");

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

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

 76:   TSSetDuration(ts,maxsteps,1.0);
 77:   TSMonitorSet(ts,MyTSMonitor,0,0);

 79:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80:      Customize nonlinear solver
 81:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 82:   TSSetType(ts,TSBEULER);
 83:   TSGetSNES(ts,&ts_snes);
 84:   PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
 85:   SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Set initial conditions
 89:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 90:   FormInitialSolution(da,x);
 91:   TSSetInitialTimeStep(ts,0.0,.0001);
 92:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
 93:   TSSetSolution(ts,x);

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:      Set runtime options
 97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   TSSetFromOptions(ts);

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Solve nonlinear system
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   TSSolve(ts,x);
104:   TSGetSolveTime(ts,&ftime);
105:   TSGetTimeStepNumber(ts,&steps);

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Free work space.  All PETSc objects should be destroyed when they
109:      are no longer needed.
110:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111:   VecDestroy(&x);
112:   VecDestroy(&r);
113:   TSDestroy(&ts);
114:   DMDestroy(&da);

116:   PetscFinalize();
117:   return(0);
118: }
119: /* ------------------------------------------------------------------- */
122: /*
123:    FormFunction - Evaluates nonlinear function, F(x).

125:    Input Parameters:
126: .  ts - the TS context
127: .  X - input vector
128: .  ptr - optional user-defined context, as set by SNESSetFunction()

130:    Output Parameter:
131: .  F - function vector
132:  */
133: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
134: {
135:   DM             da = (DM)ptr;
137:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
138:   PetscReal      hx,hy,/*hxdhy,hydhx,*/ sx,sy;
139:   PetscScalar    u,uxx,uyy,v,***x,***f;
140:   Vec            localX;

143:   DMGetLocalVector(da,&localX);
144:   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);

146:   hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
147:   hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
148:   /*hxdhy  = hx/hy;*/
149:   /*hydhx  = hy/hx;*/

151:   /*
152:      Scatter ghost points to local vector,using the 2-step process
153:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
154:      By placing code between these two statements, computations can be
155:      done while messages are in transition.
156:   */
157:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
158:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

160:   /*
161:      Get pointers to vector data
162:   */
163:   DMDAVecGetArrayDOF(da,localX,&x);
164:   DMDAVecGetArrayDOF(da,F,&f);

166:   /*
167:      Get local grid boundaries
168:   */
169:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

171:   /*
172:      Compute function over the locally owned part of the grid
173:   */
174:   for (j=ys; j<ys+ym; j++) {
175:     for (i=xs; i<xs+xm; i++) {
176:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
177:         f[j][i][0] = x[j][i][0];
178:         f[j][i][1] = x[j][i][1];
179:         continue;
180:       }
181:       u          = x[j][i][0];
182:       v          = x[j][i][1];
183:       uxx        = (-2.0*u + x[j][i-1][0] + x[j][i+1][0])*sx;
184:       uyy        = (-2.0*u + x[j-1][i][0] + x[j+1][i][0])*sy;
185:       f[j][i][0] = v;
186:       f[j][i][1] = uxx + uyy;
187:     }
188:   }

190:   /*
191:      Restore vectors
192:   */
193:   DMDAVecRestoreArrayDOF(da,localX,&x);
194:   DMDAVecRestoreArrayDOF(da,F,&f);
195:   DMRestoreLocalVector(da,&localX);
196:   PetscLogFlops(11.0*ym*xm);
197:   return(0);
198: }

200: /* ------------------------------------------------------------------- */
203: PetscErrorCode FormInitialSolution(DM da,Vec U)
204: {
206:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
207:   PetscScalar    ***u;
208:   PetscReal      hx,hy,x,y,r;

211:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
212:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

214:   hx = 1.0/(PetscReal)(Mx-1);
215:   hy = 1.0/(PetscReal)(My-1);

217:   /*
218:      Get pointers to vector data
219:   */
220:   DMDAVecGetArrayDOF(da,U,&u);

222:   /*
223:      Get local grid boundaries
224:   */
225:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

227:   /*
228:      Compute function over the locally owned part of the grid
229:   */
230:   for (j=ys; j<ys+ym; j++) {
231:     y = j*hy;
232:     for (i=xs; i<xs+xm; i++) {
233:       x = i*hx;
234:       r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
235:       if (r < .125) {
236:         u[j][i][0] = PetscExpReal(-30.0*r*r*r);
237:         u[j][i][1] = 0.0;
238:       } else {
239:         u[j][i][0] = 0.0;
240:         u[j][i][1] = 0.0;
241:       }
242:     }
243:   }

245:   /*
246:      Restore vectors
247:   */
248:   DMDAVecRestoreArrayDOF(da,U,&u);
249:   return(0);
250: }

254: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
255: {
257:   PetscReal      norm;
258:   MPI_Comm       comm;

261:   VecNorm(v,NORM_2,&norm);
262:   PetscObjectGetComm((PetscObject)ts,&comm);
263:   if (step > -1) { /* -1 is used to indicate an interpolated value */
264:     PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
265:   }
266:   return(0);
267: }

271: /*
272:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
273:    Input Parameters:
274:      snes - the SNES context
275:      its - iteration number
276:      fnorm - 2-norm function value (may be estimated)
277:      ctx - optional user-defined context for private data for the
278:          monitor routine, as set by SNESMonitorSet()
279:  */
280: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
281: {

285:   SNESMonitorDefaultShort(snes,its,fnorm,vf);
286:   return(0);
287: }