Actual source code: ex7.c


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

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

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

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

 33:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 34:      Initialize program
 35:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 36:   PetscFunctionBeginUser;
 37:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 38:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 39:      Create distributed array (DMDA) to manage parallel grid and vectors
 40:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 41:   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));
 42:   PetscCall(DMSetFromOptions(da));
 43:   PetscCall(DMSetUp(da));

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

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

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

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

 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

138:   hx = 1.0 / (PetscReal)(Mx - 1);
139:   sx = 1.0 / (hx * hx);
140:   hy = 1.0 / (PetscReal)(My - 1);
141:   sy = 1.0 / (hy * hy);

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

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

158:   /*
159:      Get local grid boundaries
160:   */
161:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

163:   /*
164:      Compute function over the locally owned part of the grid
165:   */
166:   for (j = ys; j < ys + ym; j++) {
167:     for (i = xs; i < xs + xm; i++) {
168:       if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
169:         f[j][i] = x[j][i];
170:         continue;
171:       }
172:       u   = x[j][i];
173:       uxx = (two * u - x[j][i - 1] - x[j][i + 1]) * sx;
174:       uyy = (two * u - x[j - 1][i] - x[j + 1][i]) * sy;
175:       /*      f[j][i] = -(uxx + uyy); */
176:       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);
177:     }
178:   }

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

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

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

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

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

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

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

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

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

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

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

262: /*TEST

264:     test:
265:       args: -ts_max_steps 5

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

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

275: TEST*/