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 SNESLineSearch linesearch; 111 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Initialize program 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 116 PetscFunctionBeginUser; 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 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"p-Bratu options",__FILE__); 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 PetscOptionsEnd(); 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 PetscCall(DMDASNESSetPicardLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionPicardLocal, 198 (PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,&user)); 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 = %" PetscInt_FMT "\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 PetscBool flg; 649 650 PetscFunctionBeginUser; 651 PetscOptionsBegin(precheck->comm,NULL,"PreCheck Options","none"); 652 PetscCall(PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,NULL)); 653 flg = PETSC_FALSE; 654 PetscCall(PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,NULL)); 655 if (flg) { 656 PetscCall(PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor)); 657 } 658 PetscOptionsEnd(); 659 PetscFunctionReturn(0); 660 } 661 662 /* 663 Compare the direction of the current and previous step, modify the current step accordingly 664 */ 665 PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx) 666 { 667 PreCheck precheck; 668 Vec Ylast; 669 PetscScalar dot; 670 PetscInt iter; 671 PetscReal ynorm,ylastnorm,theta,angle_radians; 672 SNES snes; 673 674 PetscFunctionBeginUser; 675 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 676 precheck = (PreCheck)ctx; 677 if (!precheck->Ylast) PetscCall(VecDuplicate(Y,&precheck->Ylast)); 678 Ylast = precheck->Ylast; 679 PetscCall(SNESGetIterationNumber(snes,&iter)); 680 if (iter < 1) { 681 PetscCall(VecCopy(Y,Ylast)); 682 *changed = PETSC_FALSE; 683 PetscFunctionReturn(0); 684 } 685 686 PetscCall(VecDot(Y,Ylast,&dot)); 687 PetscCall(VecNorm(Y,NORM_2,&ynorm)); 688 PetscCall(VecNorm(Ylast,NORM_2,&ylastnorm)); 689 /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */ 690 theta = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0)); 691 angle_radians = precheck->angle * PETSC_PI / 180.; 692 if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) { 693 /* Modify the step Y */ 694 PetscReal alpha,ydiffnorm; 695 PetscCall(VecAXPY(Ylast,-1.0,Y)); 696 PetscCall(VecNorm(Ylast,NORM_2,&ydiffnorm)); 697 alpha = ylastnorm / ydiffnorm; 698 PetscCall(VecCopy(Y,Ylast)); 699 PetscCall(VecScale(Y,alpha)); 700 if (precheck->monitor) { 701 PetscCall(PetscViewerASCIIPrintf(precheck->monitor,"Angle %g degrees less than threshold %g, corrected step by alpha=%g\n",(double)(theta*180./PETSC_PI),(double)precheck->angle,(double)alpha)); 702 } 703 } else { 704 PetscCall(VecCopy(Y,Ylast)); 705 *changed = PETSC_FALSE; 706 if (precheck->monitor) { 707 PetscCall(PetscViewerASCIIPrintf(precheck->monitor,"Angle %g degrees exceeds threshold %g, no correction applied\n",(double)(theta*180./PETSC_PI),(double)precheck->angle)); 708 } 709 } 710 PetscFunctionReturn(0); 711 } 712 713 PetscErrorCode PreCheckDestroy(PreCheck *precheck) 714 { 715 PetscFunctionBeginUser; 716 if (!*precheck) PetscFunctionReturn(0); 717 PetscCall(VecDestroy(&(*precheck)->Ylast)); 718 PetscCall(PetscViewerDestroy(&(*precheck)->monitor)); 719 PetscCall(PetscFree(*precheck)); 720 PetscFunctionReturn(0); 721 } 722 723 PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck) 724 { 725 PetscFunctionBeginUser; 726 PetscCall(PetscNew(precheck)); 727 728 (*precheck)->comm = comm; 729 (*precheck)->angle = 10.; /* only active if angle is less than 10 degrees */ 730 PetscFunctionReturn(0); 731 } 732 733 /* 734 Applies some sweeps on nonlinear Gauss-Seidel on each process 735 736 */ 737 PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx) 738 { 739 PetscInt i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m; 740 PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc; 741 PetscScalar **x,**b,bij,F,F0=0,J,y,u,eu; 742 PetscReal atol,rtol,stol; 743 DM da; 744 AppCtx *user = (AppCtx*)ctx; 745 Vec localX,localB; 746 DMDALocalInfo info; 747 748 PetscFunctionBeginUser; 749 PetscCall(SNESGetDM(snes,&da)); 750 PetscCall(DMDAGetLocalInfo(da,&info)); 751 752 hx = 1.0/(PetscReal)(info.mx-1); 753 hy = 1.0/(PetscReal)(info.my-1); 754 sc = hx*hy*user->lambda; 755 dhx = 1/hx; 756 dhy = 1/hy; 757 hxdhy = hx/hy; 758 hydhx = hy/hx; 759 760 tot_its = 0; 761 PetscCall(SNESNGSGetSweeps(snes,&sweeps)); 762 PetscCall(SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its)); 763 PetscCall(DMGetLocalVector(da,&localX)); 764 if (B) { 765 PetscCall(DMGetLocalVector(da,&localB)); 766 } 767 if (B) { 768 PetscCall(DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB)); 769 PetscCall(DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB)); 770 } 771 if (B) PetscCall(DMDAVecGetArrayRead(da,localB,&b)); 772 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 773 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 774 PetscCall(DMDAVecGetArray(da,localX,&x)); 775 for (l=0; l<sweeps; l++) { 776 /* 777 Get local grid boundaries (for 2-dimensional DMDA): 778 xs, ys - starting grid indices (no ghost points) 779 xm, ym - widths of local grid (no ghost points) 780 */ 781 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 782 for (m=0; m<2; m++) { 783 for (j=ys; j<ys+ym; j++) { 784 for (i=xs+(m+j)%2; i<xs+xm; i+=2) { 785 PetscReal xx = i*hx,yy = j*hy; 786 if (B) bij = b[j][i]; 787 else bij = 0.; 788 789 if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) { 790 /* boundary conditions are all zero Dirichlet */ 791 x[j][i] = 0.0 + bij; 792 } else { 793 const PetscScalar 794 u_E = x[j][i+1], 795 u_W = x[j][i-1], 796 u_N = x[j+1][i], 797 u_S = x[j-1][i]; 798 const PetscScalar 799 uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]), 800 uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]), 801 ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]), 802 ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]); 803 u = x[j][i]; 804 for (k=0; k<its; k++) { 805 const PetscScalar 806 ux_E = dhx*(u_E-u), 807 ux_W = dhx*(u-u_W), 808 uy_N = dhy*(u_N-u), 809 uy_S = dhy*(u-u_S), 810 e_E = eta(user,xx,yy,ux_E,uy_E), 811 e_W = eta(user,xx,yy,ux_W,uy_W), 812 e_N = eta(user,xx,yy,ux_N,uy_N), 813 e_S = eta(user,xx,yy,ux_S,uy_S), 814 de_E = deta(user,xx,yy,ux_E,uy_E), 815 de_W = deta(user,xx,yy,ux_W,uy_W), 816 de_N = deta(user,xx,yy,ux_N,uy_N), 817 de_S = deta(user,xx,yy,ux_S,uy_S), 818 newt_E = e_E+de_E*PetscSqr(ux_E), 819 newt_W = e_W+de_W*PetscSqr(ux_W), 820 newt_N = e_N+de_N*PetscSqr(uy_N), 821 newt_S = e_S+de_S*PetscSqr(uy_S), 822 uxx = -hy * (e_E*ux_E - e_W*ux_W), 823 uyy = -hx * (e_N*uy_N - e_S*uy_S); 824 825 if (sc) eu = PetscExpScalar(u); 826 else eu = 0; 827 828 F = uxx + uyy - sc*eu - bij; 829 if (k == 0) F0 = F; 830 J = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*eu; 831 y = F/J; 832 u -= y; 833 tot_its++; 834 if (atol > PetscAbsReal(PetscRealPart(F)) || 835 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) || 836 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) { 837 break; 838 } 839 } 840 x[j][i] = u; 841 } 842 } 843 } 844 } 845 /* 846 x Restore vector 847 */ 848 } 849 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 850 PetscCall(DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X)); 851 PetscCall(DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X)); 852 PetscCall(PetscLogFlops(tot_its*(118.0))); 853 PetscCall(DMRestoreLocalVector(da,&localX)); 854 if (B) { 855 PetscCall(DMDAVecRestoreArrayRead(da,localB,&b)); 856 PetscCall(DMRestoreLocalVector(da,&localB)); 857 } 858 PetscFunctionReturn(0); 859 } 860 861 /*TEST 862 863 test: 864 nsize: 2 865 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON 866 requires: !single 867 868 test: 869 suffix: 2 870 nsize: 2 871 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1 872 requires: !single 873 874 test: 875 suffix: 3 876 nsize: 2 877 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3 878 requires: !single 879 880 test: 881 suffix: 3_star 882 nsize: 2 883 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 884 output_file: output/ex15_3.out 885 requires: !single 886 887 test: 888 suffix: 4 889 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 890 requires: !single 891 892 test: 893 suffix: lag_jac 894 nsize: 4 895 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 896 requires: !single 897 898 test: 899 suffix: lag_pc 900 nsize: 4 901 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 902 requires: !single 903 904 test: 905 suffix: nleqerr 906 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 907 requires: !single 908 909 test: 910 suffix: mf 911 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator 912 requires: !single 913 914 test: 915 suffix: mf_picard 916 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator -picard 917 requires: !single 918 output_file: output/ex15_mf.out 919 920 test: 921 suffix: fd_picard 922 args: -snes_monitor_short -pc_type lu -da_refine 2 -p 3 -ksp_rtol 1.e-12 -snes_fd -picard 923 requires: !single 924 925 test: 926 suffix: fd_color_picard 927 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_fd_color -picard 928 requires: !single 929 output_file: output/ex15_mf.out 930 931 TEST*/ 932