1 static char help[] = "Copy of ex5.c\n";
2
3 /* ------------------------------------------------------------------------
4
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 ------------------------------------------------------------------------- */
9
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> /* For SNES_Solve event */
19 #include "ex55.h"
20
21 /* ------------------------------------------------------------------- */
22 /*
23 FormInitialGuess - Forms initial approximation.
24
25 Input Parameters:
26 da - The DM
27 user - user-defined application context
28
29 Output Parameter:
30 X - vector
31 */
FormInitialGuess(DM da,AppCtx * user,Vec X)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;
37
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));
40
41 lambda = user->param;
42 hx = 1.0 / (PetscReal)(Mx - 1);
43 hy = 1.0 / (PetscReal)(My - 1);
44 temp1 = lambda / (lambda + 1.0);
45
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));
54
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)
59
60 */
61 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
62
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 }
77
78 /*
79 Restore vector
80 */
81 PetscCall(DMDAVecRestoreArray(da, X, &x));
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84
85 /*
86 FormExactSolution - Forms MMS solution
87
88 Input Parameters:
89 da - The DM
90 user - user-defined application context
91
92 Output Parameter:
93 X - vector
94 */
FormExactSolution(DM da,AppCtx * user,Vec U)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;
102
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 }
116
ZeroBCSolution(AppCtx * user,const DMDACoor2d * c,PetscScalar * u)117 static PetscErrorCode ZeroBCSolution(AppCtx *user, const DMDACoor2d *c, PetscScalar *u)
118 {
119 u[0] = 0.;
120 return PETSC_SUCCESS;
121 }
122
123 /* The functions below evaluate the MMS solution u(x,y) and associated forcing
124
125 f(x,y) = -u_xx - u_yy - lambda exp(u)
126
127 such that u(x,y) is an exact solution with f(x,y) as the right-hand side forcing term.
128 */
MMSSolution1(AppCtx * user,const DMDACoor2d * c,PetscScalar * u)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 }
MMSForcing1(AppCtx * user,const DMDACoor2d * c,PetscScalar * f)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 }
142
MMSSolution2(AppCtx * user,const DMDACoor2d * c,PetscScalar * u)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 }
MMSForcing2(AppCtx * user,const DMDACoor2d * c,PetscScalar * f)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 }
156
MMSSolution3(AppCtx * user,const DMDACoor2d * c,PetscScalar * u)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 }
MMSForcing3(AppCtx * user,const DMDACoor2d * c,PetscScalar * f)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 }
171
MMSSolution4(AppCtx * user,const DMDACoor2d * c,PetscScalar * u)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 }
MMSForcing4(AppCtx * user,const DMDACoor2d * c,PetscScalar * f)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 }
187
188 /* ------------------------------------------------------------------- */
189 /*
190 FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
191
192 */
FormFunctionLocal(DMDALocalInfo * info,PetscScalar ** x,PetscScalar ** f,AppCtx * user)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;
199
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];
222
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 }
244
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 }
258
259 /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
FormObjectiveLocal(DMDALocalInfo * info,PetscScalar ** x,PetscReal * obj,AppCtx * user)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;
266
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];
289
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.;
294
295 /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */
296
297 uxux = u * (2. * u - ue - uw) * hydhx;
298 uyuy = u * (2. * u - un - us) * hxdhy;
299
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 }
308
309 /*
310 FormJacobianLocal - Evaluates Jacobian matrix on local process patch
311 */
FormJacobianLocal(DMDALocalInfo * info,PetscScalar ** x,Mat jac,Mat jacpre,AppCtx * user)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;
320
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;
333
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 }
368
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++;
373
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 }
390
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));
397
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 }
405
FormFunctionMatlab(SNES snes,Vec X,Vec F,void * ptr)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;
415
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));
423
424 lambda = user->param;
425 hx = 1.0 / (PetscReal)(Mx - 1);
426 hy = 1.0 / (PetscReal)(My - 1);
427
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));
440
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 }
453
454 /* ------------------------------------------------------------------- */
455 /*
456 Applies some sweeps on nonlinear Gauss-Seidel on each process
457
458 */
NonlinearGS(SNES snes,Vec X,Vec B,PetscCtx ctx)459 static PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, PetscCtx 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;
468
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));
475
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));
477
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;
484
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));
509
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.;
518
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];
524
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++;
535
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 }
557
main(int argc,char ** argv)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;
573
574 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
575 Initialize program
576 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
577
578 PetscFunctionBeginUser;
579 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
580
581 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_kokkos", &useKokkos, NULL));
582
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 }
598
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));
605
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));
624
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 }
652
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 }
659
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 }
666
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 }
673
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 }
682
683 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
684 Customize nonlinear solver; set runtime options
685 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
686 PetscCall(SNESSetFromOptions(snes));
687
688 PetscCall(FormInitialGuess(da, &user, x));
689
690 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
691 Solve nonlinear system
692 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
693 PetscCall(SNESSolve(snes, NULL, x));
694 PetscCall(SNESGetIterationNumber(snes, &its));
695
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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
702
703 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
704 If using MMS, check the l_2 error
705 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
706 if (setMMS) {
707 Vec e;
708 PetscReal errorl2, errorinf;
709 PetscInt N;
710
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 }
725
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 }
737
738 /*TEST
739 build:
740 requires: !windows_compilers
741 depends: ex55k.kokkos.cxx
742
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"
748
749 test:
750 suffix: asm_0
751 args: -fd {{0 1}}
752
753 test:
754 requires: kokkos_kernels
755 suffix: asm_0_kok
756 args: -use_kokkos -fd {{0 1}}
757
758 testset:
759 output_file: output/ex55_1.out
760 requires: !single
761 args: -snes_monitor
762 filter: grep -v "type"
763
764 test:
765 suffix: 1
766 args: -fd {{0 1}}
767
768 test:
769 requires: kokkos_kernels
770 suffix: 1_kok
771 args: -use_kokkos -fd {{0 1}}
772
773 test:
774 requires: h2opus
775 suffix: 1_h2opus
776 args: -pc_type h2opus -fd {{0 1}}
777
778 test:
779 requires: h2opus kokkos_kernels
780 suffix: 1_h2opus_k
781 args: -use_kokkos -pc_type h2opus -fd {{0 1}}
782
783 TEST*/
784