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