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: *obj = lobj;
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, NULL, 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, (DMDASNESFunctionVecFn *)FormFunctionLocalVec, &user));
656: } else {
657: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, &user));
658: }
660: flg = PETSC_FALSE;
661: PetscCall(PetscOptionsGetBool(NULL, NULL, "-fd", &flg, NULL));
662: if (!flg) {
663: if (useKokkos) PetscCall(DMDASNESSetJacobianLocalVec(da, (DMDASNESJacobianVecFn *)FormJacobianLocalVec, &user));
664: else PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, &user));
665: }
667: flg = PETSC_FALSE;
668: PetscCall(PetscOptionsGetBool(NULL, NULL, "-obj", &flg, NULL));
669: if (flg) {
670: if (useKokkos) PetscCall(DMDASNESSetObjectiveLocalVec(da, (DMDASNESObjectiveVecFn *)FormObjectiveLocalVec, &user));
671: else PetscCall(DMDASNESSetObjectiveLocal(da, (DMDASNESObjectiveFn *)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*/