Actual source code: ex55.c

  1: static char help[] = "Copy of ex5.c\n";

  3: /* ------------------------------------------------------------------------

  5:   Copy of ex5.c.
  6:   Once petsc test harness supports conditional linking, we can remove this duplicate.
  7:   See https://gitlab.com/petsc/petsc/-/issues/1173
  8:   ------------------------------------------------------------------------- */

 10: /*
 11:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 12:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 13: */
 14: #include <petscdm.h>
 15: #include <petscdmda.h>
 16: #include <petscsnes.h>
 17: #include <petscmatlab.h>
 18: #include <petsc/private/snesimpl.h>
 19: #include "ex55.h"

 21: /* ------------------------------------------------------------------- */
 22: /*
 23:    FormInitialGuess - Forms initial approximation.

 25:    Input Parameters:
 26:    da - The DM
 27:    user - user-defined application context

 29:    Output Parameter:
 30:    X - vector
 31:  */
 32: static PetscErrorCode FormInitialGuess(DM da, AppCtx *user, Vec X)
 33: {
 34:   PetscInt      i, j, Mx, My, xs, ys, xm, ym;
 35:   PetscReal     lambda, temp1, temp, hx, hy;
 36:   PetscScalar **x;

 38:   PetscFunctionBeginUser;
 39:   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));

 41:   lambda = user->param;
 42:   hx     = 1.0 / (PetscReal)(Mx - 1);
 43:   hy     = 1.0 / (PetscReal)(My - 1);
 44:   temp1  = lambda / (lambda + 1.0);

 46:   /*
 47:      Get a pointer to vector data.
 48:        - For default PETSc vectors, VecGetArray() returns a pointer to
 49:          the data array.  Otherwise, the routine is implementation dependent.
 50:        - You MUST call VecRestoreArray() when you no longer need access to
 51:          the array.
 52:   */
 53:   PetscCall(DMDAVecGetArray(da, X, &x));

 55:   /*
 56:      Get local grid boundaries (for 2-dimensional DMDA):
 57:        xs, ys   - starting grid indices (no ghost points)
 58:        xm, ym   - widths of local grid (no ghost points)

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

 63:   /*
 64:      Compute initial guess over the locally owned part of the grid
 65:   */
 66:   for (j = ys; j < ys + ym; j++) {
 67:     temp = (PetscReal)(PetscMin(j, My - j - 1)) * hy;
 68:     for (i = xs; i < xs + xm; i++) {
 69:       if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
 70:         /* boundary conditions are all zero Dirichlet */
 71:         x[j][i] = 0.0;
 72:       } else {
 73:         x[j][i] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, Mx - i - 1)) * hx, temp));
 74:       }
 75:     }
 76:   }

 78:   /*
 79:      Restore vector
 80:   */
 81:   PetscCall(DMDAVecRestoreArray(da, X, &x));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: /*
 86:   FormExactSolution - Forms MMS solution

 88:   Input Parameters:
 89:   da - The DM
 90:   user - user-defined application context

 92:   Output Parameter:
 93:   X - vector
 94:  */
 95: static PetscErrorCode FormExactSolution(DM da, AppCtx *user, Vec U)
 96: {
 97:   DM            coordDA;
 98:   Vec           coordinates;
 99:   DMDACoor2d  **coords;
100:   PetscScalar **u;
101:   PetscInt      xs, ys, xm, ym, i, j;

103:   PetscFunctionBeginUser;
104:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
105:   PetscCall(DMGetCoordinateDM(da, &coordDA));
106:   PetscCall(DMGetCoordinates(da, &coordinates));
107:   PetscCall(DMDAVecGetArray(coordDA, coordinates, &coords));
108:   PetscCall(DMDAVecGetArray(da, U, &u));
109:   for (j = ys; j < ys + ym; ++j) {
110:     for (i = xs; i < xs + xm; ++i) PetscCall(user->mms_solution(user, &coords[j][i], &u[j][i]));
111:   }
112:   PetscCall(DMDAVecRestoreArray(da, U, &u));
113:   PetscCall(DMDAVecRestoreArray(coordDA, coordinates, &coords));
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static PetscErrorCode ZeroBCSolution(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
118: {
119:   u[0] = 0.;
120:   return PETSC_SUCCESS;
121: }

123: /* The functions below evaluate the MMS solution u(x,y) and associated forcing

125:      f(x,y) = -u_xx - u_yy - lambda exp(u)

127:   such that u(x,y) is an exact solution with f(x,y) as the right hand side forcing term.
128:  */
129: static PetscErrorCode MMSSolution1(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
130: {
131:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
132:   u[0] = x * (1 - x) * y * (1 - y);
133:   PetscCall(PetscLogFlops(5));
134:   return PETSC_SUCCESS;
135: }
136: static PetscErrorCode MMSForcing1(AppCtx *user, const DMDACoor2d *c, PetscScalar *f)
137: {
138:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
139:   f[0] = 2 * x * (1 - x) + 2 * y * (1 - y) - user->param * PetscExpReal(x * (1 - x) * y * (1 - y));
140:   return PETSC_SUCCESS;
141: }

143: static PetscErrorCode MMSSolution2(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
144: {
145:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
146:   u[0] = PetscSinReal(PETSC_PI * x) * PetscSinReal(PETSC_PI * y);
147:   PetscCall(PetscLogFlops(5));
148:   return PETSC_SUCCESS;
149: }
150: static PetscErrorCode MMSForcing2(AppCtx *user, const DMDACoor2d *c, PetscScalar *f)
151: {
152:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
153:   f[0] = 2 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x) * PetscSinReal(PETSC_PI * y) - user->param * PetscExpReal(PetscSinReal(PETSC_PI * x) * PetscSinReal(PETSC_PI * y));
154:   return PETSC_SUCCESS;
155: }

157: static PetscErrorCode MMSSolution3(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
158: {
159:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
160:   u[0] = PetscSinReal(user->m * PETSC_PI * x * (1 - y)) * PetscSinReal(user->n * PETSC_PI * y * (1 - x));
161:   PetscCall(PetscLogFlops(5));
162:   return PETSC_SUCCESS;
163: }
164: static PetscErrorCode MMSForcing3(AppCtx *user, const DMDACoor2d *c, PetscScalar *f)
165: {
166:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
167:   PetscReal m = user->m, n = user->n, lambda = user->param;
168:   f[0] = (-(PetscExpReal(PetscSinReal(m * PETSC_PI * x * (1 - y)) * PetscSinReal(n * PETSC_PI * (1 - x) * y)) * lambda) + PetscSqr(PETSC_PI) * (-2 * m * n * ((-1 + x) * x + (-1 + y) * y) * PetscCosReal(m * PETSC_PI * x * (-1 + y)) * PetscCosReal(n * PETSC_PI * (-1 + x) * y) + (PetscSqr(m) * (PetscSqr(x) + PetscSqr(-1 + y)) + PetscSqr(n) * (PetscSqr(-1 + x) + PetscSqr(y))) * PetscSinReal(m * PETSC_PI * x * (-1 + y)) * PetscSinReal(n * PETSC_PI * (-1 + x) * y)));
169:   return PETSC_SUCCESS;
170: }

172: static PetscErrorCode MMSSolution4(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
173: {
174:   const PetscReal Lx = 1., Ly = 1.;
175:   PetscReal       x = PetscRealPart(c->x), y = PetscRealPart(c->y);
176:   u[0] = (PetscPowReal(x, 4) - PetscSqr(Lx) * PetscSqr(x)) * (PetscPowReal(y, 4) - PetscSqr(Ly) * PetscSqr(y));
177:   PetscCall(PetscLogFlops(9));
178:   return PETSC_SUCCESS;
179: }
180: static PetscErrorCode MMSForcing4(AppCtx *user, const DMDACoor2d *c, PetscScalar *f)
181: {
182:   const PetscReal Lx = 1., Ly = 1.;
183:   PetscReal       x = PetscRealPart(c->x), y = PetscRealPart(c->y);
184:   f[0] = (2 * PetscSqr(x) * (PetscSqr(x) - PetscSqr(Lx)) * (PetscSqr(Ly) - 6 * PetscSqr(y)) + 2 * PetscSqr(y) * (PetscSqr(Lx) - 6 * PetscSqr(x)) * (PetscSqr(y) - PetscSqr(Ly)) - user->param * PetscExpReal((PetscPowReal(x, 4) - PetscSqr(Lx) * PetscSqr(x)) * (PetscPowReal(y, 4) - PetscSqr(Ly) * PetscSqr(y))));
185:   return PETSC_SUCCESS;
186: }

188: /* ------------------------------------------------------------------- */
189: /*
190:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch

192:  */
193: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
194: {
195:   PetscInt    i, j;
196:   PetscReal   lambda, hx, hy, hxdhy, hydhx;
197:   PetscScalar u, ue, uw, un, us, uxx, uyy, mms_solution, mms_forcing;
198:   DMDACoor2d  c;

200:   PetscFunctionBeginUser;
201:   lambda = user->param;
202:   hx     = 1.0 / (PetscReal)(info->mx - 1);
203:   hy     = 1.0 / (PetscReal)(info->my - 1);
204:   hxdhy  = hx / hy;
205:   hydhx  = hy / hx;
206:   /*
207:      Compute function over the locally owned part of the grid
208:   */
209:   for (j = info->ys; j < info->ys + info->ym; j++) {
210:     for (i = info->xs; i < info->xs + info->xm; i++) {
211:       if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
212:         c.x = i * hx;
213:         c.y = j * hy;
214:         PetscCall(user->mms_solution(user, &c, &mms_solution));
215:         f[j][i] = 2.0 * (hydhx + hxdhy) * (x[j][i] - mms_solution);
216:       } else {
217:         u  = x[j][i];
218:         uw = x[j][i - 1];
219:         ue = x[j][i + 1];
220:         un = x[j - 1][i];
221:         us = x[j + 1][i];

223:         /* Enforce boundary conditions at neighboring points -- setting these values causes the Jacobian to be symmetric. */
224:         if (i - 1 == 0) {
225:           c.x = (i - 1) * hx;
226:           c.y = j * hy;
227:           PetscCall(user->mms_solution(user, &c, &uw));
228:         }
229:         if (i + 1 == info->mx - 1) {
230:           c.x = (i + 1) * hx;
231:           c.y = j * hy;
232:           PetscCall(user->mms_solution(user, &c, &ue));
233:         }
234:         if (j - 1 == 0) {
235:           c.x = i * hx;
236:           c.y = (j - 1) * hy;
237:           PetscCall(user->mms_solution(user, &c, &un));
238:         }
239:         if (j + 1 == info->my - 1) {
240:           c.x = i * hx;
241:           c.y = (j + 1) * hy;
242:           PetscCall(user->mms_solution(user, &c, &us));
243:         }

245:         uxx         = (2.0 * u - uw - ue) * hydhx;
246:         uyy         = (2.0 * u - un - us) * hxdhy;
247:         mms_forcing = 0;
248:         c.x         = i * hx;
249:         c.y         = j * hy;
250:         if (user->mms_forcing) PetscCall(user->mms_forcing(user, &c, &mms_forcing));
251:         f[j][i] = uxx + uyy - hx * hy * (lambda * PetscExpScalar(u) + mms_forcing);
252:       }
253:     }
254:   }
255:   PetscCall(PetscLogFlops(11.0 * info->ym * info->xm));
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

259: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
260: static PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info, PetscScalar **x, PetscReal *obj, AppCtx *user)
261: {
262:   PetscInt    i, j;
263:   PetscReal   lambda, hx, hy, hxdhy, hydhx, sc, lobj = 0;
264:   PetscScalar u, ue, uw, un, us, uxux, uyuy;
265:   MPI_Comm    comm;

267:   PetscFunctionBeginUser;
268:   *obj = 0;
269:   PetscCall(PetscObjectGetComm((PetscObject)info->da, &comm));
270:   lambda = user->param;
271:   hx     = 1.0 / (PetscReal)(info->mx - 1);
272:   hy     = 1.0 / (PetscReal)(info->my - 1);
273:   sc     = hx * hy * lambda;
274:   hxdhy  = hx / hy;
275:   hydhx  = hy / hx;
276:   /*
277:      Compute function over the locally owned part of the grid
278:   */
279:   for (j = info->ys; j < info->ys + info->ym; j++) {
280:     for (i = info->xs; i < info->xs + info->xm; i++) {
281:       if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
282:         lobj += PetscRealPart((hydhx + hxdhy) * x[j][i] * x[j][i]);
283:       } else {
284:         u  = x[j][i];
285:         uw = x[j][i - 1];
286:         ue = x[j][i + 1];
287:         un = x[j - 1][i];
288:         us = x[j + 1][i];

290:         if (i - 1 == 0) uw = 0.;
291:         if (i + 1 == info->mx - 1) ue = 0.;
292:         if (j - 1 == 0) un = 0.;
293:         if (j + 1 == info->my - 1) us = 0.;

295:         /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */

297:         uxux = u * (2. * u - ue - uw) * hydhx;
298:         uyuy = u * (2. * u - un - us) * hxdhy;

300:         lobj += PetscRealPart(0.5 * (uxux + uyuy) - sc * PetscExpScalar(u));
301:       }
302:     }
303:   }
304:   PetscCall(PetscLogFlops(12.0 * info->ym * info->xm));
305:   PetscCall(MPIU_Allreduce(&lobj, obj, 1, MPIU_REAL, MPIU_SUM, comm));
306:   PetscFunctionReturn(PETSC_SUCCESS);
307: }

309: /*
310:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
311: */
312: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **x, Mat jac, Mat jacpre, AppCtx *user)
313: {
314:   PetscInt     i, j, k;
315:   MatStencil   col[5], row;
316:   PetscScalar  lambda, v[5], hx, hy, hxdhy, hydhx, sc;
317:   DM           coordDA;
318:   Vec          coordinates;
319:   DMDACoor2d **coords;

321:   PetscFunctionBeginUser;
322:   lambda = user->param;
323:   /* Extract coordinates */
324:   PetscCall(DMGetCoordinateDM(info->da, &coordDA));
325:   PetscCall(DMGetCoordinates(info->da, &coordinates));
326:   PetscCall(DMDAVecGetArray(coordDA, coordinates, &coords));
327:   hx = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs + 1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
328:   hy = info->ym > 1 ? PetscRealPart(coords[info->ys + 1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
329:   PetscCall(DMDAVecRestoreArray(coordDA, coordinates, &coords));
330:   hxdhy = hx / hy;
331:   hydhx = hy / hx;
332:   sc    = hx * hy * lambda;

334:   /*
335:      Compute entries for the locally owned part of the Jacobian.
336:       - Currently, all PETSc parallel matrix formats are partitioned by
337:         contiguous chunks of rows across the processors.
338:       - Each processor needs to insert only elements that it owns
339:         locally (but any non-local elements will be sent to the
340:         appropriate processor during matrix assembly).
341:       - Here, we set all entries for a particular row at once.
342:       - We can set matrix entries either using either
343:         MatSetValuesLocal() or MatSetValues(), as discussed above.
344:   */
345:   for (j = info->ys; j < info->ys + info->ym; j++) {
346:     for (i = info->xs; i < info->xs + info->xm; i++) {
347:       row.j = j;
348:       row.i = i;
349:       /* boundary points */
350:       if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
351:         v[0] = 2.0 * (hydhx + hxdhy);
352:         PetscCall(MatSetValuesStencil(jacpre, 1, &row, 1, &row, v, INSERT_VALUES));
353:       } else {
354:         k = 0;
355:         /* interior grid points */
356:         if (j - 1 != 0) {
357:           v[k]     = -hxdhy;
358:           col[k].j = j - 1;
359:           col[k].i = i;
360:           k++;
361:         }
362:         if (i - 1 != 0) {
363:           v[k]     = -hydhx;
364:           col[k].j = j;
365:           col[k].i = i - 1;
366:           k++;
367:         }

369:         v[k]     = 2.0 * (hydhx + hxdhy) - sc * PetscExpScalar(x[j][i]);
370:         col[k].j = row.j;
371:         col[k].i = row.i;
372:         k++;

374:         if (i + 1 != info->mx - 1) {
375:           v[k]     = -hydhx;
376:           col[k].j = j;
377:           col[k].i = i + 1;
378:           k++;
379:         }
380:         if (j + 1 != info->mx - 1) {
381:           v[k]     = -hxdhy;
382:           col[k].j = j + 1;
383:           col[k].i = i;
384:           k++;
385:         }
386:         PetscCall(MatSetValuesStencil(jacpre, 1, &row, k, col, v, INSERT_VALUES));
387:       }
388:     }
389:   }

391:   /*
392:      Assemble matrix, using the 2-step process:
393:        MatAssemblyBegin(), MatAssemblyEnd().
394:   */
395:   PetscCall(MatAssemblyBegin(jacpre, MAT_FINAL_ASSEMBLY));
396:   PetscCall(MatAssemblyEnd(jacpre, MAT_FINAL_ASSEMBLY));

398:   /*
399:      Tell the matrix we will never add a new nonzero location to the
400:      matrix. If we do, it will generate an error.
401:   */
402:   PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
403:   PetscFunctionReturn(PETSC_SUCCESS);
404: }

406: static PetscErrorCode FormFunctionMatlab(SNES snes, Vec X, Vec F, void *ptr)
407: {
408: #if PetscDefined(HAVE_MATLAB)
409:   AppCtx   *user = (AppCtx *)ptr;
410:   PetscInt  Mx, My;
411:   PetscReal lambda, hx, hy;
412:   Vec       localX, localF;
413:   MPI_Comm  comm;
414:   DM        da;

416:   PetscFunctionBeginUser;
417:   PetscCall(SNESGetDM(snes, &da));
418:   PetscCall(DMGetLocalVector(da, &localX));
419:   PetscCall(DMGetLocalVector(da, &localF));
420:   PetscCall(PetscObjectSetName((PetscObject)localX, "localX"));
421:   PetscCall(PetscObjectSetName((PetscObject)localF, "localF"));
422:   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));

424:   lambda = user->param;
425:   hx     = 1.0 / (PetscReal)(Mx - 1);
426:   hy     = 1.0 / (PetscReal)(My - 1);

428:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
429:   /*
430:      Scatter ghost points to local vector,using the 2-step process
431:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
432:      By placing code between these two statements, computations can be
433:      done while messages are in transition.
434:   */
435:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
436:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
437:   PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm), (PetscObject)localX));
438:   PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm), "localF=ex5m(localX,%18.16e,%18.16e,%18.16e)", (double)hx, (double)hy, (double)lambda));
439:   PetscCall(PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm), (PetscObject)localF));

