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