Actual source code: ex35.c

  1: static const char help[] = "-Laplacian u = b as a nonlinear problem.\n\n";

  3: /*
  4:     The linear and nonlinear versions of these should give almost identical results on this problem

  6:     Richardson
  7:       Nonlinear:
  8:         -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor

 10:       Linear:
 11:         -snes_rtol 1.e-12 -snes_monitor -ksp_rtol 1.e-12  -ksp_monitor -ksp_type richardson -pc_type none -ksp_richardson_self_scale -info

 13:     GMRES
 14:       Nonlinear:
 15:        -snes_rtol 1.e-12 -snes_monitor  -snes_type ngmres

 17:       Linear:
 18:        -snes_rtol 1.e-12 -snes_monitor  -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none

 20:     CG
 21:        Nonlinear:
 22:             -snes_rtol 1.e-12 -snes_monitor  -snes_type ncg -snes_linesearch_monitor

 24:        Linear:
 25:              -snes_rtol 1.e-12 -snes_monitor  -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none

 27:     Multigrid
 28:        Linear:
 29:           1 level:
 30:             -snes_rtol 1.e-12 -snes_monitor  -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor
 31:             -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12  -ksp_monitor_true_residual

 33:           n levels:
 34:             -da_refine n

 36:        Nonlinear:
 37:          1 level:
 38:            -snes_rtol 1.e-12 -snes_monitor  -snes_type fas -fas_levels_snes_monitor

 40:           n levels:
 41:             -da_refine n  -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly
 42: */

 44: /*
 45:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 46:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 47: */
 48: #include <petscdm.h>
 49: #include <petscdmda.h>
 50: #include <petscsnes.h>

 52: /*
 53:    User-defined routines
 54: */
 55: extern PetscErrorCode FormMatrix(DM, Mat);
 56: extern PetscErrorCode MyComputeFunction(SNES, Vec, Vec, void *);
 57: extern PetscErrorCode MyComputeJacobian(SNES, Vec, Mat, Mat, void *);
 58: extern PetscErrorCode NonlinearGS(SNES, Vec);

 60: int main(int argc, char **argv)
 61: {
 62:   SNES      snes;  /* nonlinear solver */
 63:   SNES      psnes; /* nonlinear Gauss-Seidel approximate solver */
 64:   Vec       x, b;  /* solution vector */
 65:   PetscInt  its;   /* iterations for convergence */
 66:   DM        da;
 67:   PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */

 69:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 70:      Initialize program
 71:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 73:   PetscFunctionBeginUser;
 74:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:      Create nonlinear solver context
 78:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 79:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

 81:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_ngs_as_npc", &use_ngs_as_npc, 0));

 83:   if (use_ngs_as_npc) {
 84:     PetscCall(SNESGetNPC(snes, &psnes));
 85:     PetscCall(SNESSetType(psnes, SNESSHELL));
 86:     PetscCall(SNESShellSetSolve(psnes, NonlinearGS));
 87:   }

 89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90:      Create distributed array (DMDA) to manage parallel grid and vectors
 91:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 92:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
 93:   PetscCall(DMSetFromOptions(da));
 94:   PetscCall(DMSetUp(da));
 95:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
 96:   PetscCall(SNESSetDM(snes, da));
 97:   if (use_ngs_as_npc) PetscCall(SNESShellSetContext(psnes, da));
 98:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:      Extract global vectors from DMDA; then duplicate for remaining
100:      vectors that are the same types
101:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102:   PetscCall(DMCreateGlobalVector(da, &x));
103:   PetscCall(DMCreateGlobalVector(da, &b));
104:   PetscCall(VecSet(b, 1.0));

106:   PetscCall(SNESSetFunction(snes, NULL, MyComputeFunction, NULL));
107:   PetscCall(SNESSetJacobian(snes, NULL, NULL, MyComputeJacobian, NULL));

109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:      Customize nonlinear solver; set runtime options
111:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112:   PetscCall(SNESSetFromOptions(snes));

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Solve nonlinear system
116:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   PetscCall(SNESSolve(snes, b, x));
118:   PetscCall(SNESGetIterationNumber(snes, &its));

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

123:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124:      Free work space.  All PETSc objects should be destroyed when they
125:      are no longer needed.
126:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127:   PetscCall(VecDestroy(&x));
128:   PetscCall(VecDestroy(&b));
129:   PetscCall(SNESDestroy(&snes));
130:   PetscCall(DMDestroy(&da));
131:   PetscCall(PetscFinalize());
132:   return 0;
133: }

