Actual source code: ex12.c

petsc-master 2017-10-17
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*);

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

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 47:      Initialize program
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 49:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

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

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

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

 75:   TSSetMaxTime(ts,1.0);
 76:   TSMonitorSet(ts,MyTSMonitor,0,0);

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

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

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

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

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

115:   PetscFinalize();
116:   return ierr;
117: }
118: /* ------------------------------------------------------------------- */
119: /*
120:    FormFunction - Evaluates nonlinear function, F(x).

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

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

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

143:   hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
144:   hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
145:   /*hxdhy  = hx/hy;*/
146:   /*hydhx  = hy/hx;*/

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

157:   /*
158:      Get pointers to vector data
159:   */
160:   DMDAVecGetArrayDOF(da,localX,&x);
161:   DMDAVecGetArrayDOF(da,F,&f);

163:   /*
164:      Get local grid boundaries
165:   */
166:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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

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

197: /* ------------------------------------------------------------------- */
198: PetscErrorCode FormInitialSolution(DM da,Vec U)
199: {
201:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
202:   PetscScalar    ***u;
203:   PetscReal      hx,hy,x,y,r;

206:   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);

208:   hx = 1.0/(PetscReal)(Mx-1);
209:   hy = 1.0/(PetscReal)(My-1);

211:   /*
212:      Get pointers to vector data
213:   */
214:   DMDAVecGetArrayDOF(da,U,&u);

216:   /*
217:      Get local grid boundaries
218:   */
219:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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

239:   /*
240:      Restore vectors
241:   */
242:   DMDAVecRestoreArrayDOF(da,U,&u);
243:   return(0);
244: }

246: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
247: {
249:   PetscReal      norm;
250:   MPI_Comm       comm;

253:   VecNorm(v,NORM_2,&norm);
254:   PetscObjectGetComm((PetscObject)ts,&comm);
255:   if (step > -1) { /* -1 is used to indicate an interpolated value */
256:     PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
257:   }
258:   return(0);
259: }

261: /*
262:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
263:    Input Parameters:
264:      snes - the SNES context
265:      its - iteration number
266:      fnorm - 2-norm function value (may be estimated)
267:      ctx - optional user-defined context for private data for the
268:          monitor routine, as set by SNESMonitorSet()
269:  */
270: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
271: {

275:   SNESMonitorDefaultShort(snes,its,fnorm,vf);
276:   return(0);
277: }