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 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, (char *)0, 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 */ 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 */ 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 */ 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 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 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 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