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