441:   /*
442:      Insert values into global vector
443:   */
444:   PetscCall(DMLocalToGlobalBegin(da, localF, INSERT_VALUES, F));
445:   PetscCall(DMLocalToGlobalEnd(da, localF, INSERT_VALUES, F));
446:   PetscCall(DMRestoreLocalVector(da, &localX));
447:   PetscCall(DMRestoreLocalVector(da, &localF));
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: #else
450:   return PETSC_SUCCESS; /* Never called */
451: #endif
452: }

454: /* ------------------------------------------------------------------- */
455: /*
456:       Applies some sweeps on nonlinear Gauss-Seidel on each process

458:  */
459: static PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ctx)
460: {
461:   PetscInt      i, j, k, Mx, My, xs, ys, xm, ym, its, tot_its, sweeps, l;
462:   PetscReal     lambda, hx, hy, hxdhy, hydhx, sc;
463:   PetscScalar **x, **b, bij, F, F0 = 0, J, u, un, us, ue, eu, uw, uxx, uyy, y;
464:   PetscReal     atol, rtol, stol;
465:   DM            da;
466:   AppCtx       *user;
467:   Vec           localX, localB;

469:   PetscFunctionBeginUser;
470:   tot_its = 0;
471:   PetscCall(SNESNGSGetSweeps(snes, &sweeps));
472:   PetscCall(SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its));
473:   PetscCall(SNESGetDM(snes, &da));
474:   PetscCall(DMGetApplicationContext(da, &user));

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

