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