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