135: /* ------------------------------------------------------------------- */
136: PetscErrorCode MyComputeFunction(SNES snes, Vec x, Vec F, void *ctx)
137: {
138:   Mat J;
139:   DM  dm;

141:   PetscFunctionBeginUser;
142:   PetscCall(SNESGetDM(snes, &dm));
143:   PetscCall(PetscObjectQuery((PetscObject)dm, "_ex35_J", (PetscObject *)&J));
144:   if (!J) {
145:     PetscCall(DMSetMatType(dm, MATAIJ));
146:     PetscCall(DMCreateMatrix(dm, &J));
147:     PetscCall(MatSetDM(J, NULL));
148:     PetscCall(FormMatrix(dm, J));
149:     PetscCall(PetscObjectCompose((PetscObject)dm, "_ex35_J", (PetscObject)J));
150:     PetscCall(PetscObjectDereference((PetscObject)J));
151:   }
152:   PetscCall(MatMult(J, x, F));
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

156: PetscErrorCode MyComputeJacobian(SNES snes, Vec x, Mat J, Mat Jp, void *ctx)
157: {
158:   DM dm;

160:   PetscFunctionBeginUser;
161:   PetscCall(SNESGetDM(snes, &dm));
162:   PetscCall(FormMatrix(dm, Jp));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: PetscErrorCode FormMatrix(DM da, Mat jac)
167: {
168:   PetscInt      i, j, nrows = 0;
169:   MatStencil    col[5], row, *rows;
170:   PetscScalar   v[5], hx, hy, hxdhy, hydhx;
171:   DMDALocalInfo info;

173:   PetscFunctionBeginUser;
174:   PetscCall(DMDAGetLocalInfo(da, &info));
175:   hx    = 1.0 / (PetscReal)(info.mx - 1);
176:   hy    = 1.0 / (PetscReal)(info.my - 1);
177:   hxdhy = hx / hy;
178:   hydhx = hy / hx;

180:   PetscCall(PetscMalloc1(info.ym * info.xm, &rows));
181:   /*
182:      Compute entries for the locally owned part of the Jacobian.
183:       - Currently, all PETSc parallel matrix formats are partitioned by
184:         contiguous chunks of rows across the processors.
185:       - Each processor needs to insert only elements that it owns
186:         locally (but any non-local elements will be sent to the
187:         appropriate processor during matrix assembly).
188:       - Here, we set all entries for a particular row at once.
189:       - We can set matrix entries either using either
190:         MatSetValuesLocal() or MatSetValues(), as discussed above.
191:   */
192:   for (j = info.ys; j < info.ys + info.ym; j++) {
193:     for (i = info.xs; i < info.xs + info.xm; i++) {
194:       row.j = j;
195:       row.i = i;
196:       /* boundary points */
197:       if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) {
198:         v[0] = 2.0 * (hydhx + hxdhy);
199:         PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
200:         rows[nrows].i   = i;
201:         rows[nrows++].j = j;
202:       } else {
203:         /* interior grid points */
204:         v[0]     = -hxdhy;
205:         col[0].j = j - 1;
206:         col[0].i = i;
207:         v[1]     = -hydhx;
208:         col[1].j = j;
209:         col[1].i = i - 1;
210:         v[2]     = 2.0 * (hydhx + hxdhy);
211:         col[2].j = row.j;
212:         col[2].i = row.i;
213:         v[3]     = -hydhx;
214:         col[3].j = j;
215:         col[3].i = i + 1;
216:         v[4]     = -hxdhy;
217:         col[4].j = j + 1;
218:         col[4].i = i;
219:         PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES));
220:       }
221:     }
222:   }

