Actual source code: ex7.c

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

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

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

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

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

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

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

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

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

 69:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

132:   PetscFunctionBeginUser;
133:   PetscCall(TSGetDM(ts, &da));
134:   PetscCall(DMGetLocalVector(da, &localX));
135:   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));

137:   hx = 1.0 / (PetscReal)(Mx - 1);
138:   sx = 1.0 / (hx * hx);
139:   hy = 1.0 / (PetscReal)(My - 1);
140:   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:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
149:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

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

157:   /*
158:      Get local grid boundaries
159:   */
160:   PetscCall(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 + (x[j + 1][i] - x[j - 1][i]) * (x[j + 1][i] - x[j - 1][i]) * .25 * sy);
176:     }
177:   }

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

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

196:   PetscFunctionBeginUser;
197:   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));

199:   hx = 1.0 / (PetscReal)(Mx - 1);
200:   hy = 1.0 / (PetscReal)(My - 1);

202:   /*
203:      Get pointers to vector data
204:   */
205:   PetscCall(DMDAVecGetArray(da, U, &u));

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

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

225:   /*
226:      Restore vectors
227:   */
228:   PetscCall(DMDAVecRestoreArray(da, U, &u));
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

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

237:   PetscFunctionBeginUser;
238:   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* step of -1 indicates an interpolated solution */
239:   PetscCall(VecNorm(v, NORM_2, &norm));
240:   PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
241:   PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
242:   PetscFunctionReturn(PETSC_SUCCESS);
243: }

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

261: /*TEST

263:     test:
264:       args: -ts_max_steps 5

266:     test:
267:       suffix: 2
268:       args: -ts_max_steps 5 -snes_mf_operator

270:     test:
271:       suffix: 3
272:       args: -ts_max_steps 5 -snes_mf -pc_type none

274: TEST*/