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:   DM                    da;
 38:   PetscReal             ftime;
 39:   SNES                  ts_snes;
 40:   PetscBool             usemonitor = PETSC_TRUE;
 41:   PetscViewerAndFormat *vf;

 43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 44:      Initialize program
 45:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 46:   PetscFunctionBeginUser;
 47:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 48:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-usemonitor", &usemonitor, NULL));

 50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51:      Create distributed array (DMDA) to manage parallel grid and vectors
 52:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 53:   PetscCall(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:   PetscCall(DMSetFromOptions(da));
 55:   PetscCall(DMSetUp(da));
 56:   PetscCall(DMDASetFieldName(da, 0, "u"));
 57:   PetscCall(DMDASetFieldName(da, 1, "v"));

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

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

 74:   PetscCall(TSSetMaxTime(ts, 1.0));
 75:   if (usemonitor) PetscCall(TSMonitorSet(ts, MyTSMonitor, 0, 0));

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

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

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

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

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

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

128:    Output Parameter:
129: .  F - function vector
130:  */
131: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
132: {
133:   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;

139:   PetscFunctionBeginUser;
140:   PetscCall(DMGetLocalVector(da, &localX));
141:   PetscCall(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);
144:   sx = 1.0 / (hx * hx);
145:   hy = 1.0 / (PetscReal)(My - 1);
146:   sy = 1.0 / (hy * hy);
147:   /*hxdhy  = hx/hy;*/
148:   /*hydhx  = hy/hx;*/

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

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

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

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

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

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

206:   PetscFunctionBeginUser;
207:   PetscCall(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));

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

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

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

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

240:   /*
241:      Restore vectors
242:   */
243:   PetscCall(DMDAVecRestoreArrayDOF(da, U, &u));
244:   PetscFunctionReturn(PETSC_SUCCESS);
245: }

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

252:   PetscFunctionBeginUser;
253:   PetscCall(VecNorm(v, NORM_2, &norm));
254:   PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
255:   if (step > -1) { /* -1 is used to indicate an interpolated value */
256:     PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
257:   }
258:   PetscFunctionReturn(PETSC_SUCCESS);
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: {
272:   PetscFunctionBeginUser;
273:   PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf));
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }
276: /*TEST

278:     test:
279:       args: -da_grid_x 20 -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor_short
280:       requires: !single

282:     test:
283:       suffix: 2
284:       args: -da_grid_x 20 -ts_max_time 0.11 -ts_dt 1e-1 -ts_type glle -ts_monitor -ksp_monitor_short
285:       requires: !single

287:     test:
288:       suffix: glvis_da_2d_vect
289:       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
290:       requires: !single

292:     test:
293:       suffix: glvis_da_2d_vect_ll
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 -viewer_glvis_dm_da_ll
295:       requires: !single

297: TEST*/