Actual source code: ex12.c

  1: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
  2: /*
  3:   Solves the equation

  5:     u_tt - \Delta u = 0

  7:   which we split into two first-order systems

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

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

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

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

 44:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45:      Initialize program
 46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 47:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 48:   PetscOptionsGetBool(NULL,NULL,"-usemonitor",&usemonitor,NULL);

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

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

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

 74:   TSSetMaxTime(ts,1.0);
 75:   if (usemonitor) {
 76:     TSMonitorSet(ts,MyTSMonitor,0,0);
 77:   }

 79:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80:      Customize nonlinear solver
 81:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 82:   TSSetType(ts,TSBEULER);
 83:   TSGetSNES(ts,&ts_snes);
 84:   if (usemonitor) {
 85:     PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
 86:     SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
 87:   }
 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:      Set initial conditions
 90:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 91:   FormInitialSolution(da,x);
 92:   TSSetTimeStep(ts,.0001);
 93:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
 94:   TSSetSolution(ts,x);

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

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Solve nonlinear system
103:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104:   TSSolve(ts,x);
105:   TSGetSolveTime(ts,&ftime);
106:   TSGetStepNumber(ts,&steps);
107:   VecViewFromOptions(x,NULL,"-final_sol");

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

118:   PetscFinalize();
119:   return ierr;
120: }
121: /* ------------------------------------------------------------------- */
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: /* ------------------------------------------------------------------- */
201: PetscErrorCode FormInitialSolution(DM da,Vec U)
202: {
204:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
205:   PetscScalar    ***u;
206:   PetscReal      hx,hy,x,y,r;

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

211:   hx = 1.0/(PetscReal)(Mx-1);
212:   hy = 1.0/(PetscReal)(My-1);

214:   /*
215:      Get pointers to vector data
216:   */
217:   DMDAVecGetArrayDOF(da,U,&u);

219:   /*
220:      Get local grid boundaries
221:   */
222:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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

242:   /*
243:      Restore vectors
244:   */
245:   DMDAVecRestoreArrayDOF(da,U,&u);
246:   return(0);
247: }

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

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

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

278:   SNESMonitorDefaultShort(snes,its,fnorm,vf);
279:   return(0);
280: }
281: /*TEST

283:     test:
284:       args: -da_grid_x 20 -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor_short
285:       requires: !single

287:     test:
288:       suffix: 2
289:       args: -da_grid_x 20 -ts_max_time 0.11 -ts_dt 1e-1 -ts_type glle -ts_monitor -ksp_monitor_short
290:       requires: !single

292:     test:
293:       suffix: glvis_da_2d_vect
294:       args: -usemonitor 0 -da_grid_x 20 -ts_max_time 0.3 -ts_dt 1e-1 -ts_type glle -final_sol glvis: -viewer_glvis_dm_da_bs 2,0
295:       requires: !single

297:     test:
298:       suffix: glvis_da_2d_vect_ll
299:       args: -usemonitor 0 -da_grid_x 20 -ts_max_time 0.3 -ts_dt 1e-1 -ts_type glle -final_sol glvis: -viewer_glvis_dm_da_bs 2,0 -viewer_glvis_dm_da_ll
300:       requires: !single

302: TEST*/