478:   lambda = user->param;
479:   hx     = 1.0 / (PetscReal)(Mx - 1);
480:   hy     = 1.0 / (PetscReal)(My - 1);
481:   sc     = hx * hy * lambda;
482:   hxdhy  = hx / hy;
483:   hydhx  = hy / hx;

485:   PetscCall(DMGetLocalVector(da, &localX));
486:   if (B) PetscCall(DMGetLocalVector(da, &localB));
487:   for (l = 0; l < sweeps; l++) {
488:     PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
489:     PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
490:     if (B) {
491:       PetscCall(DMGlobalToLocalBegin(da, B, INSERT_VALUES, localB));
492:       PetscCall(DMGlobalToLocalEnd(da, B, INSERT_VALUES, localB));
493:     }
494:     /*
495:      Get a pointer to vector data.
496:      - For default PETSc vectors, VecGetArray() returns a pointer to
497:      the data array.  Otherwise, the routine is implementation dependent.
498:      - You MUST call VecRestoreArray() when you no longer need access to
499:      the array.
500:      */
501:     PetscCall(DMDAVecGetArray(da, localX, &x));
502:     if (B) PetscCall(DMDAVecGetArray(da, localB, &b));
503:     /*
504:      Get local grid boundaries (for 2-dimensional DMDA):
505:      xs, ys   - starting grid indices (no ghost points)
506:      xm, ym   - widths of local grid (no ghost points)
507:      */
508:     PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

510:     for (j = ys; j < ys + ym; j++) {
511:       for (i = xs; i < xs + xm; i++) {
512:         if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
513:           /* boundary conditions are all zero Dirichlet */
514:           x[j][i] = 0.0;
515:         } else {
516:           if (B) bij = b[j][i];
517:           else bij = 0.;

519:           u  = x[j][i];
520:           un = x[j - 1][i];
521:           us = x[j + 1][i];
522:           ue = x[j][i - 1];
523:           uw = x[j][i + 1];

525:           for (k = 0; k < its; k++) {
526:             eu  = PetscExpScalar(u);
527:             uxx = (2.0 * u - ue - uw) * hydhx;
528:             uyy = (2.0 * u - un - us) * hxdhy;
529:             F   = uxx + uyy - sc * eu - bij;
530:             if (k == 0) F0 = F;
531:             J = 2.0 * (hydhx + hxdhy) - sc * eu;
532:             y = F / J;
533:             u -= y;
534:             tot_its++;

536:             if (atol > PetscAbsReal(PetscRealPart(F)) || rtol * PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) || stol * PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) break;
537:           }
538:           x[j][i] = u;
539:         }
540:       }
541:     }
542:     /*
543:      Restore vector
544:      */
545:     PetscCall(DMDAVecRestoreArray(da, localX, &x));
546:     PetscCall(DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X));
547:     PetscCall(DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X));
548:   }
549:   PetscCall(PetscLogFlops(tot_its * (21.0)));
550:   PetscCall(DMRestoreLocalVector(da, &localX));
551:   if (B) {
552:     PetscCall(DMDAVecRestoreArray(da, localB, &b));
553:     PetscCall(DMRestoreLocalVector(da, &localB));
554:   }
555:   PetscFunctionReturn(PETSC_SUCCESS);
556: }

