1 #include <petscsnes.h>
2 #include <petscdm.h>
3 #include <petscdmda.h>
4
5 static const char help[] = "Parallel version of the minimum surface area problem in 2D using DMDA.\n\
6 It solves a system of nonlinear equations in mixed\n\
7 complementarity form.This example is based on a\n\
8 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\
9 boundary values along the edges of the domain, the objective is to find the\n\
10 surface with the minimal area that satisfies the boundary conditions.\n\
11 This application solves this problem using complementarity -- We are actually\n\
12 solving the system (grad f)_i >= 0, if x_i == l_i \n\
13 (grad f)_i = 0, if l_i < x_i < u_i \n\
14 (grad f)_i <= 0, if x_i == u_i \n\
15 where f is the function to be minimized. \n\
16 \n\
17 The command line options are:\n\
18 -da_grid_x <nx>, where <nx> = number of grid points in the 1st coordinate direction\n\
19 -da_grid_y <ny>, where <ny> = number of grid points in the 2nd coordinate direction\n\
20 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise\n\
21 -lb <value>, lower bound on the variables\n\
22 -ub <value>, upper bound on the variables\n\n";
23
24 /*
25 User-defined application context - contains data needed by the
26 application-provided call-back routines, FormJacobian() and
27 FormFunction().
28 */
29
30 /*
31 Run, for example, with the options ./ex58 -snes_vi_monitor -ksp_monitor -mg_levels_ksp_monitor -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin pmat -ksp_type fgmres
32
33 Or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
34 multigrid levels, it will be determined automatically based on the number of refinements done)
35
36 ./ex58 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3
37 -mg_levels_ksp_monitor -snes_vi_monitor -mg_levels_pc_type sor -pc_mg_type full
38
39 */
40
41 typedef struct {
42 PetscScalar *bottom, *top, *left, *right;
43 PetscScalar lb, ub;
44 } AppCtx;
45
46 /* -------- User-defined Routines --------- */
47
48 extern PetscErrorCode FormBoundaryConditions(SNES, PetscCtxRt);
49 extern PetscErrorCode DestroyBoundaryConditions(PetscCtxRt);
50 extern PetscErrorCode ComputeInitialGuess(SNES, Vec, void *);
51 extern PetscErrorCode FormGradient(SNES, Vec, Vec, void *);
52 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
53 extern PetscErrorCode FormBounds(SNES, Vec, Vec);
54
main(int argc,char ** argv)55 int main(int argc, char **argv)
56 {
57 Vec x, r; /* solution and residual vectors */
58 SNES snes; /* nonlinear solver context */
59 Mat J; /* Jacobian matrix */
60 DM da;
61
62 PetscFunctionBeginUser;
63 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
64
65 /* Create distributed array to manage the 2d grid */
66 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
67 PetscCall(DMSetFromOptions(da));
68 PetscCall(DMSetUp(da));
69
70 /* Extract global vectors from DMDA; */
71 PetscCall(DMCreateGlobalVector(da, &x));
72 PetscCall(VecDuplicate(x, &r));
73
74 PetscCall(DMSetMatType(da, MATAIJ));
75 PetscCall(DMCreateMatrix(da, &J));
76
77 /* Create nonlinear solver context */
78 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
79 PetscCall(SNESSetDM(snes, da));
80
81 /* Set function evaluation and Jacobian evaluation routines */
82 PetscCall(SNESSetFunction(snes, r, FormGradient, NULL));
83 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
84
85 PetscCall(SNESSetComputeApplicationContext(snes, FormBoundaryConditions, DestroyBoundaryConditions));
86
87 PetscCall(SNESSetComputeInitialGuess(snes, ComputeInitialGuess, NULL));
88
89 PetscCall(SNESVISetComputeVariableBounds(snes, FormBounds));
90
91 PetscCall(SNESSetFromOptions(snes));
92
93 /* Solve the application */
94 PetscCall(SNESSolve(snes, NULL, x));
95
96 /* Free memory */
97 PetscCall(VecDestroy(&x));
98 PetscCall(VecDestroy(&r));
99 PetscCall(MatDestroy(&J));
100 PetscCall(SNESDestroy(&snes));
101
102 /* Free user-created data structures */
103 PetscCall(DMDestroy(&da));
104
105 PetscCall(PetscFinalize());
106 return 0;
107 }
108
109 /* -------------------------------------------------------------------- */
110
111 /* FormBounds - sets the upper and lower bounds
112
113 Input Parameters:
114 . snes - the SNES context
115
116 Output Parameters:
117 . xl - lower bounds
118 . xu - upper bounds
119 */
FormBounds(SNES snes,Vec xl,Vec xu)120 PetscErrorCode FormBounds(SNES snes, Vec xl, Vec xu)
121 {
122 AppCtx *ctx;
123
124 PetscFunctionBeginUser;
125 PetscCall(SNESGetApplicationContext(snes, &ctx));
126 PetscCall(VecSet(xl, ctx->lb));
127 PetscCall(VecSet(xu, ctx->ub));
128 PetscFunctionReturn(PETSC_SUCCESS);
129 }
130
131 /* -------------------------------------------------------------------- */
132
133 /* FormGradient - Evaluates gradient of f.
134
135 Input Parameters:
136 . snes - the SNES context
137 . X - input vector
138 . unused - optional user-defined context, as set by SNESSetFunction()
139
140 Output Parameters:
141 . G - vector containing the newly evaluated gradient
142 */
FormGradient(SNES snes,Vec X,Vec G,void * unused)143 PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *unused)
144 {
145 AppCtx *ctx;
146 PetscInt i, j;
147 PetscInt mx, my;
148 PetscScalar hx, hy, hydhx, hxdhy;
149 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
150 PetscScalar df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
151 PetscScalar **g, **x;
152 PetscInt xs, xm, ys, ym;
153 Vec localX;
154 DM da;
155
156 PetscFunctionBeginUser;
157 PetscCall(SNESGetDM(snes, &da));
158 PetscCall(SNESGetApplicationContext(snes, &ctx));
159 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));
160 hx = 1.0 / (mx + 1);
161 hy = 1.0 / (my + 1);
162 hydhx = hy / hx;
163 hxdhy = hx / hy;
164
165 PetscCall(VecSet(G, 0.0));
166
167 /* Get local vector */
168 PetscCall(DMGetLocalVector(da, &localX));
169 /* Get ghost points */
170 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
171 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
172 /* Get pointer to local vector data */
173 PetscCall(DMDAVecGetArray(da, localX, &x));
174 PetscCall(DMDAVecGetArray(da, G, &g));
175
176 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
177 /* Compute function over the locally owned part of the mesh */
178 for (j = ys; j < ys + ym; j++) {
179 for (i = xs; i < xs + xm; i++) {
180 xc = x[j][i];
181 xlt = xrb = xl = xr = xb = xt = xc;
182
183 if (i == 0) { /* left side */
184 xl = ctx->left[j + 1];
185 xlt = ctx->left[j + 2];
186 } else xl = x[j][i - 1];
187
188 if (j == 0) { /* bottom side */
189 xb = ctx->bottom[i + 1];
190 xrb = ctx->bottom[i + 2];
191 } else xb = x[j - 1][i];
192
193 if (i + 1 == mx) { /* right side */
194 xr = ctx->right[j + 1];
195 xrb = ctx->right[j];
196 } else xr = x[j][i + 1];
197
198 if (j + 1 == 0 + my) { /* top side */
199 xt = ctx->top[i + 1];
200 xlt = ctx->top[i];
201 } else xt = x[j + 1][i];
202
203 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; /* left top side */
204 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; /* right bottom */
205
206 d1 = (xc - xl);
207 d2 = (xc - xr);
208 d3 = (xc - xt);
209 d4 = (xc - xb);
210 d5 = (xr - xrb);
211 d6 = (xrb - xb);
212 d7 = (xlt - xl);
213 d8 = (xt - xlt);
214
215 df1dxc = d1 * hydhx;
216 df2dxc = (d1 * hydhx + d4 * hxdhy);
217 df3dxc = d3 * hxdhy;
218 df4dxc = (d2 * hydhx + d3 * hxdhy);
219 df5dxc = d2 * hydhx;
220 df6dxc = d4 * hxdhy;
221
222 d1 /= hx;
223 d2 /= hx;
224 d3 /= hy;
225 d4 /= hy;
226 d5 /= hy;
227 d6 /= hx;
228 d7 /= hy;
229 d8 /= hx;
230
231 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
232 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
233 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
234 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
235 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
236 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
237
238 df1dxc /= f1;
239 df2dxc /= f2;
240 df3dxc /= f3;
241 df4dxc /= f4;
242 df5dxc /= f5;
243 df6dxc /= f6;
244
245 g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
246 }
247 }
248
249 /* Restore vectors */
250 PetscCall(DMDAVecRestoreArray(da, localX, &x));
251 PetscCall(DMDAVecRestoreArray(da, G, &g));
252 PetscCall(DMRestoreLocalVector(da, &localX));
253 PetscCall(PetscLogFlops(67.0 * mx * my));
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
257 /* ------------------------------------------------------------------- */
258 /*
259 FormJacobian - Evaluates Jacobian matrix.
260
261 Input Parameters:
262 . snes - SNES context
263 . X - input vector
264 . unused - optional user-defined context, as set by SNESSetJacobian()
265
266 Output Parameters:
267 . tH - Jacobian matrix
268
269 */
FormJacobian(SNES snes,Vec X,Mat H,Mat tHPre,void * unused)270 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat H, Mat tHPre, void *unused)
271 {
272 AppCtx *ctx;
273 PetscInt i, j, k;
274 PetscInt mx, my;
275 MatStencil row, col[7];
276 PetscScalar hx, hy, hydhx, hxdhy;
277 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
278 PetscScalar hl, hr, ht, hb, hc, htl, hbr;
279 PetscScalar **x, v[7];
280 PetscBool assembled;
281 PetscInt xs, xm, ys, ym;
282 Vec localX;
283 DM da;
284
285 PetscFunctionBeginUser;
286 PetscCall(SNESGetDM(snes, &da));
287 PetscCall(SNESGetApplicationContext(snes, &ctx));
288 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));
289 hx = 1.0 / (mx + 1);
290 hy = 1.0 / (my + 1);
291 hydhx = hy / hx;
292 hxdhy = hx / hy;
293
294 /* Set various matrix options */
295 PetscCall(MatAssembled(H, &assembled));
296 if (assembled) PetscCall(MatZeroEntries(H));
297
298 /* Get local vector */
299 PetscCall(DMGetLocalVector(da, &localX));
300 /* Get ghost points */
301 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
302 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
303
304 /* Get pointers to vector data */
305 PetscCall(DMDAVecGetArray(da, localX, &x));
306
307 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
308 /* Compute Jacobian over the locally owned part of the mesh */
309 for (j = ys; j < ys + ym; j++) {
310 for (i = xs; i < xs + xm; i++) {
311 xc = x[j][i];
312 xlt = xrb = xl = xr = xb = xt = xc;
313
314 /* Left */
315 if (i == 0) {
316 xl = ctx->left[j + 1];
317 xlt = ctx->left[j + 2];
318 } else xl = x[j][i - 1];
319
320 /* Bottom */
321 if (j == 0) {
322 xb = ctx->bottom[i + 1];
323 xrb = ctx->bottom[i + 2];
324 } else xb = x[j - 1][i];
325
326 /* Right */
327 if (i + 1 == mx) {
328 xr = ctx->right[j + 1];
329 xrb = ctx->right[j];
330 } else xr = x[j][i + 1];
331
332 /* Top */
333 if (j + 1 == my) {
334 xt = ctx->top[i + 1];
335 xlt = ctx->top[i];
336 } else xt = x[j + 1][i];
337
338 /* Top left */
339 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
340
341 /* Bottom right */
342 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
343
344 d1 = (xc - xl) / hx;
345 d2 = (xc - xr) / hx;
346 d3 = (xc - xt) / hy;
347 d4 = (xc - xb) / hy;
348 d5 = (xrb - xr) / hy;
349 d6 = (xrb - xb) / hx;
350 d7 = (xlt - xl) / hy;
351 d8 = (xlt - xt) / hx;
352
353 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
354 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
355 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
356 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
357 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
358 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
359
360 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
361 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
362 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
363 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
364
365 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
366 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
367
368 hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2.0 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2.0 * d2 * d3) / (f4 * f4 * f4);
369
370 hl /= 2.0;
371 hr /= 2.0;
372 ht /= 2.0;
373 hb /= 2.0;
374 hbr /= 2.0;
375 htl /= 2.0;
376 hc /= 2.0;
377
378 k = 0;
379 row.i = i;
380 row.j = j;
381 /* Bottom */
382 if (j > 0) {
383 v[k] = hb;
384 col[k].i = i;
385 col[k].j = j - 1;
386 k++;
387 }
388
389 /* Bottom right */
390 if (j > 0 && i < mx - 1) {
391 v[k] = hbr;
392 col[k].i = i + 1;
393 col[k].j = j - 1;
394 k++;
395 }
396
397 /* left */
398 if (i > 0) {
399 v[k] = hl;
400 col[k].i = i - 1;
401 col[k].j = j;
402 k++;
403 }
404
405 /* Centre */
406 v[k] = hc;
407 col[k].i = row.i;
408 col[k].j = row.j;
409 k++;
410
411 /* Right */
412 if (i < mx - 1) {
413 v[k] = hr;
414 col[k].i = i + 1;
415 col[k].j = j;
416 k++;
417 }
418
419 /* Top left */
420 if (i > 0 && j < my - 1) {
421 v[k] = htl;
422 col[k].i = i - 1;
423 col[k].j = j + 1;
424 k++;
425 }
426
427 /* Top */
428 if (j < my - 1) {
429 v[k] = ht;
430 col[k].i = i;
431 col[k].j = j + 1;
432 k++;
433 }
434
435 PetscCall(MatSetValuesStencil(H, 1, &row, k, col, v, INSERT_VALUES));
436 }
437 }
438
439 /* Assemble the matrix */
440 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
441 PetscCall(DMDAVecRestoreArray(da, localX, &x));
442 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
443 PetscCall(DMRestoreLocalVector(da, &localX));
444
445 PetscCall(PetscLogFlops(199.0 * mx * my));
446 PetscFunctionReturn(PETSC_SUCCESS);
447 }
448
449 /* ------------------------------------------------------------------- */
450 /*
451 FormBoundaryConditions - Calculates the boundary conditions for
452 the region.
453
454 Input Parameter:
455 . user - user-defined application context
456
457 Output Parameter:
458 . user - user-defined application context
459 */
FormBoundaryConditions(SNES snes,PetscCtxRt Ctx)460 PetscErrorCode FormBoundaryConditions(SNES snes, PetscCtxRt Ctx)
461 {
462 PetscInt i, j, k, limit = 0, maxits = 5;
463 PetscInt mx, my;
464 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
465 PetscScalar one = 1.0, two = 2.0, three = 3.0;
466 PetscScalar det, hx, hy, xt = 0, yt = 0;
467 PetscReal fnorm, tol = 1e-10;
468 PetscScalar u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
469 PetscScalar b = -0.5, t = 0.5, l = -0.5, r = 0.5;
470 PetscScalar *boundary;
471 AppCtx *ctx;
472 DM da;
473
474 PetscFunctionBeginUser;
475 PetscCall(SNESGetDM(snes, &da));
476 PetscCall(PetscNew(&ctx));
477 *(void **)Ctx = ctx;
478 ctx->lb = .05;
479 ctx->ub = PETSC_INFINITY;
480 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));
481
482 /* Check if lower and upper bounds are set */
483 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lb", &ctx->lb, 0));
484 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ub", &ctx->ub, 0));
485 bsize = mx + 2;
486 lsize = my + 2;
487 rsize = my + 2;
488 tsize = mx + 2;
489
490 PetscCall(PetscMalloc1(bsize, &ctx->bottom));
491 PetscCall(PetscMalloc1(tsize, &ctx->top));
492 PetscCall(PetscMalloc1(lsize, &ctx->left));
493 PetscCall(PetscMalloc1(rsize, &ctx->right));
494
495 hx = (r - l) / (mx + 1.0);
496 hy = (t - b) / (my + 1.0);
497
498 for (j = 0; j < 4; j++) {
499 if (j == 0) {
500 yt = b;
501 xt = l;
502 limit = bsize;
503 boundary = ctx->bottom;
504 } else if (j == 1) {
505 yt = t;
506 xt = l;
507 limit = tsize;
508 boundary = ctx->top;
509 } else if (j == 2) {
510 yt = b;
511 xt = l;
512 limit = lsize;
513 boundary = ctx->left;
514 } else { /* if (j==3) */
515 yt = b;
516 xt = r;
517 limit = rsize;
518 boundary = ctx->right;
519 }
520
521 for (i = 0; i < limit; i++) {
522 u1 = xt;
523 u2 = -yt;
524 for (k = 0; k < maxits; k++) {
525 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
526 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
527 fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2));
528 if (fnorm <= tol) break;
529 njac11 = one + u2 * u2 - u1 * u1;
530 njac12 = two * u1 * u2;
531 njac21 = -two * u1 * u2;
532 njac22 = -one - u1 * u1 + u2 * u2;
533 det = njac11 * njac22 - njac21 * njac12;
534 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
535 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
536 }
537
538 boundary[i] = u1 * u1 - u2 * u2;
539 if (j == 0 || j == 1) xt = xt + hx;
540 else yt = yt + hy; /* if (j==2 || j==3) */
541 }
542 }
543 PetscFunctionReturn(PETSC_SUCCESS);
544 }
545
DestroyBoundaryConditions(PetscCtxRt Ctx)546 PetscErrorCode DestroyBoundaryConditions(PetscCtxRt Ctx)
547 {
548 AppCtx *ctx = *(AppCtx **)Ctx;
549
550 PetscFunctionBeginUser;
551 PetscCall(PetscFree(ctx->bottom));
552 PetscCall(PetscFree(ctx->top));
553 PetscCall(PetscFree(ctx->left));
554 PetscCall(PetscFree(ctx->right));
555 PetscCall(PetscFree(ctx));
556 PetscFunctionReturn(PETSC_SUCCESS);
557 }
558
559 /* ------------------------------------------------------------------- */
560 /*
561 ComputeInitialGuess - Calculates the initial guess
562
563 Input Parameters:
564 . unused - user-defined application context
565 . X - vector for initial guess
566
567 Output Parameters:
568 . X - newly computed initial guess
569 */
ComputeInitialGuess(SNES snes,Vec X,void * unused)570 PetscErrorCode ComputeInitialGuess(SNES snes, Vec X, void *unused)
571 {
572 PetscInt i, j, mx, my;
573 DM da;
574 AppCtx *ctx;
575 PetscScalar **x;
576 PetscInt xs, xm, ys, ym;
577
578 PetscFunctionBeginUser;
579 PetscCall(SNESGetDM(snes, &da));
580 PetscCall(SNESGetApplicationContext(snes, &ctx));
581
582 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
583 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));
584
585 /* Get pointers to vector data */
586 PetscCall(DMDAVecGetArray(da, X, &x));
587 /* Perform local computations */
588 for (j = ys; j < ys + ym; j++) {
589 for (i = xs; i < xs + xm; i++) x[j][i] = (((j + 1.0) * ctx->bottom[i + 1] + (my - j + 1.0) * ctx->top[i + 1]) / (my + 2.0) + ((i + 1.0) * ctx->left[j + 1] + (mx - i + 1.0) * ctx->right[j + 1]) / (mx + 2.0)) / 2.0;
590 }
591 /* Restore vectors */
592 PetscCall(DMDAVecRestoreArray(da, X, &x));
593 PetscFunctionReturn(PETSC_SUCCESS);
594 }
595
596 /*TEST
597
598 test:
599 args: -snes_type vinewtonrsls -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason
600 requires: !single
601
602 test:
603 suffix: 2
604 args: -snes_type vinewtonssls -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason
605 requires: !single
606
607 test:
608 suffix: 3
609 args: -snes_type vinewtonrsls -snes_linesearch_type bisection -snes_linesearch_monitor -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason
610 requires: !single
611
612 test:
613 suffix: 4
614 args: -snes_type vinewtonssls -snes_linesearch_type bisection -snes_linesearch_monitor -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason
615 requires: !single
616
617 test:
618 suffix: 5
619 args: -snes_type vinewtonssls -snes_linesearch_type cp -snes_linesearch_order 3 -snes_linesearch_monitor -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason -snes_linesearch_maxlambda 2.0
620 requires: !single
621
622 test:
623 suffix: 6
624 args: -snes_type vinewtonssls -snes_linesearch_type secant -snes_linesearch_max_it 10 -snes_linesearch_ltol 1e-2 -snes_linesearch_atol 1e-32 -snes_linesearch_monitor -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason -snes_linesearch_maxlambda 1.0
625 requires: !single
626
627 TEST*/
628