1 #include <petscsnes.h>
2 #include <petscdmda.h>
3
4 static const char help[] = "Minimum surface area problem in 2D using DMDA.\n\
5 It solves an unconstrained minimization problem. This example is based on a \n\
6 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
7 boundary values along the edges of the domain, the objective is to find the\n\
8 surface with the minimal area that satisfies the boundary conditions.\n\
9 \n\
10 The command line options are:\n\
11 -da_grid_x <nx>, where <nx> = number of grid points in the 1st coordinate direction\n\
12 -da_grid_y <ny>, where <ny> = number of grid points in the 2nd coordinate direction\n\
13 \n";
14
15 /*
16 User-defined application context - contains data needed by the
17 application-provided call-back routines, FormJacobianLocal() and
18 FormFunctionLocal().
19 */
20
21 typedef enum {
22 PROBLEM_ENNEPER,
23 PROBLEM_SINS,
24 } ProblemType;
25 static const char *const ProblemTypes[] = {"ENNEPER", "SINS", "ProblemType", "PROBLEM_", 0};
26
27 typedef struct {
28 PetscScalar *bottom, *top, *left, *right;
29 } AppCtx;
30
31 /* -------- User-defined Routines --------- */
32
33 static PetscErrorCode FormBoundaryConditions_Enneper(SNES, AppCtx **);
34 static PetscErrorCode FormBoundaryConditions_Sins(SNES, AppCtx **);
35 static PetscErrorCode DestroyBoundaryConditions(AppCtx **);
36 static PetscErrorCode FormObjectiveLocal(DMDALocalInfo *, PetscScalar **, PetscReal *, void *);
37 static PetscErrorCode FormFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, void *);
38 static PetscErrorCode FormJacobianLocal(DMDALocalInfo *, PetscScalar **, Mat, Mat, void *);
39
main(int argc,char ** argv)40 int main(int argc, char **argv)
41 {
42 Vec x;
43 SNES snes;
44 DM da;
45 ProblemType ptype = PROBLEM_ENNEPER;
46 PetscBool use_obj = PETSC_TRUE;
47 PetscReal bbox[4] = {0.};
48 AppCtx *user;
49 PetscErrorCode (*form_bc)(SNES, AppCtx **) = NULL;
50
51 PetscFunctionBeginUser;
52 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
53
54 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Minimal surface options", __FILE__);
55 PetscCall(PetscOptionsEnum("-problem_type", "Problem type", NULL, ProblemTypes, (PetscEnum)ptype, (PetscEnum *)&ptype, NULL));
56 PetscCall(PetscOptionsBool("-use_objective", "Use objective function", NULL, use_obj, &use_obj, NULL));
57 PetscOptionsEnd();
58 switch (ptype) {
59 case PROBLEM_ENNEPER:
60 bbox[0] = -0.5;
61 bbox[1] = 0.5;
62 bbox[2] = -0.5;
63 bbox[3] = 0.5;
64 form_bc = FormBoundaryConditions_Enneper;
65 break;
66 case PROBLEM_SINS:
67 bbox[0] = 0.0;
68 bbox[1] = 1.0;
69 bbox[2] = 0.0;
70 bbox[3] = 1.0;
71 form_bc = FormBoundaryConditions_Sins;
72 break;
73 }
74
75 /* Create distributed array to manage the 2d grid */
76 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));
77 PetscCall(DMSetFromOptions(da));
78 PetscCall(DMSetUp(da));
79 PetscCall(DMDASetUniformCoordinates(da, bbox[0], bbox[1], bbox[2], bbox[3], PETSC_DECIDE, PETSC_DECIDE));
80
81 /* Extract global vectors from DMDA; */
82 PetscCall(DMCreateGlobalVector(da, &x));
83
84 /* Create nonlinear solver context */
85 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
86 PetscCall(SNESSetDM(snes, da));
87 PetscCall((*form_bc)(snes, &user));
88 PetscCall(SNESSetApplicationContext(snes, user));
89
90 /* Set local callbacks */
91 if (use_obj) PetscCall(DMDASNESSetObjectiveLocal(da, (DMDASNESObjectiveFn *)FormObjectiveLocal, NULL));
92 PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, NULL));
93 PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, NULL));
94
95 /* Customize from command line */
96 PetscCall(SNESSetFromOptions(snes));
97
98 /* Solve the application */
99 PetscCall(SNESSolve(snes, NULL, x));
100
101 /* Free user-created data structures */
102 PetscCall(VecDestroy(&x));
103 PetscCall(SNESDestroy(&snes));
104 PetscCall(DMDestroy(&da));
105 PetscCall(DestroyBoundaryConditions(&user));
106
107 PetscCall(PetscFinalize());
108 return 0;
109 }
110
111 /* Compute objective function over the locally owned part of the mesh */
FormObjectiveLocal(DMDALocalInfo * info,PetscScalar ** x,PetscReal * v,void * ptr)112 PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info, PetscScalar **x, PetscReal *v, void *ptr)
113 {
114 AppCtx *user = (AppCtx *)ptr;
115 PetscInt mx = info->mx, my = info->my;
116 PetscInt xs = info->xs, xm = info->xm, ys = info->ys, ym = info->ym;
117 PetscInt i, j;
118 PetscScalar hx, hy;
119 PetscScalar f2, f4, d1, d2, d3, d4, xc, xl, xr, xt, xb;
120 PetscReal ft = 0, area;
121
122 PetscFunctionBeginUser;
123 PetscCheck(user, PetscObjectComm((PetscObject)info->da), PETSC_ERR_PLIB, "Missing application context");
124 hx = 1.0 / (mx + 1);
125 hy = 1.0 / (my + 1);
126 area = 0.5 * hx * hy;
127 for (j = ys; j < ys + ym; j++) {
128 for (i = xs; i < xs + xm; i++) {
129 xc = x[j][i];
130 xl = xr = xb = xt = xc;
131
132 if (i == 0) { /* left side */
133 xl = user->left[j + 1];
134 } else xl = x[j][i - 1];
135
136 if (j == 0) { /* bottom side */
137 xb = user->bottom[i + 1];
138 } else xb = x[j - 1][i];
139
140 if (i + 1 == mx) { /* right side */
141 xr = user->right[j + 1];
142 } else xr = x[j][i + 1];
143
144 if (j + 1 == 0 + my) { /* top side */
145 xt = user->top[i + 1];
146 } else xt = x[j + 1][i];
147
148 d1 = (xc - xl);
149 d2 = (xc - xr);
150 d3 = (xc - xt);
151 d4 = (xc - xb);
152
153 d1 /= hx;
154 d2 /= hx;
155 d3 /= hy;
156 d4 /= hy;
157
158 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
159 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
160
161 ft += PetscRealPart(f2 + f4);
162 }
163 }
164
165 /* Compute triangular areas along the border of the domain. */
166 if (xs == 0) { /* left side */
167 for (j = ys; j < ys + ym; j++) {
168 d3 = (user->left[j + 1] - user->left[j + 2]) / hy;
169 d2 = (user->left[j + 1] - x[j][0]) / hx;
170 ft += PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
171 }
172 }
173 if (ys == 0) { /* bottom side */
174 for (i = xs; i < xs + xm; i++) {
175 d2 = (user->bottom[i + 1] - user->bottom[i + 2]) / hx;
176 d3 = (user->bottom[i + 1] - x[0][i]) / hy;
177 ft += PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
178 }
179 }
180 if (xs + xm == mx) { /* right side */
181 for (j = ys; j < ys + ym; j++) {
182 d1 = (x[j][mx - 1] - user->right[j + 1]) / hx;
183 d4 = (user->right[j] - user->right[j + 1]) / hy;
184 ft += PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
185 }
186 }
187 if (ys + ym == my) { /* top side */
188 for (i = xs; i < xs + xm; i++) {
189 d1 = (x[my - 1][i] - user->top[i + 1]) / hy;
190 d4 = (user->top[i + 1] - user->top[i]) / hx;
191 ft += PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
192 }
193 }
194 if (ys == 0 && xs == 0) {
195 d1 = (user->left[0] - user->left[1]) / hy;
196 d2 = (user->bottom[0] - user->bottom[1]) / hx;
197 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
198 }
199 if (ys + ym == my && xs + xm == mx) {
200 d1 = (user->right[ym + 1] - user->right[ym]) / hy;
201 d2 = (user->top[xm + 1] - user->top[xm]) / hx;
202 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
203 }
204 ft *= area;
205 *v = ft;
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
209 /* Compute gradient over the locally owned part of the mesh */
FormFunctionLocal(DMDALocalInfo * info,PetscScalar ** x,PetscScalar ** g,void * ptr)210 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **g, void *ptr)
211 {
212 AppCtx *user = (AppCtx *)ptr;
213 PetscInt mx = info->mx, my = info->my;
214 PetscInt xs = info->xs, xm = info->xm, ys = info->ys, ym = info->ym;
215 PetscInt i, j;
216 PetscScalar hx, hy, hydhx, hxdhy;
217 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
218 PetscScalar df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
219
220 PetscFunctionBeginUser;
221 PetscCheck(user, PetscObjectComm((PetscObject)info->da), PETSC_ERR_PLIB, "Missing application context");
222 hx = 1.0 / (mx + 1);
223 hy = 1.0 / (my + 1);
224 hydhx = hy / hx;
225 hxdhy = hx / hy;
226
227 for (j = ys; j < ys + ym; j++) {
228 for (i = xs; i < xs + xm; i++) {
229 xc = x[j][i];
230 xlt = xrb = xl = xr = xb = xt = xc;
231
232 if (i == 0) { /* left side */
233 xl = user->left[j + 1];
234 xlt = user->left[j + 2];
235 } else xl = x[j][i - 1];
236
237 if (j == 0) { /* bottom side */
238 xb = user->bottom[i + 1];
239 xrb = user->bottom[i + 2];
240 } else xb = x[j - 1][i];
241
242 if (i + 1 == mx) { /* right side */
243 xr = user->right[j + 1];
244 xrb = user->right[j];
245 } else xr = x[j][i + 1];
246
247 if (j + 1 == 0 + my) { /* top side */
248 xt = user->top[i + 1];
249 xlt = user->top[i];
250 } else xt = x[j + 1][i];
251
252 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; /* left top side */
253 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; /* right bottom */
254
255 d1 = (xc - xl);
256 d2 = (xc - xr);
257 d3 = (xc - xt);
258 d4 = (xc - xb);
259 d5 = (xr - xrb);
260 d6 = (xrb - xb);
261 d7 = (xlt - xl);
262 d8 = (xt - xlt);
263
264 df1dxc = d1 * hydhx;
265 df2dxc = (d1 * hydhx + d4 * hxdhy);
266 df3dxc = d3 * hxdhy;
267 df4dxc = (d2 * hydhx + d3 * hxdhy);
268 df5dxc = d2 * hydhx;
269 df6dxc = d4 * hxdhy;
270
271 d1 /= hx;
272 d2 /= hx;
273 d3 /= hy;
274 d4 /= hy;
275 d5 /= hy;
276 d6 /= hx;
277 d7 /= hy;
278 d8 /= hx;
279
280 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
281 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
282 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
283 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
284 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
285 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
286
287 df1dxc /= f1;
288 df2dxc /= f2;
289 df3dxc /= f3;
290 df4dxc /= f4;
291 df5dxc /= f5;
292 df6dxc /= f6;
293
294 g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
295 }
296 }
297 PetscFunctionReturn(PETSC_SUCCESS);
298 }
299
300 /* Compute Hessian over the locally owned part of the mesh */
FormJacobianLocal(DMDALocalInfo * info,PetscScalar ** x,Mat H,Mat Hp,void * ptr)301 PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **x, Mat H, Mat Hp, void *ptr)
302 {
303 AppCtx *user = (AppCtx *)ptr;
304 PetscInt mx = info->mx, my = info->my;
305 PetscInt xs = info->xs, xm = info->xm, ys = info->ys, ym = info->ym;
306 PetscInt i, j, k;
307 MatStencil row, col[7];
308 PetscScalar hx, hy, hydhx, hxdhy;
309 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
310 PetscScalar hl, hr, ht, hb, hc, htl, hbr;
311 PetscScalar v[7];
312
313 PetscFunctionBeginUser;
314 PetscCheck(user, PetscObjectComm((PetscObject)info->da), PETSC_ERR_PLIB, "Missing application context");
315 hx = 1.0 / (mx + 1);
316 hy = 1.0 / (my + 1);
317 hydhx = hy / hx;
318 hxdhy = hx / hy;
319
320 for (j = ys; j < ys + ym; j++) {
321 for (i = xs; i < xs + xm; i++) {
322 xc = x[j][i];
323 xlt = xrb = xl = xr = xb = xt = xc;
324
325 /* Left */
326 if (i == 0) {
327 xl = user->left[j + 1];
328 xlt = user->left[j + 2];
329 } else xl = x[j][i - 1];
330
331 /* Bottom */
332 if (j == 0) {
333 xb = user->bottom[i + 1];
334 xrb = user->bottom[i + 2];
335 } else xb = x[j - 1][i];
336
337 /* Right */
338 if (i + 1 == mx) {
339 xr = user->right[j + 1];
340 xrb = user->right[j];
341 } else xr = x[j][i + 1];
342
343 /* Top */
344 if (j + 1 == my) {
345 xt = user->top[i + 1];
346 xlt = user->top[i];
347 } else xt = x[j + 1][i];
348
349 /* Top left */
350 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
351
352 /* Bottom right */
353 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
354
355 d1 = (xc - xl) / hx;
356 d2 = (xc - xr) / hx;
357 d3 = (xc - xt) / hy;
358 d4 = (xc - xb) / hy;
359 d5 = (xrb - xr) / hy;
360 d6 = (xrb - xb) / hx;
361 d7 = (xlt - xl) / hy;
362 d8 = (xlt - xt) / hx;
363
364 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
365 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
366 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
367 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
368 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
369 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
370
371 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
372 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
373 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
374 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
375
376 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
377 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
378
379 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);
380
381 hl /= 2.0;
382 hr /= 2.0;
383 ht /= 2.0;
384 hb /= 2.0;
385 hbr /= 2.0;
386 htl /= 2.0;
387 hc /= 2.0;
388
389 k = 0;
390 row.i = i;
391 row.j = j;
392 /* Bottom */
393 if (j > 0) {
394 v[k] = hb;
395 col[k].i = i;
396 col[k].j = j - 1;
397 k++;
398 }
399
400 /* Bottom right */
401 if (j > 0 && i < mx - 1) {
402 v[k] = hbr;
403 col[k].i = i + 1;
404 col[k].j = j - 1;
405 k++;
406 }
407
408 /* left */
409 if (i > 0) {
410 v[k] = hl;
411 col[k].i = i - 1;
412 col[k].j = j;
413 k++;
414 }
415
416 /* Centre */
417 v[k] = hc;
418 col[k].i = row.i;
419 col[k].j = row.j;
420 k++;
421
422 /* Right */
423 if (i < mx - 1) {
424 v[k] = hr;
425 col[k].i = i + 1;
426 col[k].j = j;
427 k++;
428 }
429
430 /* Top left */
431 if (i > 0 && j < my - 1) {
432 v[k] = htl;
433 col[k].i = i - 1;
434 col[k].j = j + 1;
435 k++;
436 }
437
438 /* Top */
439 if (j < my - 1) {
440 v[k] = ht;
441 col[k].i = i;
442 col[k].j = j + 1;
443 k++;
444 }
445
446 PetscCall(MatSetValuesStencil(Hp, 1, &row, k, col, v, INSERT_VALUES));
447 }
448 }
449
450 /* Assemble the matrix */
451 PetscCall(MatAssemblyBegin(Hp, MAT_FINAL_ASSEMBLY));
452 PetscCall(MatAssemblyEnd(Hp, MAT_FINAL_ASSEMBLY));
453 PetscFunctionReturn(PETSC_SUCCESS);
454 }
455
FormBoundaryConditions_Enneper(SNES snes,AppCtx ** ouser)456 PetscErrorCode FormBoundaryConditions_Enneper(SNES snes, AppCtx **ouser)
457 {
458 PetscInt i, j, k, limit = 0, maxits = 5;
459 PetscInt mx, my;
460 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
461 PetscScalar one = 1.0, two = 2.0, three = 3.0;
462 PetscScalar det, hx, hy, xt = 0, yt = 0;
463 PetscReal fnorm, tol = 1e-10;
464 PetscScalar u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
465 PetscScalar b = -0.5, t = 0.5, l = -0.5, r = 0.5;
466 PetscScalar *boundary;
467 AppCtx *user;
468 DM da;
469
470 PetscFunctionBeginUser;
471 PetscCall(SNESGetDM(snes, &da));
472 PetscCall(PetscNew(&user));
473 *ouser = user;
474 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));
475 bsize = mx + 2;
476 lsize = my + 2;
477 rsize = my + 2;
478 tsize = mx + 2;
479
480 PetscCall(PetscMalloc1(bsize, &user->bottom));
481 PetscCall(PetscMalloc1(tsize, &user->top));
482 PetscCall(PetscMalloc1(lsize, &user->left));
483 PetscCall(PetscMalloc1(rsize, &user->right));
484
485 hx = 1.0 / (mx + 1.0);
486 hy = 1.0 / (my + 1.0);
487
488 for (j = 0; j < 4; j++) {
489 if (j == 0) {
490 yt = b;
491 xt = l;
492 limit = bsize;
493 boundary = user->bottom;
494 } else if (j == 1) {
495 yt = t;
496 xt = l;
497 limit = tsize;
498 boundary = user->top;
499 } else if (j == 2) {
500 yt = b;
501 xt = l;
502 limit = lsize;
503 boundary = user->left;
504 } else { /* if (j==3) */
505 yt = b;
506 xt = r;
507 limit = rsize;
508 boundary = user->right;
509 }
510
511 for (i = 0; i < limit; i++) {
512 u1 = xt;
513 u2 = -yt;
514 for (k = 0; k < maxits; k++) {
515 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
516 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
517 fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2));
518 if (fnorm <= tol) break;
519 njac11 = one + u2 * u2 - u1 * u1;
520 njac12 = two * u1 * u2;
521 njac21 = -two * u1 * u2;
522 njac22 = -one - u1 * u1 + u2 * u2;
523 det = njac11 * njac22 - njac21 * njac12;
524 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
525 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
526 }
527
528 boundary[i] = u1 * u1 - u2 * u2;
529 if (j == 0 || j == 1) xt = xt + hx;
530 else yt = yt + hy; /* if (j==2 || j==3) */
531 }
532 }
533 PetscFunctionReturn(PETSC_SUCCESS);
534 }
535
DestroyBoundaryConditions(AppCtx ** ouser)536 PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser)
537 {
538 AppCtx *user = *ouser;
539
540 PetscFunctionBeginUser;
541 PetscCall(PetscFree(user->bottom));
542 PetscCall(PetscFree(user->top));
543 PetscCall(PetscFree(user->left));
544 PetscCall(PetscFree(user->right));
545 PetscCall(PetscFree(*ouser));
546 PetscFunctionReturn(PETSC_SUCCESS);
547 }
548
FormBoundaryConditions_Sins(SNES snes,AppCtx ** ouser)549 PetscErrorCode FormBoundaryConditions_Sins(SNES snes, AppCtx **ouser)
550 {
551 PetscInt i, j;
552 PetscInt mx, my;
553 PetscInt limit, bsize = 0, lsize = 0, tsize = 0, rsize = 0;
554 PetscScalar hx, hy, xt = 0, yt = 0;
555 PetscScalar b = 0.0, t = 1.0, l = 0.0, r = 1.0;
556 PetscScalar *boundary;
557 AppCtx *user;
558 DM da;
559 PetscReal pi2 = 2 * PETSC_PI;
560
561 PetscFunctionBeginUser;
562 PetscCall(SNESGetDM(snes, &da));
563 PetscCall(PetscNew(&user));
564 *ouser = user;
565 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));
566 bsize = mx + 2;
567 lsize = my + 2;
568 rsize = my + 2;
569 tsize = mx + 2;
570
571 PetscCall(PetscMalloc1(bsize, &user->bottom));
572 PetscCall(PetscMalloc1(tsize, &user->top));
573 PetscCall(PetscMalloc1(lsize, &user->left));
574 PetscCall(PetscMalloc1(rsize, &user->right));
575
576 hx = 1.0 / (mx + 1.0);
577 hy = 1.0 / (my + 1.0);
578
579 for (j = 0; j < 4; j++) {
580 if (j == 0) {
581 yt = b;
582 xt = l;
583 limit = bsize;
584 boundary = user->bottom;
585 } else if (j == 1) {
586 yt = t;
587 xt = l;
588 limit = tsize;
589 boundary = user->top;
590 } else if (j == 2) {
591 yt = b;
592 xt = l;
593 limit = lsize;
594 boundary = user->left;
595 } else { /* if (j==3) */
596 yt = b;
597 xt = r;
598 limit = rsize;
599 boundary = user->right;
600 }
601
602 for (i = 0; i < limit; i++) {
603 if (j == 0) { /* bottom */
604 boundary[i] = -0.5 * PetscSinReal(pi2 * xt);
605 } else if (j == 1) { /* top */
606 boundary[i] = 0.5 * PetscSinReal(pi2 * xt);
607 } else if (j == 2) { /* left */
608 boundary[i] = -0.5 * PetscSinReal(pi2 * yt);
609 } else { /* right */
610 boundary[i] = 0.5 * PetscSinReal(pi2 * yt);
611 }
612 if (j == 0 || j == 1) xt = xt + hx;
613 else yt = yt + hy;
614 }
615 }
616 PetscFunctionReturn(PETSC_SUCCESS);
617 }
618
619 /*TEST
620
621 build:
622 requires: !complex
623
624 test:
625 requires: !single
626 filter: sed -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"
627 suffix: qn_nasm
628 args: -snes_type qn -snes_npc_side {{left right}separate output} -npc_snes_type nasm -snes_converged_reason -da_local_subdomains 4 -use_objective {{0 1}separate output}
629
630 TEST*/
631