558: int main(int argc, char **argv)
559: {
560:   SNES      snes; /* nonlinear solver */
561:   Vec       x;    /* solution vector */
562:   AppCtx    user; /* user-defined work context */
563:   PetscInt  its;  /* iterations for convergence */
564:   PetscReal bratu_lambda_max = 6.81;
565:   PetscReal bratu_lambda_min = 0.;
566:   PetscInt  MMS              = 1;
567:   PetscBool flg, setMMS;
568:   DM        da;
569:   Vec       r = NULL;
570:   KSP       ksp;
571:   PetscInt  lits, slits;
572:   PetscBool useKokkos = PETSC_FALSE;

574:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
575:      Initialize program
576:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

578:   PetscFunctionBeginUser;
579:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

581:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_kokkos", &useKokkos, NULL));

583:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
584:      Initialize problem parameters
585:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
586:   user.ncoo  = 0;
587:   user.param = 6.0;
588:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL));
589:   PetscCheck(user.param <= bratu_lambda_max && user.param >= bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda, %g, is out of range, [%g, %g]", (double)user.param, (double)bratu_lambda_min, (double)bratu_lambda_max);
590:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mms", &MMS, &setMMS));
591:   if (MMS == 3) {
592:     PetscInt mPar = 2, nPar = 1;
593:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_par", &mPar, NULL));
594:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-n_par", &nPar, NULL));
595:     user.m = PetscPowInt(2, mPar);
596:     user.n = PetscPowInt(2, nPar);
597:   }

599:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
600:      Create nonlinear solver context
601:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
602:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
603:   PetscCall(SNESSetCountersReset(snes, PETSC_FALSE));
604:   PetscCall(SNESSetNGS(snes, NonlinearGS, NULL));

606:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
607:      Create distributed array (DMDA) to manage parallel grid and vectors
608:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
609:   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));
610:   if (useKokkos) {
611:     PetscCall(DMSetVecType(da, VECKOKKOS));
612:     PetscCall(DMSetMatType(da, MATAIJKOKKOS));
613:   }
614:   PetscCall(DMSetFromOptions(da));
615:   PetscCall(DMSetUp(da));
616:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
617:   PetscCall(DMSetApplicationContext(da, &user));
618:   PetscCall(SNESSetDM(snes, da));
619:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
620:      Extract global vectors from DMDA; then duplicate for remaining
621:      vectors that are the same types
622:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
623:   PetscCall(DMCreateGlobalVector(da, &x));

625:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
626:      Set local function evaluation routine
627:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
628:   switch (MMS) {
629:   case 0:
630:     user.mms_solution = ZeroBCSolution;
631:     user.mms_forcing  = NULL;
632:     break;
633:   case 1:
634:     user.mms_solution = MMSSolution1;
635:     user.mms_forcing  = MMSForcing1;
636:     break;
637:   case 2:
638:     user.mms_solution = MMSSolution2;
639:     user.mms_forcing  = MMSForcing2;
640:     break;
641:   case 3:
642:     user.mms_solution = MMSSolution3;
643:     user.mms_forcing  = MMSForcing3;
644:     break;
645:   case 4:
646:     user.mms_solution = MMSSolution4;
647:     user.mms_forcing  = MMSForcing4;
648:     break;
649:   default:
650:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Unknown MMS type %" PetscInt_FMT, MMS);
651:   }

653:   if (useKokkos) {
654:     PetscCheck(MMS == 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "FormFunctionLocalVec_Kokkos only works with MMS 1");
655:     PetscCall(DMDASNESSetFunctionLocalVec(da, INSERT_VALUES, (DMDASNESFunctionVec)FormFunctionLocalVec, &user));
656:   } else {
657:     PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunction)FormFunctionLocal, &user));
658:   }

660:   flg = PETSC_FALSE;
661:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-fd", &flg, NULL));
662:   if (!flg) {
663:     if (useKokkos) PetscCall(DMDASNESSetJacobianLocalVec(da, (DMDASNESJacobianVec)FormJacobianLocalVec, &user));
664:     else PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobian)FormJacobianLocal, &user));
665:   }

