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