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;
39: 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: 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: 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: DMDAVecRestoreArray(da, X, &x);
82: return 0;
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;
104: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
105: DMGetCoordinateDM(da, &coordDA);
106: DMGetCoordinates(da, &coordinates);
107: DMDAVecGetArray(coordDA, coordinates, &coords);
108: DMDAVecGetArray(da, U, &u);
109: for (j = ys; j < ys + ym; ++j) {
110: for (i = xs; i < xs + xm; ++i) user->mms_solution(user, &coords[j][i], &u[j][i]);
111: }
112: DMDAVecRestoreArray(da, U, &u);
113: DMDAVecRestoreArray(coordDA, coordinates, &coords);
114: return 0;
115: }
117: static PetscErrorCode ZeroBCSolution(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
118: {
119: u[0] = 0.;
120: return 0;
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: PetscLogFlops(5);
134: return 0;
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 0;
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: PetscLogFlops(5);
148: return 0;
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 0;
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: PetscLogFlops(5);
162: return 0;
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 0;
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: PetscLogFlops(9);
178: return 0;
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 0;
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;
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: 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: 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: user->mms_solution(user, &c, &ue);
233: }
234: if (j - 1 == 0) {
235: c.x = i * hx;
236: c.y = (j - 1) * hy;
237: 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: 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) user->mms_forcing(user, &c, &mms_forcing);
251: f[j][i] = uxx + uyy - hx * hy * (lambda * PetscExpScalar(u) + mms_forcing);
252: }
253: }
254: }
255: PetscLogFlops(11.0 * info->ym * info->xm);
256: return 0;
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;
268: *obj = 0;
269: 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: PetscLogFlops(12.0 * info->ym * info->xm);
305: MPI_Allreduce(&lobj, obj, 1, MPIU_REAL, MPIU_SUM, comm);
306: return 0;
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;
322: lambda = user->param;
323: /* Extract coordinates */
324: DMGetCoordinateDM(info->da, &coordDA);
325: DMGetCoordinates(info->da, &coordinates);
326: 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: 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: 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: 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: MatAssemblyBegin(jacpre, MAT_FINAL_ASSEMBLY);
396: 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: MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
403: return 0;
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;
417: SNESGetDM(snes, &da);
418: DMGetLocalVector(da, &localX);
419: DMGetLocalVector(da, &localF);
420: PetscObjectSetName((PetscObject)localX, "localX");
421: PetscObjectSetName((PetscObject)localF, "localF");
422: 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: 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: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
436: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
437: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm), (PetscObject)localX);
438: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm), "localF=ex5m(localX,%18.16e,%18.16e,%18.16e)", (double)hx, (double)hy, (double)lambda);
439: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm), (PetscObject)localF);
441: /*
442: Insert values into global vector
443: */
444: DMLocalToGlobalBegin(da, localF, INSERT_VALUES, F);
445: DMLocalToGlobalEnd(da, localF, INSERT_VALUES, F);
446: DMRestoreLocalVector(da, &localX);
447: DMRestoreLocalVector(da, &localF);
448: return 0;
449: #else
450: return 0; /* 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;
470: tot_its = 0;
471: SNESNGSGetSweeps(snes, &sweeps);
472: SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its);
473: SNESGetDM(snes, &da);
474: DMGetApplicationContext(da, &user);
476: 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: DMGetLocalVector(da, &localX);
486: if (B) DMGetLocalVector(da, &localB);
487: for (l = 0; l < sweeps; l++) {
488: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
489: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
490: if (B) {
491: DMGlobalToLocalBegin(da, B, INSERT_VALUES, localB);
492: 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: DMDAVecGetArray(da, localX, &x);
502: if (B) 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: 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: DMDAVecRestoreArray(da, localX, &x);
546: DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X);
547: DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X);
548: }
549: PetscLogFlops(tot_its * (21.0));
550: DMRestoreLocalVector(da, &localX);
551: if (B) {
552: DMDAVecRestoreArray(da, localB, &b);
553: DMRestoreLocalVector(da, &localB);
554: }
555: return 0;
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 = PETSC_FALSE, 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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
579: PetscInitialize(&argc, &argv, (char *)0, help);
581: PetscOptionsGetBool(NULL, NULL, "-use_kokkos", &useKokkos, NULL);
583: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
584: Initialize problem parameters
585: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
586: user.param = 6.0;
587: PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL);
589: PetscOptionsGetInt(NULL, NULL, "-mms", &MMS, &setMMS);
590: if (MMS == 3) {
591: PetscInt mPar = 2, nPar = 1;
592: PetscOptionsGetInt(NULL, NULL, "-m_par", &mPar, NULL);
593: PetscOptionsGetInt(NULL, NULL, "-n_par", &nPar, NULL);
594: user.m = PetscPowInt(2, mPar);
595: user.n = PetscPowInt(2, nPar);
596: }
598: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
599: Create nonlinear solver context
600: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
601: SNESCreate(PETSC_COMM_WORLD, &snes);
602: SNESSetCountersReset(snes, PETSC_FALSE);
603: SNESSetNGS(snes, NonlinearGS, NULL);
605: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
606: Create distributed array (DMDA) to manage parallel grid and vectors
607: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
608: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);
609: DMSetFromOptions(da);
610: DMSetUp(da);
611: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
612: DMSetApplicationContext(da, &user);
613: SNESSetDM(snes, da);
614: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
615: Extract global vectors from DMDA; then duplicate for remaining
616: vectors that are the same types
617: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
618: DMCreateGlobalVector(da, &x);
620: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
621: Set local function evaluation routine
622: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
623: switch (MMS) {
624: case 0:
625: user.mms_solution = ZeroBCSolution;
626: user.mms_forcing = NULL;
627: break;
628: case 1:
629: user.mms_solution = MMSSolution1;
630: user.mms_forcing = MMSForcing1;
631: break;
632: case 2:
633: user.mms_solution = MMSSolution2;
634: user.mms_forcing = MMSForcing2;
635: break;
636: case 3:
637: user.mms_solution = MMSSolution3;
638: user.mms_forcing = MMSForcing3;
639: break;
640: case 4:
641: user.mms_solution = MMSSolution4;
642: user.mms_forcing = MMSForcing4;
643: break;
644: default:
645: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Unknown MMS type %" PetscInt_FMT, MMS);
646: }
648: if (useKokkos) {
650: DMDASNESSetFunctionLocalVec(da, INSERT_VALUES, (DMDASNESFunctionVec)FormFunctionLocalVec, &user);
651: } else {
652: DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunction)FormFunctionLocal, &user);
653: }
655: PetscOptionsGetBool(NULL, NULL, "-fd", &flg, NULL);
656: if (!flg) {
657: if (useKokkos) DMDASNESSetJacobianLocalVec(da, (DMDASNESJacobianVec)FormJacobianLocalVec, &user);
658: else DMDASNESSetJacobianLocal(da, (DMDASNESJacobian)FormJacobianLocal, &user);
659: }
661: PetscOptionsGetBool(NULL, NULL, "-obj", &flg, NULL);
662: if (flg) {
663: if (useKokkos) DMDASNESSetObjectiveLocalVec(da, (DMDASNESObjectiveVec)FormObjectiveLocalVec, &user);
664: else DMDASNESSetObjectiveLocal(da, (DMDASNESObjective)FormObjectiveLocal, &user);
665: }
667: if (PetscDefined(HAVE_MATLAB)) {
668: PetscBool matlab_function = PETSC_FALSE;
669: PetscOptionsGetBool(NULL, NULL, "-matlab_function", &matlab_function, 0);
670: if (matlab_function) {
671: VecDuplicate(x, &r);
672: SNESSetFunction(snes, r, FormFunctionMatlab, &user);
673: }
674: }
676: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
677: Customize nonlinear solver; set runtime options
678: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
679: SNESSetFromOptions(snes);
681: FormInitialGuess(da, &user, x);
683: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
684: Solve nonlinear system
685: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
686: SNESSolve(snes, NULL, x);
687: SNESGetIterationNumber(snes, &its);
689: SNESGetLinearSolveIterations(snes, &slits);
690: SNESGetKSP(snes, &ksp);
691: KSPGetTotalIterations(ksp, &lits);
693: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
694: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
696: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
697: If using MMS, check the l_2 error
698: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
699: if (setMMS) {
700: Vec e;
701: PetscReal errorl2, errorinf;
702: PetscInt N;
704: VecDuplicate(x, &e);
705: PetscObjectViewFromOptions((PetscObject)x, NULL, "-sol_view");
706: FormExactSolution(da, &user, e);
707: PetscObjectViewFromOptions((PetscObject)e, NULL, "-exact_view");
708: VecAXPY(e, -1.0, x);
709: PetscObjectViewFromOptions((PetscObject)e, NULL, "-error_view");
710: VecNorm(e, NORM_2, &errorl2);
711: VecNorm(e, NORM_INFINITY, &errorinf);
712: VecGetSize(e, &N);
713: PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " error L2 %g inf %g\n", N, (double)(errorl2 / PetscSqrtReal((PetscReal)N)), (double)errorinf);
714: VecDestroy(&e);
715: PetscLogEventSetDof(SNES_Solve, 0, N);
716: PetscLogEventSetError(SNES_Solve, 0, errorl2 / PetscSqrtReal(N));
717: }
719: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
720: Free work space. All PETSc objects should be destroyed when they
721: are no longer needed.
722: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
723: VecDestroy(&r);
724: VecDestroy(&x);
725: SNESDestroy(&snes);
726: DMDestroy(&da);
727: PetscFinalize();
728: return 0;
729: }
731: /*TEST
732: build:
733: requires: kokkos_kernels
734: depends: ex55k.kokkos.cxx
736: testset:
737: output_file: output/ex55_asm_0.out
738: requires: !single
739: 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
740: filter: grep -v "type"
742: test:
743: suffix: asm_0
744: test:
745: suffix: asm_0_kok
746: args: -use_kokkos 1 -dm_mat_type aijkokkos -dm_vec_type kokkos
748: TEST*/