667:   flg = PETSC_FALSE;
668:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-obj", &flg, NULL));
669:   if (flg) {
670:     if (useKokkos) PetscCall(DMDASNESSetObjectiveLocalVec(da, (DMDASNESObjectiveVec)FormObjectiveLocalVec, &user));
671:     else PetscCall(DMDASNESSetObjectiveLocal(da, (DMDASNESObjective)FormObjectiveLocal, &user));
672:   }

674:   if (PetscDefined(HAVE_MATLAB)) {
675:     PetscBool matlab_function = PETSC_FALSE;
676:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-matlab_function", &matlab_function, 0));
677:     if (matlab_function) {
678:       PetscCall(VecDuplicate(x, &r));
679:       PetscCall(SNESSetFunction(snes, r, FormFunctionMatlab, &user));
680:     }
681:   }

683:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
684:      Customize nonlinear solver; set runtime options
685:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
686:   PetscCall(SNESSetFromOptions(snes));

688:   PetscCall(FormInitialGuess(da, &user, x));

690:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
691:      Solve nonlinear system
692:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
693:   PetscCall(SNESSolve(snes, NULL, x));
694:   PetscCall(SNESGetIterationNumber(snes, &its));

696:   PetscCall(SNESGetLinearSolveIterations(snes, &slits));
697:   PetscCall(SNESGetKSP(snes, &ksp));
698:   PetscCall(KSPGetTotalIterations(ksp, &lits));
699:   PetscCheck(lits == slits, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Number of total linear iterations reported by SNES %" PetscInt_FMT " does not match reported by KSP %" PetscInt_FMT, slits, lits);
700:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
701:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

703:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
704:      If using MMS, check the l_2 error
705:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
706:   if (setMMS) {
707:     Vec       e;
708:     PetscReal errorl2, errorinf;
709:     PetscInt  N;

711:     PetscCall(VecDuplicate(x, &e));
712:     PetscCall(PetscObjectViewFromOptions((PetscObject)x, NULL, "-sol_view"));
713:     PetscCall(FormExactSolution(da, &user, e));
714:     PetscCall(PetscObjectViewFromOptions((PetscObject)e, NULL, "-exact_view"));
715:     PetscCall(VecAXPY(e, -1.0, x));
716:     PetscCall(PetscObjectViewFromOptions((PetscObject)e, NULL, "-error_view"));
717:     PetscCall(VecNorm(e, NORM_2, &errorl2));
718:     PetscCall(VecNorm(e, NORM_INFINITY, &errorinf));
719:     PetscCall(VecGetSize(e, &N));
720:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " error L2 %g inf %g\n", N, (double)(errorl2 / PetscSqrtReal((PetscReal)N)), (double)errorinf));
721:     PetscCall(VecDestroy(&e));
722:     PetscCall(PetscLogEventSetDof(SNES_Solve, 0, N));
723:     PetscCall(PetscLogEventSetError(SNES_Solve, 0, errorl2 / PetscSqrtReal(N)));
724:   }

726:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
727:      Free work space.  All PETSc objects should be destroyed when they
728:      are no longer needed.
729:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
730:   PetscCall(VecDestroy(&r));
731:   PetscCall(VecDestroy(&x));
732:   PetscCall(SNESDestroy(&snes));
733:   PetscCall(DMDestroy(&da));
734:   PetscCall(PetscFinalize());
735:   return 0;
736: }

738: /*TEST
739:   build:
740:     requires: !windows_compilers
741:     depends: ex55k.kokkos.cxx

743:   testset:
744:     output_file: output/ex55_asm_0.out
745:     requires: !single
746:     args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu
747:     filter: grep -v "type"

749:     test:
750:       suffix: asm_0
751:       args: -fd {{0 1}}

753:     test:
754:       requires: kokkos_kernels
755:       suffix: asm_0_kok
756:       args: -use_kokkos -fd {{0 1}}

758:   testset:
759:     output_file: output/ex55_1.out
760:     requires: !single
761:     args: -snes_monitor
762:     filter: grep -v "type"

764:     test:
765:       suffix: 1
766:       args: -fd {{0 1}}

768:     test:
769:       requires: kokkos_kernels
770:       suffix: 1_kok
771:       args: -use_kokkos -fd {{0 1}}

773:     test:
774:       requires: h2opus
775:       suffix: 1_h2opus
776:       args: -pc_type h2opus -fd {{0 1}}

778:     test:
779:       requires: h2opus kokkos_kernels
780:       suffix: 1_h2opus_k
781:       args: -use_kokkos -pc_type h2opus -fd {{0 1}}

783: TEST*/