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