1 static const char help[] = "p-Bratu nonlinear PDE in 2d.\n\ 2 We solve the p-Laplacian (nonlinear diffusion) combined with\n\ 3 the Bratu (solid fuel ignition) nonlinearity in a 2D rectangular\n\ 4 domain, using distributed arrays (DAs) to partition the parallel grid.\n\ 5 The command line options include:\n\ 6 -p <2>: `p' in p-Laplacian term\n\ 7 -epsilon <1e-05>: Strain-regularization in p-Laplacian\n\ 8 -lambda <6>: Bratu parameter\n\ 9 -blocks <bx,by>: number of coefficient interfaces in x and y direction\n\ 10 -kappa <1e-3>: diffusivity in odd regions\n\ 11 \n"; 12 13 /*F 14 The $p$-Bratu problem is a combination of the $p$-Laplacian (nonlinear diffusion) and the Brutu solid fuel ignition problem. 15 This problem is modeled by the partial differential equation 16 17 \begin{equation*} 18 -\nabla\cdot (\eta \nabla u) - \lambda \exp(u) = 0 19 \end{equation*} 20 21 on $\Omega = (-1,1)^2$ with closure 22 23 \begin{align*} 24 \eta(\gamma) &= (\epsilon^2 + \gamma)^{(p-2)/2} & \gamma &= \frac 1 2 |\nabla u|^2 25 \end{align*} 26 27 and boundary conditions $u = 0$ for $(x,y) \in \partial \Omega$ 28 29 A 9-point finite difference stencil is used to discretize 30 the boundary value problem to obtain a nonlinear system of equations. 31 This would be a 5-point stencil if not for the $p$-Laplacian's nonlinearity. 32 F*/ 33 34 /* 35 mpiexec -n 2 ./ex15 -snes_monitor -ksp_monitor log_summary 36 */ 37 38 /* 39 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 40 Include "petscsnes.h" so that we can use SNES solvers. Note that this 41 file automatically includes: 42 petsc.h - base PETSc routines petscvec.h - vectors 43 petscsys.h - system routines petscmat.h - matrices 44 petscis.h - index sets petscksp.h - Krylov subspace methods 45 petscviewer.h - viewers petscpc.h - preconditioners 46 petscksp.h - linear solvers 47 */ 48 49 #include <petscdm.h> 50 #include <petscdmda.h> 51 #include <petscsnes.h> 52 53 typedef enum {JAC_BRATU,JAC_PICARD,JAC_STAR,JAC_NEWTON} JacType; 54 static const char *const JacTypes[] = {"BRATU","PICARD","STAR","NEWTON","JacType","JAC_",0}; 55 56 /* 57 User-defined application context - contains data needed by the 58 application-provided call-back routines, FormJacobianLocal() and 59 FormFunctionLocal(). 60 */ 61 typedef struct { 62 PetscReal lambda; /* Bratu parameter */ 63 PetscReal p; /* Exponent in p-Laplacian */ 64 PetscReal epsilon; /* Regularization */ 65 PetscReal source; /* Source term */ 66 JacType jtype; /* What type of Jacobian to assemble */ 67 PetscBool picard; 68 PetscInt blocks[2]; 69 PetscReal kappa; 70 PetscInt initial; /* initial conditions type */ 71 } AppCtx; 72 73 /* 74 User-defined routines 75 */ 76 static PetscErrorCode FormRHS(AppCtx*,DM,Vec); 77 static PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec); 78 static PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*); 79 static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*); 80 static PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*); 81 static PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*); 82 83 typedef struct _n_PreCheck *PreCheck; 84 struct _n_PreCheck { 85 MPI_Comm comm; 86 PetscReal angle; 87 Vec Ylast; 88 PetscViewer monitor; 89 }; 90 PetscErrorCode PreCheckCreate(MPI_Comm,PreCheck*); 91 PetscErrorCode PreCheckDestroy(PreCheck*); 92 PetscErrorCode PreCheckFunction(SNESLineSearch,Vec,Vec,PetscBool*,void*); 93 PetscErrorCode PreCheckSetFromOptions(PreCheck); 94 95 int main(int argc,char **argv) 96 { 97 SNES snes; /* nonlinear solver */ 98 Vec x,r,b; /* solution, residual, rhs vectors */ 99 AppCtx user; /* user-defined work context */ 100 PetscInt its; /* iterations for convergence */ 101 SNESConvergedReason reason; /* Check convergence */ 102 PetscBool alloc_star; /* Only allocate for the STAR stencil */ 103 PetscBool write_output; 104 char filename[PETSC_MAX_PATH_LEN] = "ex15.vts"; 105 PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.; 106 DM da; /* distributed array data structure */ 107 PreCheck precheck = NULL; /* precheck context for version in this file */ 108 PetscInt use_precheck; /* 0=none, 1=version in this file, 2=SNES-provided version */ 109 PetscReal precheck_angle; /* When manually setting the SNES-provided precheck function */ 110 PetscErrorCode ierr; 111 SNESLineSearch linesearch; 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Initialize program 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 117 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 118 119 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 120 Initialize problem parameters 121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 122 user.lambda = 0.0; user.p = 2.0; user.epsilon = 1e-5; user.source = 0.1; user.jtype = JAC_NEWTON;user.initial=-1; 123 user.blocks[0] = 1; user.blocks[1] = 1; user.kappa = 1e-3; 124 alloc_star = PETSC_FALSE; 125 use_precheck = 0; precheck_angle = 10.; 126 user.picard = PETSC_FALSE; 127 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"p-Bratu options",__FILE__);PetscCall(ierr); 128 { 129 PetscInt two=2; 130 PetscCall(PetscOptionsReal("-lambda","Bratu parameter","",user.lambda,&user.lambda,NULL)); 131 PetscCall(PetscOptionsReal("-p","Exponent `p' in p-Laplacian","",user.p,&user.p,NULL)); 132 PetscCall(PetscOptionsReal("-epsilon","Strain-regularization in p-Laplacian","",user.epsilon,&user.epsilon,NULL)); 133 PetscCall(PetscOptionsReal("-source","Constant source term","",user.source,&user.source,NULL)); 134 PetscCall(PetscOptionsEnum("-jtype","Jacobian approximation to assemble","",JacTypes,(PetscEnum)user.jtype,(PetscEnum*)&user.jtype,NULL)); 135 PetscCall(PetscOptionsName("-picard","Solve with defect-correction Picard iteration","",&user.picard)); 136 if (user.picard) { 137 user.jtype = JAC_PICARD; 138 PetscCheck(user.p == 3,PETSC_COMM_WORLD,PETSC_ERR_SUP,"Picard iteration is only supported for p == 3"); 139 /* the Picard linearization only requires a 5 point stencil, while the Newton linearization requires a 9 point stencil */ 140 /* hence allocating the 5 point stencil gives the same convergence as the 9 point stencil since the extra stencil points are not used */ 141 PetscCall(PetscOptionsBool("-alloc_star","Allocate for STAR stencil (5-point)","",alloc_star,&alloc_star,NULL)); 142 } 143 PetscCall(PetscOptionsInt("-precheck","Use a pre-check correction intended for use with Picard iteration 1=this version, 2=library","",use_precheck,&use_precheck,NULL)); 144 PetscCall(PetscOptionsInt("-initial","Initial conditions type (-1: default, 0: zero-valued, 1: peaked guess)","",user.initial,&user.initial,NULL)); 145 if (use_precheck == 2) { /* Using library version, get the angle */ 146 PetscCall(PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck_angle,&precheck_angle,NULL)); 147 } 148 PetscCall(PetscOptionsIntArray("-blocks","number of coefficient interfaces in x and y direction","",user.blocks,&two,NULL)); 149 PetscCall(PetscOptionsReal("-kappa","diffusivity in odd regions","",user.kappa,&user.kappa,NULL)); 150 PetscCall(PetscOptionsString("-o","Output solution in vts format","",filename,filename,sizeof(filename),&write_output)); 151 } 152 ierr = PetscOptionsEnd();PetscCall(ierr); 153 if (user.lambda > bratu_lambda_max || user.lambda < bratu_lambda_min) { 154 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"WARNING: lambda %g out of range for p=2\n",(double)user.lambda)); 155 } 156 157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 Create nonlinear solver context 159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 160 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 161 162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 163 Create distributed array (DMDA) to manage parallel grid and vectors 164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 165 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,alloc_star ? DMDA_STENCIL_STAR : DMDA_STENCIL_BOX,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da)); 166 PetscCall(DMSetFromOptions(da)); 167 PetscCall(DMSetUp(da)); 168 169 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 170 Extract global vectors from DM; then duplicate for remaining 171 vectors that are the same types 172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 173 PetscCall(DMCreateGlobalVector(da,&x)); 174 PetscCall(VecDuplicate(x,&r)); 175 PetscCall(VecDuplicate(x,&b)); 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 User can override with: 179 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 180 (unless user explicitly sets preconditioner) 181 -snes_mf_operator : form preconditioning matrix as set by the user, 182 but use matrix-free approx for Jacobian-vector 183 products within Newton-Krylov method 184 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 187 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 188 Set local function evaluation routine 189 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 190 PetscCall(DMSetApplicationContext(da, &user)); 191 PetscCall(SNESSetDM(snes,da)); 192 if (user.picard) { 193 /* 194 This is not really right requiring the user to call SNESSetFunction/Jacobian but the DMDASNESSetPicardLocal() cannot access 195 the SNES to set it 196 */ 197 ierr = DMDASNESSetPicardLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionPicardLocal, 198 (PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user);PetscCall(ierr); 199 PetscCall(SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user)); 200 PetscCall(SNESSetJacobian(snes,NULL,NULL,SNESPicardComputeJacobian,&user)); 201 } else { 202 PetscCall(DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user)); 203 PetscCall(DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user)); 204 } 205 206 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207 Customize nonlinear solver; set runtime options 208 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 209 PetscCall(SNESSetFromOptions(snes)); 210 PetscCall(SNESSetNGS(snes,NonlinearGS,&user)); 211 PetscCall(SNESGetLineSearch(snes, &linesearch)); 212 /* Set up the precheck context if requested */ 213 if (use_precheck == 1) { /* Use the precheck routines in this file */ 214 PetscCall(PreCheckCreate(PETSC_COMM_WORLD,&precheck)); 215 PetscCall(PreCheckSetFromOptions(precheck)); 216 PetscCall(SNESLineSearchSetPreCheck(linesearch,PreCheckFunction,precheck)); 217 } else if (use_precheck == 2) { /* Use the version provided by the library */ 218 PetscCall(SNESLineSearchSetPreCheck(linesearch,SNESLineSearchPreCheckPicard,&precheck_angle)); 219 } 220 221 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 222 Evaluate initial guess 223 Note: The user should initialize the vector, x, with the initial guess 224 for the nonlinear solver prior to calling SNESSolve(). In particular, 225 to employ an initial guess of zero, the user should explicitly set 226 this vector to zero by calling VecSet(). 227 */ 228 229 PetscCall(FormInitialGuess(&user,da,x)); 230 PetscCall(FormRHS(&user,da,b)); 231 232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233 Solve nonlinear system 234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 235 PetscCall(SNESSolve(snes,b,x)); 236 PetscCall(SNESGetIterationNumber(snes,&its)); 237 PetscCall(SNESGetConvergedReason(snes,&reason)); 238 239 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s Number of nonlinear iterations = %D\n",SNESConvergedReasons[reason],its)); 240 241 if (write_output) { 242 PetscViewer viewer; 243 PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer)); 244 PetscCall(VecView(x,viewer)); 245 PetscCall(PetscViewerDestroy(&viewer)); 246 } 247 248 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 249 Free work space. All PETSc objects should be destroyed when they 250 are no longer needed. 251 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 252 253 PetscCall(VecDestroy(&x)); 254 PetscCall(VecDestroy(&r)); 255 PetscCall(VecDestroy(&b)); 256 PetscCall(SNESDestroy(&snes)); 257 PetscCall(DMDestroy(&da)); 258 PetscCall(PreCheckDestroy(&precheck)); 259 PetscCall(PetscFinalize()); 260 return 0; 261 } 262 263 /* ------------------------------------------------------------------- */ 264 /* 265 FormInitialGuess - Forms initial approximation. 266 267 Input Parameters: 268 user - user-defined application context 269 X - vector 270 271 Output Parameter: 272 X - vector 273 */ 274 static PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X) 275 { 276 PetscInt i,j,Mx,My,xs,ys,xm,ym; 277 PetscReal temp1,temp,hx,hy; 278 PetscScalar **x; 279 280 PetscFunctionBeginUser; 281 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)); 282 283 hx = 1.0/(PetscReal)(Mx-1); 284 hy = 1.0/(PetscReal)(My-1); 285 temp1 = user->lambda / (user->lambda + 1.); 286 287 /* 288 Get a pointer to vector data. 289 - For default PETSc vectors, VecGetArray() returns a pointer to 290 the data array. Otherwise, the routine is implementation dependent. 291 - You MUST call VecRestoreArray() when you no longer need access to 292 the array. 293 */ 294 PetscCall(DMDAVecGetArray(da,X,&x)); 295 296 /* 297 Get local grid boundaries (for 2-dimensional DA): 298 xs, ys - starting grid indices (no ghost points) 299 xm, ym - widths of local grid (no ghost points) 300 301 */ 302 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 303 304 /* 305 Compute initial guess over the locally owned part of the grid 306 */ 307 for (j=ys; j<ys+ym; j++) { 308 temp = (PetscReal)(PetscMin(j,My-j-1))*hy; 309 for (i=xs; i<xs+xm; i++) { 310 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 311 /* boundary conditions are all zero Dirichlet */ 312 x[j][i] = 0.0; 313 } else { 314 if (user->initial == -1) { 315 if (user->lambda != 0) { 316 x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp)); 317 } else { 318 /* The solution above is an exact solution for lambda=0, this avoids "accidentally" starting 319 * with an exact solution. */ 320 const PetscReal 321 xx = 2*(PetscReal)i/(Mx-1) - 1, 322 yy = 2*(PetscReal)j/(My-1) - 1; 323 x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy; 324 } 325 } else if (user->initial == 0) { 326 x[j][i] = 0.; 327 } else if (user->initial == 1) { 328 const PetscReal 329 xx = 2*(PetscReal)i/(Mx-1) - 1, 330 yy = 2*(PetscReal)j/(My-1) - 1; 331 x[j][i] = (1 - xx*xx) * (1-yy*yy) * xx * yy; 332 } else { 333 if (user->lambda != 0) { 334 x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp)); 335 } else { 336 x[j][i] = 0.5*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp)); 337 } 338 } 339 } 340 } 341 } 342 /* 343 Restore vector 344 */ 345 PetscCall(DMDAVecRestoreArray(da,X,&x)); 346 PetscFunctionReturn(0); 347 } 348 349 /* 350 FormRHS - Forms constant RHS for the problem. 351 352 Input Parameters: 353 user - user-defined application context 354 B - RHS vector 355 356 Output Parameter: 357 B - vector 358 */ 359 static PetscErrorCode FormRHS(AppCtx *user,DM da,Vec B) 360 { 361 PetscInt i,j,Mx,My,xs,ys,xm,ym; 362 PetscReal hx,hy; 363 PetscScalar **b; 364 365 PetscFunctionBeginUser; 366 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)); 367 368 hx = 1.0/(PetscReal)(Mx-1); 369 hy = 1.0/(PetscReal)(My-1); 370 PetscCall(DMDAVecGetArray(da,B,&b)); 371 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 372 for (j=ys; j<ys+ym; j++) { 373 for (i=xs; i<xs+xm; i++) { 374 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 375 b[j][i] = 0.0; 376 } else { 377 b[j][i] = hx*hy*user->source; 378 } 379 } 380 } 381 PetscCall(DMDAVecRestoreArray(da,B,&b)); 382 PetscFunctionReturn(0); 383 } 384 385 static inline PetscReal kappa(const AppCtx *ctx,PetscReal x,PetscReal y) 386 { 387 return (((PetscInt)(x*ctx->blocks[0])) + ((PetscInt)(y*ctx->blocks[1]))) % 2 ? ctx->kappa : 1.0; 388 } 389 /* p-Laplacian diffusivity */ 390 static inline PetscScalar eta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy) 391 { 392 return kappa(ctx,x,y) * PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-2.)); 393 } 394 static inline PetscScalar deta(const AppCtx *ctx,PetscReal x,PetscReal y,PetscScalar ux,PetscScalar uy) 395 { 396 return (ctx->p == 2) 397 ? 0 398 : kappa(ctx,x,y)*PetscPowScalar(PetscSqr(ctx->epsilon)+0.5*(ux*ux + uy*uy),0.5*(ctx->p-4)) * 0.5 * (ctx->p-2.); 399 } 400 401 /* ------------------------------------------------------------------- */ 402 /* 403 FormFunctionLocal - Evaluates nonlinear function, F(x). 404 */ 405 static PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user) 406 { 407 PetscReal hx,hy,dhx,dhy,sc; 408 PetscInt i,j; 409 PetscScalar eu; 410 411 PetscFunctionBeginUser; 412 hx = 1.0/(PetscReal)(info->mx-1); 413 hy = 1.0/(PetscReal)(info->my-1); 414 sc = hx*hy*user->lambda; 415 dhx = 1/hx; 416 dhy = 1/hy; 417 /* 418 Compute function over the locally owned part of the grid 419 */ 420 for (j=info->ys; j<info->ys+info->ym; j++) { 421 for (i=info->xs; i<info->xs+info->xm; i++) { 422 PetscReal xx = i*hx,yy = j*hy; 423 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { 424 f[j][i] = x[j][i]; 425 } else { 426 const PetscScalar 427 u = x[j][i], 428 ux_E = dhx*(x[j][i+1]-x[j][i]), 429 uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]), 430 ux_W = dhx*(x[j][i]-x[j][i-1]), 431 uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]), 432 ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]), 433 uy_N = dhy*(x[j+1][i]-x[j][i]), 434 ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]), 435 uy_S = dhy*(x[j][i]-x[j-1][i]), 436 e_E = eta(user,xx,yy,ux_E,uy_E), 437 e_W = eta(user,xx,yy,ux_W,uy_W), 438 e_N = eta(user,xx,yy,ux_N,uy_N), 439 e_S = eta(user,xx,yy,ux_S,uy_S), 440 uxx = -hy * (e_E*ux_E - e_W*ux_W), 441 uyy = -hx * (e_N*uy_N - e_S*uy_S); 442 if (sc) eu = PetscExpScalar(u); 443 else eu = 0.; 444 /* For p=2, these terms decay to: 445 uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx 446 uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy 447 */ 448 f[j][i] = uxx + uyy - sc*eu; 449 } 450 } 451 } 452 PetscCall(PetscLogFlops(info->xm*info->ym*(72.0))); 453 PetscFunctionReturn(0); 454 } 455 456 /* 457 This is the opposite sign of the part of FormFunctionLocal that excludes the A(x) x part of the operation, 458 that is FormFunction applies A(x) x - b(x) while this applies b(x) because for Picard we think of it as solving A(x) x = b(x) 459 460 */ 461 static PetscErrorCode FormFunctionPicardLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user) 462 { 463 PetscReal hx,hy,sc; 464 PetscInt i,j; 465 466 PetscFunctionBeginUser; 467 hx = 1.0/(PetscReal)(info->mx-1); 468 hy = 1.0/(PetscReal)(info->my-1); 469 sc = hx*hy*user->lambda; 470 /* 471 Compute function over the locally owned part of the grid 472 */ 473 for (j=info->ys; j<info->ys+info->ym; j++) { 474 for (i=info->xs; i<info->xs+info->xm; i++) { 475 if (!(i == 0 || j == 0 || i == info->mx-1 || j == info->my-1)) { 476 const PetscScalar u = x[j][i]; 477 f[j][i] = sc*PetscExpScalar(u); 478 } else { 479 f[j][i] = 0.0; /* this is zero because the A(x) x term forces the x to be zero on the boundary */ 480 } 481 } 482 } 483 PetscCall(PetscLogFlops(info->xm*info->ym*2.0)); 484 PetscFunctionReturn(0); 485 } 486 487 /* 488 FormJacobianLocal - Evaluates Jacobian matrix. 489 */ 490 static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat J,Mat B,AppCtx *user) 491 { 492 PetscInt i,j; 493 MatStencil col[9],row; 494 PetscScalar v[9]; 495 PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc; 496 497 PetscFunctionBeginUser; 498 PetscCall(MatZeroEntries(B)); 499 hx = 1.0/(PetscReal)(info->mx-1); 500 hy = 1.0/(PetscReal)(info->my-1); 501 sc = hx*hy*user->lambda; 502 dhx = 1/hx; 503 dhy = 1/hy; 504 hxdhy = hx/hy; 505 hydhx = hy/hx; 506 507 /* 508 Compute entries for the locally owned part of the Jacobian. 509 - PETSc parallel matrix formats are partitioned by 510 contiguous chunks of rows across the processors. 511 - Each processor needs to insert only elements that it owns 512 locally (but any non-local elements will be sent to the 513 appropriate processor during matrix assembly). 514 - Here, we set all entries for a particular row at once. 515 */ 516 for (j=info->ys; j<info->ys+info->ym; j++) { 517 for (i=info->xs; i<info->xs+info->xm; i++) { 518 PetscReal xx = i*hx,yy = j*hy; 519 row.j = j; row.i = i; 520 /* boundary points */ 521 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { 522 v[0] = 1.0; 523 PetscCall(MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES)); 524 } else { 525 /* interior grid points */ 526 const PetscScalar 527 ux_E = dhx*(x[j][i+1]-x[j][i]), 528 uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]), 529 ux_W = dhx*(x[j][i]-x[j][i-1]), 530 uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]), 531 ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]), 532 uy_N = dhy*(x[j+1][i]-x[j][i]), 533 ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]), 534 uy_S = dhy*(x[j][i]-x[j-1][i]), 535 u = x[j][i], 536 e_E = eta(user,xx,yy,ux_E,uy_E), 537 e_W = eta(user,xx,yy,ux_W,uy_W), 538 e_N = eta(user,xx,yy,ux_N,uy_N), 539 e_S = eta(user,xx,yy,ux_S,uy_S), 540 de_E = deta(user,xx,yy,ux_E,uy_E), 541 de_W = deta(user,xx,yy,ux_W,uy_W), 542 de_N = deta(user,xx,yy,ux_N,uy_N), 543 de_S = deta(user,xx,yy,ux_S,uy_S), 544 skew_E = de_E*ux_E*uy_E, 545 skew_W = de_W*ux_W*uy_W, 546 skew_N = de_N*ux_N*uy_N, 547 skew_S = de_S*ux_S*uy_S, 548 cross_EW = 0.25*(skew_E - skew_W), 549 cross_NS = 0.25*(skew_N - skew_S), 550 newt_E = e_E+de_E*PetscSqr(ux_E), 551 newt_W = e_W+de_W*PetscSqr(ux_W), 552 newt_N = e_N+de_N*PetscSqr(uy_N), 553 newt_S = e_S+de_S*PetscSqr(uy_S); 554 /* interior grid points */ 555 switch (user->jtype) { 556 case JAC_BRATU: 557 /* Jacobian from p=2 */ 558 v[0] = -hxdhy; col[0].j = j-1; col[0].i = i; 559 v[1] = -hydhx; col[1].j = j; col[1].i = i-1; 560 v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u); col[2].j = row.j; col[2].i = row.i; 561 v[3] = -hydhx; col[3].j = j; col[3].i = i+1; 562 v[4] = -hxdhy; col[4].j = j+1; col[4].i = i; 563 PetscCall(MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES)); 564 break; 565 case JAC_PICARD: 566 /* Jacobian arising from Picard linearization */ 567 v[0] = -hxdhy*e_S; col[0].j = j-1; col[0].i = i; 568 v[1] = -hydhx*e_W; col[1].j = j; col[1].i = i-1; 569 v[2] = (e_W+e_E)*hydhx + (e_S+e_N)*hxdhy; col[2].j = row.j; col[2].i = row.i; 570 v[3] = -hydhx*e_E; col[3].j = j; col[3].i = i+1; 571 v[4] = -hxdhy*e_N; col[4].j = j+1; col[4].i = i; 572 PetscCall(MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES)); 573 break; 574 case JAC_STAR: 575 /* Full Jacobian, but only a star stencil */ 576 col[0].j = j-1; col[0].i = i; 577 col[1].j = j; col[1].i = i-1; 578 col[2].j = j; col[2].i = i; 579 col[3].j = j; col[3].i = i+1; 580 col[4].j = j+1; col[4].i = i; 581 v[0] = -hxdhy*newt_S + cross_EW; 582 v[1] = -hydhx*newt_W + cross_NS; 583 v[2] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u); 584 v[3] = -hydhx*newt_E - cross_NS; 585 v[4] = -hxdhy*newt_N - cross_EW; 586 PetscCall(MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES)); 587 break; 588 case JAC_NEWTON: 589 /* The Jacobian is 590 591 -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u 592 593 */ 594 col[0].j = j-1; col[0].i = i-1; 595 col[1].j = j-1; col[1].i = i; 596 col[2].j = j-1; col[2].i = i+1; 597 col[3].j = j; col[3].i = i-1; 598 col[4].j = j; col[4].i = i; 599 col[5].j = j; col[5].i = i+1; 600 col[6].j = j+1; col[6].i = i-1; 601 col[7].j = j+1; col[7].i = i; 602 col[8].j = j+1; col[8].i = i+1; 603 v[0] = -0.25*(skew_S + skew_W); 604 v[1] = -hxdhy*newt_S + cross_EW; 605 v[2] = 0.25*(skew_S + skew_E); 606 v[3] = -hydhx*newt_W + cross_NS; 607 v[4] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u); 608 v[5] = -hydhx*newt_E - cross_NS; 609 v[6] = 0.25*(skew_N + skew_W); 610 v[7] = -hxdhy*newt_N - cross_EW; 611 v[8] = -0.25*(skew_N + skew_E); 612 PetscCall(MatSetValuesStencil(B,1,&row,9,col,v,INSERT_VALUES)); 613 break; 614 default: 615 SETERRQ(PetscObjectComm((PetscObject)info->da),PETSC_ERR_SUP,"Jacobian type %d not implemented",user->jtype); 616 } 617 } 618 } 619 } 620 621 /* 622 Assemble matrix, using the 2-step process: 623 MatAssemblyBegin(), MatAssemblyEnd(). 624 */ 625 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 626 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 627 628 if (J != B) { 629 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 630 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 631 } 632 633 /* 634 Tell the matrix we will never add a new nonzero location to the 635 matrix. If we do, it will generate an error. 636 */ 637 if (user->jtype == JAC_NEWTON) { 638 PetscCall(PetscLogFlops(info->xm*info->ym*(131.0))); 639 } 640 PetscFunctionReturn(0); 641 } 642 643 /*********************************************************** 644 * PreCheck implementation 645 ***********************************************************/ 646 PetscErrorCode PreCheckSetFromOptions(PreCheck precheck) 647 { 648 PetscErrorCode ierr; 649 PetscBool flg; 650 651 PetscFunctionBeginUser; 652 ierr = PetscOptionsBegin(precheck->comm,NULL,"PreCheck Options","none");PetscCall(ierr); 653 PetscCall(PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,NULL)); 654 flg = PETSC_FALSE; 655 PetscCall(PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,NULL)); 656 if (flg) { 657 PetscCall(PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor)); 658 } 659 ierr = PetscOptionsEnd();PetscCall(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 /* 664 Compare the direction of the current and previous step, modify the current step accordingly 665 */ 666 PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx) 667 { 668 PreCheck precheck; 669 Vec Ylast; 670 PetscScalar dot; 671 PetscInt iter; 672 PetscReal ynorm,ylastnorm,theta,angle_radians; 673 SNES snes; 674 675 PetscFunctionBeginUser; 676 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 677 precheck = (PreCheck)ctx; 678 if (!precheck->Ylast) PetscCall(VecDuplicate(Y,&precheck->Ylast)); 679 Ylast = precheck->Ylast; 680 PetscCall(SNESGetIterationNumber(snes,&iter)); 681 if (iter < 1) { 682 PetscCall(VecCopy(Y,Ylast)); 683 *changed = PETSC_FALSE; 684 PetscFunctionReturn(0); 685 } 686 687 PetscCall(VecDot(Y,Ylast,&dot)); 688 PetscCall(VecNorm(Y,NORM_2,&ynorm)); 689 PetscCall(VecNorm(Ylast,NORM_2,&ylastnorm)); 690 /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */ 691 theta = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0)); 692 angle_radians = precheck->angle * PETSC_PI / 180.; 693 if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) { 694 /* Modify the step Y */ 695 PetscReal alpha,ydiffnorm; 696 PetscCall(VecAXPY(Ylast,-1.0,Y)); 697 PetscCall(VecNorm(Ylast,NORM_2,&ydiffnorm)); 698 alpha = ylastnorm / ydiffnorm; 699 PetscCall(VecCopy(Y,Ylast)); 700 PetscCall(VecScale(Y,alpha)); 701 if (precheck->monitor) { 702 PetscCall(PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees less than threshold %g, corrected step by alpha=%g\n",(double)(theta*180./PETSC_PI),(double)precheck->angle,(double)alpha)); 703 } 704 } else { 705 PetscCall(VecCopy(Y,Ylast)); 706 *changed = PETSC_FALSE; 707 if (precheck->monitor) { 708 PetscCall(PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees exceeds threshold %g, no correction applied\n",(double)(theta*180./PETSC_PI),(double)precheck->angle)); 709 } 710 } 711 PetscFunctionReturn(0); 712 } 713 714 PetscErrorCode PreCheckDestroy(PreCheck *precheck) 715 { 716 PetscFunctionBeginUser; 717 if (!*precheck) PetscFunctionReturn(0); 718 PetscCall(VecDestroy(&(*precheck)->Ylast)); 719 PetscCall(PetscViewerDestroy(&(*precheck)->monitor)); 720 PetscCall(PetscFree(*precheck)); 721 PetscFunctionReturn(0); 722 } 723 724 PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck) 725 { 726 PetscFunctionBeginUser; 727 PetscCall(PetscNew(precheck)); 728 729 (*precheck)->comm = comm; 730 (*precheck)->angle = 10.; /* only active if angle is less than 10 degrees */ 731 PetscFunctionReturn(0); 732 } 733 734 /* 735 Applies some sweeps on nonlinear Gauss-Seidel on each process 736 737 */ 738 PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx) 739 { 740 PetscInt i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m; 741 PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc; 742 PetscScalar **x,**b,bij,F,F0=0,J,y,u,eu; 743 PetscReal atol,rtol,stol; 744 DM da; 745 AppCtx *user = (AppCtx*)ctx; 746 Vec localX,localB; 747 DMDALocalInfo info; 748 749 PetscFunctionBeginUser; 750 PetscCall(SNESGetDM(snes,&da)); 751 PetscCall(DMDAGetLocalInfo(da,&info)); 752 753 hx = 1.0/(PetscReal)(info.mx-1); 754 hy = 1.0/(PetscReal)(info.my-1); 755 sc = hx*hy*user->lambda; 756 dhx = 1/hx; 757 dhy = 1/hy; 758 hxdhy = hx/hy; 759 hydhx = hy/hx; 760 761 tot_its = 0; 762 PetscCall(SNESNGSGetSweeps(snes,&sweeps)); 763 PetscCall(SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its)); 764 PetscCall(DMGetLocalVector(da,&localX)); 765 if (B) { 766 PetscCall(DMGetLocalVector(da,&localB)); 767 } 768 if (B) { 769 PetscCall(DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB)); 770 PetscCall(DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB)); 771 } 772 if (B) PetscCall(DMDAVecGetArrayRead(da,localB,&b)); 773 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 774 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 775 PetscCall(DMDAVecGetArray(da,localX,&x)); 776 for (l=0; l<sweeps; l++) { 777 /* 778 Get local grid boundaries (for 2-dimensional DMDA): 779 xs, ys - starting grid indices (no ghost points) 780 xm, ym - widths of local grid (no ghost points) 781 */ 782 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 783 for (m=0; m<2; m++) { 784 for (j=ys; j<ys+ym; j++) { 785 for (i=xs+(m+j)%2; i<xs+xm; i+=2) { 786 PetscReal xx = i*hx,yy = j*hy; 787 if (B) bij = b[j][i]; 788 else bij = 0.; 789 790 if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) { 791 /* boundary conditions are all zero Dirichlet */ 792 x[j][i] = 0.0 + bij; 793 } else { 794 const PetscScalar 795 u_E = x[j][i+1], 796 u_W = x[j][i-1], 797 u_N = x[j+1][i], 798 u_S = x[j-1][i]; 799 const PetscScalar 800 uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]), 801 uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]), 802 ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]), 803 ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]); 804 u = x[j][i]; 805 for (k=0; k<its; k++) { 806 const PetscScalar 807 ux_E = dhx*(u_E-u), 808 ux_W = dhx*(u-u_W), 809 uy_N = dhy*(u_N-u), 810 uy_S = dhy*(u-u_S), 811 e_E = eta(user,xx,yy,ux_E,uy_E), 812 e_W = eta(user,xx,yy,ux_W,uy_W), 813 e_N = eta(user,xx,yy,ux_N,uy_N), 814 e_S = eta(user,xx,yy,ux_S,uy_S), 815 de_E = deta(user,xx,yy,ux_E,uy_E), 816 de_W = deta(user,xx,yy,ux_W,uy_W), 817 de_N = deta(user,xx,yy,ux_N,uy_N), 818 de_S = deta(user,xx,yy,ux_S,uy_S), 819 newt_E = e_E+de_E*PetscSqr(ux_E), 820 newt_W = e_W+de_W*PetscSqr(ux_W), 821 newt_N = e_N+de_N*PetscSqr(uy_N), 822 newt_S = e_S+de_S*PetscSqr(uy_S), 823 uxx = -hy * (e_E*ux_E - e_W*ux_W), 824 uyy = -hx * (e_N*uy_N - e_S*uy_S); 825 826 if (sc) eu = PetscExpScalar(u); 827 else eu = 0; 828 829 F = uxx + uyy - sc*eu - bij; 830 if (k == 0) F0 = F; 831 J = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*eu; 832 y = F/J; 833 u -= y; 834 tot_its++; 835 if (atol > PetscAbsReal(PetscRealPart(F)) || 836 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) || 837 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) { 838 break; 839 } 840 } 841 x[j][i] = u; 842 } 843 } 844 } 845 } 846 /* 847 x Restore vector 848 */ 849 } 850 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 851 PetscCall(DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X)); 852 PetscCall(DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X)); 853 PetscCall(PetscLogFlops(tot_its*(118.0))); 854 PetscCall(DMRestoreLocalVector(da,&localX)); 855 if (B) { 856 PetscCall(DMDAVecRestoreArrayRead(da,localB,&b)); 857 PetscCall(DMRestoreLocalVector(da,&localB)); 858 } 859 PetscFunctionReturn(0); 860 } 861 862 /*TEST 863 864 test: 865 nsize: 2 866 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON 867 requires: !single 868 869 test: 870 suffix: 2 871 nsize: 2 872 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1 873 requires: !single 874 875 test: 876 suffix: 3 877 nsize: 2 878 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3 879 requires: !single 880 881 test: 882 suffix: 3_star 883 nsize: 2 884 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3 -alloc_star 885 output_file: output/ex15_3.out 886 requires: !single 887 888 test: 889 suffix: 4 890 args: -snes_monitor_short -snes_type newtonls -npc_snes_type ngs -snes_npc_side left -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -ksp_monitor_short -pc_type none 891 requires: !single 892 893 test: 894 suffix: lag_jac 895 nsize: 4 896 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 6.0 -lambda 0 -jtype NEWTON -snes_type ngmres -npc_snes_type newtonls -npc_snes_lag_jacobian 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_jacobian_persists 897 requires: !single 898 899 test: 900 suffix: lag_pc 901 nsize: 4 902 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 6.0 -lambda 0 -jtype NEWTON -snes_type ngmres -npc_snes_type newtonls -npc_snes_lag_preconditioner 5 -npc_pc_type asm -npc_ksp_converged_reason -npc_snes_lag_preconditioner_persists 903 requires: !single 904 905 test: 906 suffix: nleqerr 907 args: -snes_monitor_short -snes_type newtonls -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -snes_linesearch_monitor -pc_type lu -snes_linesearch_type nleqerr 908 requires: !single 909 910 test: 911 suffix: mf 912 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator 913 requires: !single 914 915 test: 916 suffix: mf_picard 917 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator -picard 918 requires: !single 919 output_file: output/ex15_mf.out 920 921 test: 922 suffix: fd_picard 923 args: -snes_monitor_short -pc_type lu -da_refine 2 -p 3 -ksp_rtol 1.e-12 -snes_fd -picard 924 requires: !single 925 926 test: 927 suffix: fd_color_picard 928 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_fd_color -picard 929 requires: !single 930 output_file: output/ex15_mf.out 931 932 TEST*/ 933