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 } else { 487 f[j][i] = 0.0; /* this is zero because the A(x) x term forces the x to be zero on the boundary */ 488 } 489 } 490 } 491 ierr = PetscLogFlops(info->xm*info->ym*2.0);CHKERRQ(ierr); 492 PetscFunctionReturn(0); 493 } 494 495 /* 496 FormJacobianLocal - Evaluates Jacobian matrix. 497 */ 498 static PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat J,Mat B,AppCtx *user) 499 { 500 PetscErrorCode ierr; 501 PetscInt i,j; 502 MatStencil col[9],row; 503 PetscScalar v[9]; 504 PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc; 505 506 PetscFunctionBeginUser; 507 ierr = MatZeroEntries(B);CHKERRQ(ierr); 508 hx = 1.0/(PetscReal)(info->mx-1); 509 hy = 1.0/(PetscReal)(info->my-1); 510 sc = hx*hy*user->lambda; 511 dhx = 1/hx; 512 dhy = 1/hy; 513 hxdhy = hx/hy; 514 hydhx = hy/hx; 515 516 /* 517 Compute entries for the locally owned part of the Jacobian. 518 - PETSc parallel matrix formats are partitioned by 519 contiguous chunks of rows across the processors. 520 - Each processor needs to insert only elements that it owns 521 locally (but any non-local elements will be sent to the 522 appropriate processor during matrix assembly). 523 - Here, we set all entries for a particular row at once. 524 */ 525 for (j=info->ys; j<info->ys+info->ym; j++) { 526 for (i=info->xs; i<info->xs+info->xm; i++) { 527 PetscReal xx = i*hx,yy = j*hy; 528 row.j = j; row.i = i; 529 /* boundary points */ 530 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { 531 v[0] = 1.0; 532 ierr = MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr); 533 } else { 534 /* interior grid points */ 535 const PetscScalar 536 ux_E = dhx*(x[j][i+1]-x[j][i]), 537 uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]), 538 ux_W = dhx*(x[j][i]-x[j][i-1]), 539 uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]), 540 ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]), 541 uy_N = dhy*(x[j+1][i]-x[j][i]), 542 ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]), 543 uy_S = dhy*(x[j][i]-x[j-1][i]), 544 u = x[j][i], 545 e_E = eta(user,xx,yy,ux_E,uy_E), 546 e_W = eta(user,xx,yy,ux_W,uy_W), 547 e_N = eta(user,xx,yy,ux_N,uy_N), 548 e_S = eta(user,xx,yy,ux_S,uy_S), 549 de_E = deta(user,xx,yy,ux_E,uy_E), 550 de_W = deta(user,xx,yy,ux_W,uy_W), 551 de_N = deta(user,xx,yy,ux_N,uy_N), 552 de_S = deta(user,xx,yy,ux_S,uy_S), 553 skew_E = de_E*ux_E*uy_E, 554 skew_W = de_W*ux_W*uy_W, 555 skew_N = de_N*ux_N*uy_N, 556 skew_S = de_S*ux_S*uy_S, 557 cross_EW = 0.25*(skew_E - skew_W), 558 cross_NS = 0.25*(skew_N - skew_S), 559 newt_E = e_E+de_E*PetscSqr(ux_E), 560 newt_W = e_W+de_W*PetscSqr(ux_W), 561 newt_N = e_N+de_N*PetscSqr(uy_N), 562 newt_S = e_S+de_S*PetscSqr(uy_S); 563 /* interior grid points */ 564 switch (user->jtype) { 565 case JAC_BRATU: 566 /* Jacobian from p=2 */ 567 v[0] = -hxdhy; col[0].j = j-1; col[0].i = i; 568 v[1] = -hydhx; col[1].j = j; col[1].i = i-1; 569 v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(u); col[2].j = row.j; col[2].i = row.i; 570 v[3] = -hydhx; col[3].j = j; col[3].i = i+1; 571 v[4] = -hxdhy; col[4].j = j+1; col[4].i = i; 572 ierr = MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 573 break; 574 case JAC_PICARD: 575 /* Jacobian arising from Picard linearization */ 576 v[0] = -hxdhy*e_S; col[0].j = j-1; col[0].i = i; 577 v[1] = -hydhx*e_W; col[1].j = j; col[1].i = i-1; 578 v[2] = (e_W+e_E)*hydhx + (e_S+e_N)*hxdhy; col[2].j = row.j; col[2].i = row.i; 579 v[3] = -hydhx*e_E; col[3].j = j; col[3].i = i+1; 580 v[4] = -hxdhy*e_N; col[4].j = j+1; col[4].i = i; 581 ierr = MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 582 break; 583 case JAC_STAR: 584 /* Full Jacobian, but only a star stencil */ 585 col[0].j = j-1; col[0].i = i; 586 col[1].j = j; col[1].i = i-1; 587 col[2].j = j; col[2].i = i; 588 col[3].j = j; col[3].i = i+1; 589 col[4].j = j+1; col[4].i = i; 590 v[0] = -hxdhy*newt_S + cross_EW; 591 v[1] = -hydhx*newt_W + cross_NS; 592 v[2] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u); 593 v[3] = -hydhx*newt_E - cross_NS; 594 v[4] = -hxdhy*newt_N - cross_EW; 595 ierr = MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 596 break; 597 case JAC_NEWTON: 598 /** The Jacobian is 599 * 600 * -div [ eta (grad u) + deta (grad u0 . grad u) grad u0 ] - (eE u0) u 601 * 602 **/ 603 col[0].j = j-1; col[0].i = i-1; 604 col[1].j = j-1; col[1].i = i; 605 col[2].j = j-1; col[2].i = i+1; 606 col[3].j = j; col[3].i = i-1; 607 col[4].j = j; col[4].i = i; 608 col[5].j = j; col[5].i = i+1; 609 col[6].j = j+1; col[6].i = i-1; 610 col[7].j = j+1; col[7].i = i; 611 col[8].j = j+1; col[8].i = i+1; 612 v[0] = -0.25*(skew_S + skew_W); 613 v[1] = -hxdhy*newt_S + cross_EW; 614 v[2] = 0.25*(skew_S + skew_E); 615 v[3] = -hydhx*newt_W + cross_NS; 616 v[4] = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*PetscExpScalar(u); 617 v[5] = -hydhx*newt_E - cross_NS; 618 v[6] = 0.25*(skew_N + skew_W); 619 v[7] = -hxdhy*newt_N - cross_EW; 620 v[8] = -0.25*(skew_N + skew_E); 621 ierr = MatSetValuesStencil(B,1,&row,9,col,v,INSERT_VALUES);CHKERRQ(ierr); 622 break; 623 default: 624 SETERRQ1(PetscObjectComm((PetscObject)info->da),PETSC_ERR_SUP,"Jacobian type %d not implemented",user->jtype); 625 } 626 } 627 } 628 } 629 630 /* 631 Assemble matrix, using the 2-step process: 632 MatAssemblyBegin(), MatAssemblyEnd(). 633 */ 634 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 635 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 636 637 if (J != B) { 638 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 639 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 640 } 641 642 /* 643 Tell the matrix we will never add a new nonzero location to the 644 matrix. If we do, it will generate an error. 645 */ 646 if (user->jtype == JAC_NEWTON) { 647 ierr = PetscLogFlops(info->xm*info->ym*(131.0));CHKERRQ(ierr); 648 } 649 PetscFunctionReturn(0); 650 } 651 652 /*********************************************************** 653 * PreCheck implementation 654 ***********************************************************/ 655 PetscErrorCode PreCheckSetFromOptions(PreCheck precheck) 656 { 657 PetscErrorCode ierr; 658 PetscBool flg; 659 660 PetscFunctionBeginUser; 661 ierr = PetscOptionsBegin(precheck->comm,NULL,"PreCheck Options","none");CHKERRQ(ierr); 662 ierr = PetscOptionsReal("-precheck_angle","Angle in degrees between successive search directions necessary to activate step correction","",precheck->angle,&precheck->angle,NULL);CHKERRQ(ierr); 663 flg = PETSC_FALSE; 664 ierr = PetscOptionsBool("-precheck_monitor","Monitor choices made by precheck routine","",flg,&flg,NULL);CHKERRQ(ierr); 665 if (flg) { 666 ierr = PetscViewerASCIIOpen(precheck->comm,"stdout",&precheck->monitor);CHKERRQ(ierr); 667 } 668 ierr = PetscOptionsEnd();CHKERRQ(ierr); 669 PetscFunctionReturn(0); 670 } 671 672 /* 673 Compare the direction of the current and previous step, modify the current step accordingly 674 */ 675 PetscErrorCode PreCheckFunction(SNESLineSearch linesearch,Vec X,Vec Y,PetscBool *changed, void *ctx) 676 { 677 PetscErrorCode ierr; 678 PreCheck precheck; 679 Vec Ylast; 680 PetscScalar dot; 681 PetscInt iter; 682 PetscReal ynorm,ylastnorm,theta,angle_radians; 683 SNES snes; 684 685 PetscFunctionBeginUser; 686 ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 687 precheck = (PreCheck)ctx; 688 if (!precheck->Ylast) {ierr = VecDuplicate(Y,&precheck->Ylast);CHKERRQ(ierr);} 689 Ylast = precheck->Ylast; 690 ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr); 691 if (iter < 1) { 692 ierr = VecCopy(Y,Ylast);CHKERRQ(ierr); 693 *changed = PETSC_FALSE; 694 PetscFunctionReturn(0); 695 } 696 697 ierr = VecDot(Y,Ylast,&dot);CHKERRQ(ierr); 698 ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr); 699 ierr = VecNorm(Ylast,NORM_2,&ylastnorm);CHKERRQ(ierr); 700 /* Compute the angle between the vectors Y and Ylast, clip to keep inside the domain of acos() */ 701 theta = PetscAcosReal((PetscReal)PetscClipInterval(PetscAbsScalar(dot) / (ynorm * ylastnorm),-1.0,1.0)); 702 angle_radians = precheck->angle * PETSC_PI / 180.; 703 if (PetscAbsReal(theta) < angle_radians || PetscAbsReal(theta - PETSC_PI) < angle_radians) { 704 /* Modify the step Y */ 705 PetscReal alpha,ydiffnorm; 706 ierr = VecAXPY(Ylast,-1.0,Y);CHKERRQ(ierr); 707 ierr = VecNorm(Ylast,NORM_2,&ydiffnorm);CHKERRQ(ierr); 708 alpha = ylastnorm / ydiffnorm; 709 ierr = VecCopy(Y,Ylast);CHKERRQ(ierr); 710 ierr = VecScale(Y,alpha);CHKERRQ(ierr); 711 if (precheck->monitor) { 712 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); 713 } 714 } else { 715 ierr = VecCopy(Y,Ylast);CHKERRQ(ierr); 716 *changed = PETSC_FALSE; 717 if (precheck->monitor) { 718 ierr = PetscViewerASCIIPrintf(precheck->monitor,"Angle %E degrees exceeds threshold %g, no correction applied\n",(double)(theta*180./PETSC_PI),(double)precheck->angle);CHKERRQ(ierr); 719 } 720 } 721 PetscFunctionReturn(0); 722 } 723 724 PetscErrorCode PreCheckDestroy(PreCheck *precheck) 725 { 726 PetscErrorCode ierr; 727 728 PetscFunctionBeginUser; 729 if (!*precheck) PetscFunctionReturn(0); 730 ierr = VecDestroy(&(*precheck)->Ylast);CHKERRQ(ierr); 731 ierr = PetscViewerDestroy(&(*precheck)->monitor);CHKERRQ(ierr); 732 ierr = PetscFree(*precheck);CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 PetscErrorCode PreCheckCreate(MPI_Comm comm,PreCheck *precheck) 737 { 738 PetscErrorCode ierr; 739 740 PetscFunctionBeginUser; 741 ierr = PetscNew(precheck);CHKERRQ(ierr); 742 743 (*precheck)->comm = comm; 744 (*precheck)->angle = 10.; /* only active if angle is less than 10 degrees */ 745 PetscFunctionReturn(0); 746 } 747 748 /* 749 Applies some sweeps on nonlinear Gauss-Seidel on each process 750 751 */ 752 PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx) 753 { 754 PetscInt i,j,k,xs,ys,xm,ym,its,tot_its,sweeps,l,m; 755 PetscErrorCode ierr; 756 PetscReal hx,hy,hxdhy,hydhx,dhx,dhy,sc; 757 PetscScalar **x,**b,bij,F,F0=0,J,y,u,eu; 758 PetscReal atol,rtol,stol; 759 DM da; 760 AppCtx *user = (AppCtx*)ctx; 761 Vec localX,localB; 762 DMDALocalInfo info; 763 764 PetscFunctionBeginUser; 765 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 766 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 767 768 hx = 1.0/(PetscReal)(info.mx-1); 769 hy = 1.0/(PetscReal)(info.my-1); 770 sc = hx*hy*user->lambda; 771 dhx = 1/hx; 772 dhy = 1/hy; 773 hxdhy = hx/hy; 774 hydhx = hy/hx; 775 776 tot_its = 0; 777 ierr = SNESNGSGetSweeps(snes,&sweeps);CHKERRQ(ierr); 778 ierr = SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr); 779 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 780 if (B) { 781 ierr = DMGetLocalVector(da,&localB);CHKERRQ(ierr); 782 } 783 if (B) { 784 ierr = DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);CHKERRQ(ierr); 785 ierr = DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);CHKERRQ(ierr); 786 } 787 if (B) ierr = DMDAVecGetArrayRead(da,localB,&b);CHKERRQ(ierr); 788 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 789 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 790 ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr); 791 for (l=0; l<sweeps; l++) { 792 /* 793 Get local grid boundaries (for 2-dimensional DMDA): 794 xs, ys - starting grid indices (no ghost points) 795 xm, ym - widths of local grid (no ghost points) 796 */ 797 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 798 for (m=0; m<2; m++) { 799 for (j=ys; j<ys+ym; j++) { 800 for (i=xs+(m+j)%2; i<xs+xm; i+=2) { 801 PetscReal xx = i*hx,yy = j*hy; 802 if (B) bij = b[j][i]; 803 else bij = 0.; 804 805 if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) { 806 /* boundary conditions are all zero Dirichlet */ 807 x[j][i] = 0.0 + bij; 808 } else { 809 const PetscScalar 810 u_E = x[j][i+1], 811 u_W = x[j][i-1], 812 u_N = x[j+1][i], 813 u_S = x[j-1][i]; 814 const PetscScalar 815 uy_E = 0.25*dhy*(x[j+1][i]+x[j+1][i+1]-x[j-1][i]-x[j-1][i+1]), 816 uy_W = 0.25*dhy*(x[j+1][i-1]+x[j+1][i]-x[j-1][i-1]-x[j-1][i]), 817 ux_N = 0.25*dhx*(x[j][i+1]+x[j+1][i+1]-x[j][i-1]-x[j+1][i-1]), 818 ux_S = 0.25*dhx*(x[j-1][i+1]+x[j][i+1]-x[j-1][i-1]-x[j][i-1]); 819 u = x[j][i]; 820 for (k=0; k<its; k++) { 821 const PetscScalar 822 ux_E = dhx*(u_E-u), 823 ux_W = dhx*(u-u_W), 824 uy_N = dhy*(u_N-u), 825 uy_S = dhy*(u-u_S), 826 e_E = eta(user,xx,yy,ux_E,uy_E), 827 e_W = eta(user,xx,yy,ux_W,uy_W), 828 e_N = eta(user,xx,yy,ux_N,uy_N), 829 e_S = eta(user,xx,yy,ux_S,uy_S), 830 de_E = deta(user,xx,yy,ux_E,uy_E), 831 de_W = deta(user,xx,yy,ux_W,uy_W), 832 de_N = deta(user,xx,yy,ux_N,uy_N), 833 de_S = deta(user,xx,yy,ux_S,uy_S), 834 newt_E = e_E+de_E*PetscSqr(ux_E), 835 newt_W = e_W+de_W*PetscSqr(ux_W), 836 newt_N = e_N+de_N*PetscSqr(uy_N), 837 newt_S = e_S+de_S*PetscSqr(uy_S), 838 uxx = -hy * (e_E*ux_E - e_W*ux_W), 839 uyy = -hx * (e_N*uy_N - e_S*uy_S); 840 841 if (sc) eu = PetscExpScalar(u); 842 else eu = 0; 843 844 F = uxx + uyy - sc*eu - bij; 845 if (k == 0) F0 = F; 846 J = hxdhy*(newt_N + newt_S) + hydhx*(newt_E + newt_W) - sc*eu; 847 y = F/J; 848 u -= y; 849 tot_its++; 850 if (atol > PetscAbsReal(PetscRealPart(F)) || 851 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) || 852 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) { 853 break; 854 } 855 } 856 x[j][i] = u; 857 } 858 } 859 } 860 } 861 /* 862 x Restore vector 863 */ 864 } 865 ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr); 866 ierr = DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);CHKERRQ(ierr); 867 ierr = DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);CHKERRQ(ierr); 868 ierr = PetscLogFlops(tot_its*(118.0));CHKERRQ(ierr); 869 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 870 if (B) { 871 ierr = DMDAVecRestoreArrayRead(da,localB,&b);CHKERRQ(ierr); 872 ierr = DMRestoreLocalVector(da,&localB);CHKERRQ(ierr); 873 } 874 PetscFunctionReturn(0); 875 } 876 877 878 /*TEST 879 880 test: 881 nsize: 2 882 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype NEWTON 883 requires: !single 884 885 test: 886 suffix: 2 887 nsize: 2 888 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -precheck 1 889 requires: !single 890 891 test: 892 suffix: 3 893 nsize: 2 894 args: -snes_monitor_short -da_grid_x 20 -da_grid_y 20 -p 1.3 -lambda 1 -jtype PICARD -picard -precheck 1 -p 3 895 requires: !single 896 897 test: 898 suffix: 3_star 899 nsize: 2 900 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 901 output_file: output/ex15_3.out 902 requires: !single 903 904 test: 905 suffix: 4 906 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 907 requires: !single 908 909 test: 910 suffix: lag_jac 911 nsize: 4 912 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 913 requires: !single 914 915 test: 916 suffix: lag_pc 917 nsize: 4 918 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 919 requires: !single 920 921 test: 922 suffix: nleqerr 923 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 924 requires: !single 925 926 test: 927 suffix: mf 928 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator 929 requires: !single 930 931 test: 932 suffix: mf_picard 933 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_mf_operator -picard 934 requires: !single 935 output_file: output/ex15_mf.out 936 937 test: 938 suffix: fd_picard 939 args: -snes_monitor_short -pc_type lu -da_refine 2 -p 3 -ksp_rtol 1.e-12 -snes_fd -picard 940 requires: !single 941 942 test: 943 suffix: fd_color_picard 944 args: -snes_monitor_short -pc_type lu -da_refine 4 -p 3 -ksp_rtol 1.e-12 -snes_fd_color -picard 945 requires: !single 946 output_file: output/ex15_mf.out 947 948 TEST*/ 949