224:   /*
225:      Assemble matrix, using the 2-step process:
226:        MatAssemblyBegin(), MatAssemblyEnd().
227:   */
228:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
229:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
230:   PetscCall(MatZeroRowsColumnsStencil(jac, nrows, rows, 2.0 * (hydhx + hxdhy), NULL, NULL));
231:   PetscCall(PetscFree(rows));
232:   /*
233:      Tell the matrix we will never add a new nonzero location to the
234:      matrix. If we do, it will generate an error.
235:   */
236:   PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: /* ------------------------------------------------------------------- */
241: /*
242:       Applies some sweeps on nonlinear Gauss-Seidel on each process

244:  */
245: PetscErrorCode NonlinearGS(SNES snes, Vec X)
246: {
247:   PetscInt      i, j, Mx, My, xs, ys, xm, ym, its, l;
248:   PetscReal     hx, hy, hxdhy, hydhx;
249:   PetscScalar **x, F, J, u, uxx, uyy;
250:   DM            da;
251:   Vec           localX;

253:   PetscFunctionBeginUser;
254:   PetscCall(SNESGetTolerances(snes, NULL, NULL, NULL, &its, NULL));
255:   PetscCall(SNESShellGetContext(snes, &da));

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

259:   hx    = 1.0 / (PetscReal)(Mx - 1);
260:   hy    = 1.0 / (PetscReal)(My - 1);
261:   hxdhy = hx / hy;
262:   hydhx = hy / hx;

264:   PetscCall(DMGetLocalVector(da, &localX));

266:   for (l = 0; l < its; l++) {
267:     PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
268:     PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
269:     /*
270:      Get a pointer to vector data.
271:      - For default PETSc vectors, VecGetArray() returns a pointer to
272:      the data array.  Otherwise, the routine is implementation dependent.
273:      - You MUST call VecRestoreArray() when you no longer need access to
274:      the array.
275:      */
276:     PetscCall(DMDAVecGetArray(da, localX, &x));

278:     /*
279:      Get local grid boundaries (for 2-dimensional DMDA):
280:      xs, ys   - starting grid indices (no ghost points)
281:      xm, ym   - widths of local grid (no ghost points)

283:      */
284:     PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

286:     for (j = ys; j < ys + ym; j++) {
287:       for (i = xs; i < xs + xm; i++) {
288:         if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
289:           /* boundary conditions are all zero Dirichlet */
290:           x[j][i] = 0.0;
291:         } else {
292:           u   = x[j][i];
293:           uxx = (2.0 * u - x[j][i - 1] - x[j][i + 1]) * hydhx;
294:           uyy = (2.0 * u - x[j - 1][i] - x[j + 1][i]) * hxdhy;
295:           F   = uxx + uyy;
296:           J   = 2.0 * (hydhx + hxdhy);
297:           u   = u - F / J;

299:           x[j][i] = u;
300:         }
301:       }
302:     }

304:     /*
305:      Restore vector
306:      */
307:     PetscCall(DMDAVecRestoreArray(da, localX, &x));
308:     PetscCall(DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X));
309:     PetscCall(DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X));
310:   }
311:   PetscCall(DMRestoreLocalVector(da, &localX));
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: /*TEST

317:    test:
318:       args: -snes_monitor_short -snes_type nrichardson
319:       requires: !single

321:    test:
322:       suffix: 2
323:       args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale
324:       requires: !single

326:    test:
327:       suffix: 3
328:       args: -snes_monitor_short -snes_type ngmres

330:    test:
331:       suffix: 4
332:       args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none

334:    test:
335:       suffix: 5
336:       args: -snes_monitor_short -snes_type ncg

338:    test:
339:       suffix: 6
340:       args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none

342:    test:
343:       suffix: 7
344:       args: -da_refine 2 -snes_monitor_short -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor_short -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor_short
345:       requires: !single

347:    test:
348:       suffix: 8
349:       args: -da_refine 2 -snes_monitor_short -snes_type fas -fas_levels_snes_monitor_short -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_type fas -snes_rtol 1.e-5

351:    test:
352:       suffix: 9
353:       args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc

355:    test:
356:       suffix: 10
357:       args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc -snes_trdc_use_cauchy false

359: TEST*/