Actual source code: ex5.c

petsc-master 2018-01-14
Report Typos and Errors

  2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
  3: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
  4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
  5: The command line options include:\n\
  6:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  7:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n\
  8:   -m_par/n_par <parameter>, where <parameter> indicates an integer\n \
  9:       that MMS3 will be evaluated with 2^m_par, 2^n_par";

 11: /*T
 12:    Concepts: SNES^parallel Bratu example
 13:    Concepts: DMDA^using distributed arrays;
 14:    Concepts: IS coloirng types;
 15:    Processors: n
 16: T*/



 20: /* ------------------------------------------------------------------------

 22:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 23:     the partial differential equation

 25:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,

 27:     with boundary conditions

 29:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

 31:     A finite difference approximation with the usual 5-point stencil
 32:     is used to discretize the boundary value problem to obtain a nonlinear
 33:     system of equations.


 36:       This example shows how geometric multigrid can be run transparently with a nonlinear solver so long
 37:       as SNESSetDM() is provided. Example usage

 39:       ./ex5 -pc_type mg -ksp_monitor  -snes_view -pc_mg_levels 3 -pc_mg_galerkin pmat -da_grid_x 17 -da_grid_y 17
 40:              -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full

 42:       or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
 43:          multigrid levels, it will be determined automatically based on the number of refinements done)

 45:       ./ex5 -pc_type mg -ksp_monitor  -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3
 46:              -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full

 48:   ------------------------------------------------------------------------- */

 50: /*
 51:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 52:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 53: */
 54:  #include <petscdm.h>
 55:  #include <petscdmda.h>
 56:  #include <petscsnes.h>
 57:  #include <petscmatlab.h>

 59: /*
 60:    User-defined application context - contains data needed by the
 61:    application-provided call-back routines, FormJacobianLocal() and
 62:    FormFunctionLocal().
 63: */
 64: typedef struct AppCtx AppCtx;
 65: struct AppCtx {
 66:   PetscReal param;          /* test problem parameter */
 67:   PetscInt  m,n;            /* MMS3 parameters */
 68:   PetscErrorCode (*mms_solution)(AppCtx*,const DMDACoor2d*,PetscScalar*);
 69:   PetscErrorCode (*mms_forcing)(AppCtx*,const DMDACoor2d*,PetscScalar*);
 70: };

 72: /*
 73:    User-defined routines
 74: */
 75: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
 76: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 77: extern PetscErrorCode FormExactSolution(DM,AppCtx*,Vec);
 78: extern PetscErrorCode ZeroBCSolution(AppCtx*,const DMDACoor2d*,PetscScalar*);
 79: extern PetscErrorCode MMSSolution1(AppCtx*,const DMDACoor2d*,PetscScalar*);
 80: extern PetscErrorCode MMSForcing1(AppCtx*,const DMDACoor2d*,PetscScalar*);
 81: extern PetscErrorCode MMSSolution2(AppCtx*,const DMDACoor2d*,PetscScalar*);
 82: extern PetscErrorCode MMSForcing2(AppCtx*,const DMDACoor2d*,PetscScalar*);
 83: extern PetscErrorCode MMSSolution3(AppCtx*,const DMDACoor2d*,PetscScalar*);
 84: extern PetscErrorCode MMSForcing3(AppCtx*,const DMDACoor2d*,PetscScalar*);
 85: extern PetscErrorCode MMSSolution4(AppCtx*,const DMDACoor2d*,PetscScalar*);
 86: extern PetscErrorCode MMSForcing4(AppCtx*,const DMDACoor2d*,PetscScalar*);
 87: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
 88: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
 89: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 90: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
 91: #endif
 92: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);

 94: int main(int argc,char **argv)
 95: {
 96:   SNES           snes;                         /* nonlinear solver */
 97:   Vec            x;                            /* solution vector */
 98:   AppCtx         user;                         /* user-defined work context */
 99:   PetscInt       its;                          /* iterations for convergence */
101:   PetscReal      bratu_lambda_max = 6.81;
102:   PetscReal      bratu_lambda_min = 0.;
103:   PetscInt       MMS              = 0;
104:   PetscBool      flg              = PETSC_FALSE;
105:   DM             da;
106: #if defined(PETSC_HAVE_MATLAB_ENGINE)
107:   Vec            r               = NULL;
108:   PetscBool      matlab_function = PETSC_FALSE;
109: #endif
110:   KSP            ksp;
111:   PetscInt       lits,slits;

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Initialize program
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

117:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Initialize problem parameters
121:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122:   user.param = 6.0;
123:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
124:   if (user.param > bratu_lambda_max || user.param < bratu_lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda, %g, is out of range, [%g, %g]", user.param, bratu_lambda_min, bratu_lambda_max);
125:   PetscOptionsGetInt(NULL,NULL,"-mms",&MMS,NULL);
126:   if (MMS == 3) {
127:     PetscInt mPar = 2, nPar = 1;
128:     PetscOptionsGetInt(NULL,NULL,"-m_par",&mPar,NULL);
129:     PetscOptionsGetInt(NULL,NULL,"-n_par",&nPar,NULL);
130:     user.m = PetscPowInt(2,mPar);
131:     user.n = PetscPowInt(2,nPar);
132:   }

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Create nonlinear solver context
136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   SNESCreate(PETSC_COMM_WORLD,&snes);
138:   SNESSetCountersReset(snes,PETSC_FALSE);
139:   SNESSetNGS(snes, NonlinearGS, NULL);

141:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142:      Create distributed array (DMDA) to manage parallel grid and vectors
143:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
145:   DMSetFromOptions(da);
146:   DMSetUp(da);
147:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
148:   DMSetApplicationContext(da,&user);
149:   SNESSetDM(snes,da);
150:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:      Extract global vectors from DMDA; then duplicate for remaining
152:      vectors that are the same types
153:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154:   DMCreateGlobalVector(da,&x);

156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157:      Set local function evaluation routine
158:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159:   user.mms_solution = ZeroBCSolution;
160:   switch (MMS) {
161:   case 0: user.mms_solution = NULL; user.mms_forcing = NULL;
162:   case 1: user.mms_solution = MMSSolution1; user.mms_forcing = MMSForcing1; break;
163:   case 2: user.mms_solution = MMSSolution2; user.mms_forcing = MMSForcing2; break;
164:   case 3: user.mms_solution = MMSSolution3; user.mms_forcing = MMSForcing3; break;
165:   case 4: user.mms_solution = MMSSolution4; user.mms_forcing = MMSForcing4; break;
166:   default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unknown MMS type %d",MMS);
167:   }
168:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
169:   PetscOptionsGetBool(NULL,NULL,"-fd",&flg,NULL);
170:   if (!flg) {
171:     DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
172:   }

174:   PetscOptionsGetBool(NULL,NULL,"-obj",&flg,NULL);
175:   if (flg) {
176:     DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
177:   }

179: #if defined(PETSC_HAVE_MATLAB_ENGINE)
180:   PetscOptionsGetBool(NULL,NULL,"-matlab_function",&matlab_function,0);
181:   if (matlab_function) {
182:     VecDuplicate(x,&r);
183:     SNESSetFunction(snes,r,FormFunctionMatlab,&user);
184:   }
185: #endif

187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188:      Customize nonlinear solver; set runtime options
189:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190:   SNESSetFromOptions(snes);

192:   FormInitialGuess(da,&user,x);

194:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195:      Solve nonlinear system
196:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197:   SNESSolve(snes,NULL,x);
198:   SNESGetIterationNumber(snes,&its);

200:   SNESGetLinearSolveIterations(snes,&slits);
201:   SNESGetKSP(snes,&ksp);
202:   KSPGetTotalIterations(ksp,&lits);
203:   if (lits != slits) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Number of total linear iterations reported by SNES %D does not match reported by KSP %D",slits,lits);
204:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

207:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208:      If using MMS, check the l_2 error
209:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
210:   if (MMS) {
211:     Vec       e;
212:     PetscReal errorl2, errorinf;
213:     PetscInt  N;

215:     VecDuplicate(x, &e);
216:     PetscObjectViewFromOptions((PetscObject) x, NULL, "-sol_view");
217:     FormExactSolution(da, &user, e);
218:     PetscObjectViewFromOptions((PetscObject) e, NULL, "-exact_view");
219:     VecAXPY(e, -1.0, x);
220:     PetscObjectViewFromOptions((PetscObject) e, NULL, "-error_view");
221:     VecNorm(e, NORM_2, &errorl2);
222:     VecNorm(e, NORM_INFINITY, &errorinf);
223:     VecGetSize(e, &N);
224:     PetscPrintf(PETSC_COMM_WORLD, "N: %D error L2 %g inf %g\n", N, (double) errorl2/PetscSqrtReal(N), (double) errorinf);
225:     VecDestroy(&e);
226:   }

228:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229:      Free work space.  All PETSc objects should be destroyed when they
230:      are no longer needed.
231:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232: #if defined(PETSC_HAVE_MATLAB_ENGINE)
233:   VecDestroy(&r);
234: #endif
235:   VecDestroy(&x);
236:   SNESDestroy(&snes);
237:   DMDestroy(&da);
238:   PetscFinalize();
239:   return ierr;
240: }
241: /* ------------------------------------------------------------------- */
242: /*
243:    FormInitialGuess - Forms initial approximation.

245:    Input Parameters:
246:    da - The DM
247:    user - user-defined application context

249:    Output Parameter:
250:    X - vector
251:  */
252: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
253: {
254:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
256:   PetscReal      lambda,temp1,temp,hx,hy;
257:   PetscScalar    **x;

260:   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);

262:   lambda = user->param;
263:   hx     = 1.0/(PetscReal)(Mx-1);
264:   hy     = 1.0/(PetscReal)(My-1);
265:   temp1  = lambda/(lambda + 1.0);

267:   /*
268:      Get a pointer to vector data.
269:        - For default PETSc vectors, VecGetArray() returns a pointer to
270:          the data array.  Otherwise, the routine is implementation dependent.
271:        - You MUST call VecRestoreArray() when you no longer need access to
272:          the array.
273:   */
274:   DMDAVecGetArray(da,X,&x);

276:   /*
277:      Get local grid boundaries (for 2-dimensional DMDA):
278:        xs, ys   - starting grid indices (no ghost points)
279:        xm, ym   - widths of local grid (no ghost points)

281:   */
282:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

284:   /*
285:      Compute initial guess over the locally owned part of the grid
286:   */
287:   for (j=ys; j<ys+ym; j++) {
288:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
289:     for (i=xs; i<xs+xm; i++) {
290:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
291:         /* boundary conditions are all zero Dirichlet */
292:         x[j][i] = 0.0;
293:       } else {
294:         x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
295:       }
296:     }
297:   }

299:   /*
300:      Restore vector
301:   */
302:   DMDAVecRestoreArray(da,X,&x);
303:   return(0);
304: }

306: /*
307:   FormExactSolution - Forms MMS solution

309:   Input Parameters:
310:   da - The DM
311:   user - user-defined application context

313:   Output Parameter:
314:   X - vector
315:  */
316: PetscErrorCode FormExactSolution(DM da, AppCtx *user, Vec U)
317: {
318:   DM             coordDA;
319:   Vec            coordinates;
320:   DMDACoor2d   **coords;
321:   PetscScalar  **u;
322:   PetscInt       xs, ys, xm, ym, i, j;

326:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
327:   DMGetCoordinateDM(da, &coordDA);
328:   DMGetCoordinates(da, &coordinates);
329:   DMDAVecGetArray(coordDA, coordinates, &coords);
330:   DMDAVecGetArray(da, U, &u);
331:   for (j = ys; j < ys+ym; ++j) {
332:     for (i = xs; i < xs+xm; ++i) {
333:       user->mms_solution(user,&coords[j][i],&u[j][i]);
334:     }
335:   }
336:   DMDAVecRestoreArray(da, U, &u);
337:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
338:   return(0);
339: }

341: PetscErrorCode ZeroBCSolution(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
342: {
343:   u[0] = 0.;
344:   return 0;
345: }

347: /* The functions below evaluate the MMS solution u(x,y) and associated forcing

349:      f(x,y) = -u_xx - u_yy - lambda exp(u)

351:   such that u(x,y) is an exact solution with f(x,y) as the right hand side forcing term.
352:  */
353: PetscErrorCode MMSSolution1(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
354: {
355:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
356:   u[0] = x*(1 - x)*y*(1 - y);
357:   PetscLogFlops(5);
358:   return 0;
359: }
360: PetscErrorCode MMSForcing1(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
361: {
362:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
363:   f[0] = 2*x*(1 - x) + 2*y*(1 - y) - user->param*PetscExpReal(x*(1 - x)*y*(1 - y));
364:   return 0;
365: }

367: PetscErrorCode MMSSolution2(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
368: {
369:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
370:   u[0] = PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y);
371:   PetscLogFlops(5);
372:   return 0;
373: }
374: PetscErrorCode MMSForcing2(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
375: {
376:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
377:   f[0] = 2*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y) - user->param*PetscExpReal(PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y));
378:   return 0;
379: }

381: PetscErrorCode MMSSolution3(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
382: {
383:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
384:   u[0] = PetscSinReal(user->m*PETSC_PI*x*(1-y))*PetscSinReal(user->n*PETSC_PI*y*(1-x));
385:   PetscLogFlops(5);
386:   return 0;
387: }
388: PetscErrorCode MMSForcing3(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
389: {
390:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
391:   PetscReal m = user->m, n = user->n, lambda = user->param;
392:   f[0] = (-(PetscExpReal(PetscSinReal(m*PETSC_PI*x*(1 - y))*PetscSinReal(n*PETSC_PI*(1 - x)*y))*lambda)
393:           + PetscSqr(PETSC_PI)*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*PetscCosReal(m*PETSC_PI*x*(-1 + y))*PetscCosReal(n*PETSC_PI*(-1 + x)*y)
394:                                 + (PetscSqr(m)*(PetscSqr(x) + PetscSqr(-1 + y)) + PetscSqr(n)*(PetscSqr(-1 + x) + PetscSqr(y)))
395:                                 *PetscSinReal(m*PETSC_PI*x*(-1 + y))*PetscSinReal(n*PETSC_PI*(-1 + x)*y)));
396:   return 0;
397: }

399: PetscErrorCode MMSSolution4(AppCtx *user,const DMDACoor2d *c,PetscScalar *u)
400: {
401:   const PetscReal Lx = 1.,Ly = 1.;
402:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
403:   u[0] = (PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y));
404:   PetscLogFlops(9);
405:   return 0;
406: }
407: PetscErrorCode MMSForcing4(AppCtx *user,const DMDACoor2d *c,PetscScalar *f)
408: {
409:   const PetscReal Lx = 1.,Ly = 1.;
410:   PetscReal x = PetscRealPart(c->x), y = PetscRealPart(c->y);
411:   f[0] = (2*PetscSqr(x)*(PetscSqr(x)-PetscSqr(Lx))*(PetscSqr(Ly)-6*PetscSqr(y))
412:           + 2*PetscSqr(y)*(PetscSqr(Lx)-6*PetscSqr(x))*(PetscSqr(y)-PetscSqr(Ly))
413:           - user->param*PetscExpReal((PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y))));
414:   return 0;
415: }

417: /* ------------------------------------------------------------------- */
418: /*
419:    FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch


422:  */
423: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
424: {
426:   PetscInt       i,j;
427:   PetscReal      lambda,hx,hy,hxdhy,hydhx;
428:   PetscScalar    u,ue,uw,un,us,uxx,uyy,mms_solution,mms_forcing;
429:   DMDACoor2d     c;

432:   lambda = user->param;
433:   hx     = 1.0/(PetscReal)(info->mx-1);
434:   hy     = 1.0/(PetscReal)(info->my-1);
435:   hxdhy  = hx/hy;
436:   hydhx  = hy/hx;
437:   /*
438:      Compute function over the locally owned part of the grid
439:   */
440:   for (j=info->ys; j<info->ys+info->ym; j++) {
441:     for (i=info->xs; i<info->xs+info->xm; i++) {
442:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
443:         c.x = i*hx; c.y = j*hy;
444:         user->mms_solution(user,&c,&mms_solution);
445:         f[j][i] = 2.0*(hydhx+hxdhy)*(x[j][i] - mms_solution);
446:       } else {
447:         u  = x[j][i];
448:         uw = x[j][i-1];
449:         ue = x[j][i+1];
450:         un = x[j-1][i];
451:         us = x[j+1][i];

453:         /* Enforce boundary conditions at neighboring points -- setting these values causes the Jacobian to be symmetric. */
454:         if (i-1 == 0) {c.x = (i-1)*hx; c.y = j*hy; user->mms_solution(user,&c,&uw);}
455:         if (i+1 == info->mx-1) {c.x = (i+1)*hx; c.y = j*hy; user->mms_solution(user,&c,&ue);}
456:         if (j-1 == 0) {c.x = i*hx; c.y = (j-1)*hy; user->mms_solution(user,&c,&un);}
457:         if (j+1 == info->my-1) {c.x = i*hx; c.y = (j+1)*hy; user->mms_solution(user,&c,&us);}

459:         uxx     = (2.0*u - uw - ue)*hydhx;
460:         uyy     = (2.0*u - un - us)*hxdhy;
461:         mms_forcing = 0;
462:         c.x = i*hx; c.y = j*hy;
463:         if (user->mms_forcing) {user->mms_forcing(user,&c,&mms_forcing);}
464:         f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + mms_forcing);
465:       }
466:     }
467:   }
468:   PetscLogFlops(11.0*info->ym*info->xm);
469:   return(0);
470: }

472: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
473: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
474: {
476:   PetscInt       i,j;
477:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
478:   PetscScalar    u,ue,uw,un,us,uxux,uyuy;
479:   MPI_Comm       comm;

482:   *obj   = 0;
483:   PetscObjectGetComm((PetscObject)info->da,&comm);
484:   lambda = user->param;
485:   hx     = 1.0/(PetscReal)(info->mx-1);
486:   hy     = 1.0/(PetscReal)(info->my-1);
487:   sc     = hx*hy*lambda;
488:   hxdhy  = hx/hy;
489:   hydhx  = hy/hx;
490:   /*
491:      Compute function over the locally owned part of the grid
492:   */
493:   for (j=info->ys; j<info->ys+info->ym; j++) {
494:     for (i=info->xs; i<info->xs+info->xm; i++) {
495:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
496:         lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
497:       } else {
498:         u  = x[j][i];
499:         uw = x[j][i-1];
500:         ue = x[j][i+1];
501:         un = x[j-1][i];
502:         us = x[j+1][i];

504:         if (i-1 == 0) uw = 0.;
505:         if (i+1 == info->mx-1) ue = 0.;
506:         if (j-1 == 0) un = 0.;
507:         if (j+1 == info->my-1) us = 0.;

509:         /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */

511:         uxux = u*(2.*u - ue - uw)*hydhx;
512:         uyuy = u*(2.*u - un - us)*hxdhy;

514:         lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
515:       }
516:     }
517:   }
518:   PetscLogFlops(12.0*info->ym*info->xm);
519:   MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
520:   return(0);
521: }

523: /*
524:    FormJacobianLocal - Evaluates Jacobian matrix on local process patch
525: */
526: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
527: {
529:   PetscInt       i,j,k;
530:   MatStencil     col[5],row;
531:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;
532:   DM             coordDA;
533:   Vec            coordinates;
534:   DMDACoor2d   **coords;

537:   lambda = user->param;
538:   /* Extract coordinates */
539:   DMGetCoordinateDM(info->da, &coordDA);
540:   DMGetCoordinates(info->da, &coordinates);
541:   DMDAVecGetArray(coordDA, coordinates, &coords);
542:   hx     = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
543:   hy     = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
544:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
545:   hxdhy  = hx/hy;
546:   hydhx  = hy/hx;
547:   sc     = hx*hy*lambda;


550:   /*
551:      Compute entries for the locally owned part of the Jacobian.
552:       - Currently, all PETSc parallel matrix formats are partitioned by
553:         contiguous chunks of rows across the processors.
554:       - Each processor needs to insert only elements that it owns
555:         locally (but any non-local elements will be sent to the
556:         appropriate processor during matrix assembly).
557:       - Here, we set all entries for a particular row at once.
558:       - We can set matrix entries either using either
559:         MatSetValuesLocal() or MatSetValues(), as discussed above.
560:   */
561:   for (j=info->ys; j<info->ys+info->ym; j++) {
562:     for (i=info->xs; i<info->xs+info->xm; i++) {
563:       row.j = j; row.i = i;
564:       /* boundary points */
565:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
566:         v[0] =  2.0*(hydhx + hxdhy);
567:         MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
568:       } else {
569:         k = 0;
570:         /* interior grid points */
571:         if (j-1 != 0) {
572:           v[k]     = -hxdhy;
573:           col[k].j = j - 1; col[k].i = i;
574:           k++;
575:         }
576:         if (i-1 != 0) {
577:           v[k]     = -hydhx;
578:           col[k].j = j;     col[k].i = i-1;
579:           k++;
580:         }

582:         v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;

584:         if (i+1 != info->mx-1) {
585:           v[k]     = -hydhx;
586:           col[k].j = j;     col[k].i = i+1;
587:           k++;
588:         }
589:         if (j+1 != info->mx-1) {
590:           v[k]     = -hxdhy;
591:           col[k].j = j + 1; col[k].i = i;
592:           k++;
593:         }
594:         MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
595:       }
596:     }
597:   }

599:   /*
600:      Assemble matrix, using the 2-step process:
601:        MatAssemblyBegin(), MatAssemblyEnd().
602:   */
603:   MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
604:   MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
605:   /*
606:      Tell the matrix we will never add a new nonzero location to the
607:      matrix. If we do, it will generate an error.
608:   */
609:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
610:   return(0);
611: }

613: #if defined(PETSC_HAVE_MATLAB_ENGINE)
614: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
615: {
616:   AppCtx         *user = (AppCtx*)ptr;
618:   PetscInt       Mx,My;
619:   PetscReal      lambda,hx,hy;
620:   Vec            localX,localF;
621:   MPI_Comm       comm;
622:   DM             da;

625:   SNESGetDM(snes,&da);
626:   DMGetLocalVector(da,&localX);
627:   DMGetLocalVector(da,&localF);
628:   PetscObjectSetName((PetscObject)localX,"localX");
629:   PetscObjectSetName((PetscObject)localF,"localF");
630:   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);

632:   lambda = user->param;
633:   hx     = 1.0/(PetscReal)(Mx-1);
634:   hy     = 1.0/(PetscReal)(My-1);

636:   PetscObjectGetComm((PetscObject)snes,&comm);
637:   /*
638:      Scatter ghost points to local vector,using the 2-step process
639:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
640:      By placing code between these two statements, computations can be
641:      done while messages are in transition.
642:   */
643:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
644:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
645:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
646:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
647:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

649:   /*
650:      Insert values into global vector
651:   */
652:   DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
653:   DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
654:   DMRestoreLocalVector(da,&localX);
655:   DMRestoreLocalVector(da,&localF);
656:   return(0);
657: }
658: #endif

660: /* ------------------------------------------------------------------- */
661: /*
662:       Applies some sweeps on nonlinear Gauss-Seidel on each process

664:  */
665: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
666: {
667:   PetscInt       i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
669:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
670:   PetscScalar    **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
671:   PetscReal      atol,rtol,stol;
672:   DM             da;
673:   AppCtx         *user;
674:   Vec            localX,localB;

677:   tot_its = 0;
678:   SNESNGSGetSweeps(snes,&sweeps);
679:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
680:   SNESGetDM(snes,&da);
681:   DMGetApplicationContext(da,(void**)&user);

683:   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);

685:   lambda = user->param;
686:   hx     = 1.0/(PetscReal)(Mx-1);
687:   hy     = 1.0/(PetscReal)(My-1);
688:   sc     = hx*hy*lambda;
689:   hxdhy  = hx/hy;
690:   hydhx  = hy/hx;


693:   DMGetLocalVector(da,&localX);
694:   if (B) {
695:     DMGetLocalVector(da,&localB);
696:   }
697:   for (l=0; l<sweeps; l++) {

699:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
700:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
701:     if (B) {
702:       DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
703:       DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
704:     }
705:     /*
706:      Get a pointer to vector data.
707:      - For default PETSc vectors, VecGetArray() returns a pointer to
708:      the data array.  Otherwise, the routine is implementation dependent.
709:      - You MUST call VecRestoreArray() when you no longer need access to
710:      the array.
711:      */
712:     DMDAVecGetArray(da,localX,&x);
713:     if (B) DMDAVecGetArray(da,localB,&b);
714:     /*
715:      Get local grid boundaries (for 2-dimensional DMDA):
716:      xs, ys   - starting grid indices (no ghost points)
717:      xm, ym   - widths of local grid (no ghost points)
718:      */
719:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

721:     for (j=ys; j<ys+ym; j++) {
722:       for (i=xs; i<xs+xm; i++) {
723:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
724:           /* boundary conditions are all zero Dirichlet */
725:           x[j][i] = 0.0;
726:         } else {
727:           if (B) bij = b[j][i];
728:           else   bij = 0.;

730:           u  = x[j][i];
731:           un = x[j-1][i];
732:           us = x[j+1][i];
733:           ue = x[j][i-1];
734:           uw = x[j][i+1];

736:           for (k=0; k<its; k++) {
737:             eu  = PetscExpScalar(u);
738:             uxx = (2.0*u - ue - uw)*hydhx;
739:             uyy = (2.0*u - un - us)*hxdhy;
740:             F   = uxx + uyy - sc*eu - bij;
741:             if (k == 0) F0 = F;
742:             J  = 2.0*(hydhx + hxdhy) - sc*eu;
743:             y  = F/J;
744:             u -= y;
745:             tot_its++;

747:             if (atol > PetscAbsReal(PetscRealPart(F)) ||
748:                 rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
749:                 stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
750:               break;
751:             }
752:           }
753:           x[j][i] = u;
754:         }
755:       }
756:     }
757:     /*
758:      Restore vector
759:      */
760:     DMDAVecRestoreArray(da,localX,&x);
761:     DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
762:     DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
763:   }
764:   PetscLogFlops(tot_its*(21.0));
765:   DMRestoreLocalVector(da,&localX);
766:   if (B) {
767:     DMDAVecRestoreArray(da,localB,&b);
768:     DMRestoreLocalVector(da,&localB);
769:   }
770:   return(0);
771: }

773: /*TEST

775:    test:
776:      suffix: asm_0
777:      requires: !single
778:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu

780:    test:
781:      suffix: msm_0
782:      requires: !single
783:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu

785:    test:
786:      suffix: asm_1
787:      requires: !single
788:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8

790:    test:
791:      suffix: msm_1
792:      requires: !single
793:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8

795:    test:
796:      suffix: asm_2
797:      requires: !single
798:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 3 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8

800:    test:
801:      suffix: msm_2
802:      requires: !single
803:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 3 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8

805:    test:
806:      suffix: asm_3
807:      requires: !single
808:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8

810:    test:
811:      suffix: msm_3
812:      requires: !single
813:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8

815:    test:
816:      suffix: asm_4
817:      requires: !single
818:      nsize: 2
819:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8

821:    test:
822:      suffix: msm_4
823:      requires: !single
824:      nsize: 2
825:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 2 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8

827:    test:
828:      suffix: asm_5
829:      requires: !single
830:      nsize: 2
831:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type additive -sub_pc_type lu -da_grid_x 8

833:    test:
834:      suffix: msm_5
835:      requires: !single
836:      nsize: 2
837:      args: -mms 1 -par 0.0 -snes_monitor_short -snes_converged_reason -snes_view -ksp_rtol 1.0e-9 -ksp_monitor_short -ksp_type richardson -pc_type asm -pc_asm_blocks 4 -pc_asm_overlap 0 -pc_asm_local_type multiplicative -sub_pc_type lu -da_grid_x 8

839:    test:
840:      args: -snes_rtol 1.e-5 -pc_type mg -ksp_monitor_short -snes_view -pc_mg_levels 3 -pc_mg_galerkin pmat -da_grid_x 17 -da_grid_y 17 -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -snes_monitor_short -mg_levels_ksp_chebyshev_esteig 0.5,1.1 -mg_levels_pc_type sor -pc_mg_type full
841:      requires: !single

843:    test:
844:      suffix: 2
845:      args: -pc_type mg -ksp_converged_reason -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3 -mg_levels_ksp_norm_type unpreconditioned -snes_monitor_short -mg_levels_ksp_chebyshev_esteig 0.5,1.1 -mg_levels_pc_type sor -pc_mg_type full -ksp_atol -1.

847:    test:
848:      suffix: 3
849:      nsize: 2
850:      args: -snes_grid_sequence 2 -snes_mf_operator -snes_converged_reason -snes_view -pc_type mg -snes_atol -1 -snes_rtol 1.e-2
851:      filter: grep -v "otal number of function evaluations"

853:    test:
854:      suffix: 4
855:      nsize: 2
856:      args: -snes_grid_sequence 2 -snes_monitor_short -ksp_converged_reason -snes_converged_reason -snes_view -pc_type mg -snes_atol -1 -ksp_atol -1

858:    test:
859:      suffix: 5_anderson
860:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type anderson

862:    test:
863:      suffix: 5_aspin
864:      nsize: 4
865:      args: -snes_monitor_short -ksp_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type aspin

867:    test:
868:      suffix: 5_broyden
869:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type qn -snes_qn_type broyden -snes_qn_m 10

871:    test:
872:      suffix: 5_fas
873:      args: -fas_coarse_snes_max_it 1 -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_monitor_short -snes_type fas -fas_coarse_ksp_type richardson -da_refine 6
874:      requires: !single

876:    test:
877:      suffix: 5_fas_additive
878:      args: -fas_coarse_snes_max_it 1 -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_monitor_short -snes_type fas -fas_coarse_ksp_type richardson -da_refine 6 -snes_fas_type additive -snes_max_it 50

880:    test:
881:      suffix: 5_fas_monitor
882:      args: -da_refine 1 -snes_type fas -snes_fas_monitor
883:      requires: !single

885:    test:
886:      suffix: 5_ls
887:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type newtonls

889:    test:
890:      suffix: 5_nasm
891:      nsize: 4
892:      args: -snes_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type nasm -snes_nasm_type restrict -snes_max_it 10

894:    test:
895:      suffix: 5_ncg
896:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ncg -snes_ncg_type fr

898:    test:
899:      suffix: 5_newton_asm_dmda
900:      nsize: 4
901:      args: -snes_monitor_short -ksp_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type newtonls -pc_type asm -pc_asm_dm_subdomains -malloc_dump
902:      requires: !single

904:    test:
905:      suffix: 5_newton_gasm_dmda
906:      nsize: 4
907:      args: -snes_monitor_short -ksp_monitor_short -snes_converged_reason -da_refine 4 -da_overlap 3 -snes_type newtonls -pc_type gasm -malloc_dump
908:      requires: !single

910:    test:
911:      suffix: 5_ngmres
912:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ngmres -snes_ngmres_m 10

914:    test:
915:      suffix: 5_ngmres_fas
916:      args: -snes_rtol 1.e-4 -snes_type ngmres -npc_fas_coarse_snes_max_it 1 -npc_fas_coarse_snes_type newtonls -npc_fas_coarse_pc_type lu -npc_fas_coarse_ksp_type preonly -snes_ngmres_m 10 -snes_monitor_short -npc_snes_max_it 1 -npc_snes_type fas -npc_fas_coarse_ksp_type richardson -da_refine 6

918:    test:
919:      suffix: 5_ngmres_ngs
920:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ngmres -npc_snes_type ngs -npc_snes_max_it 1

922:    test:
923:      suffix: 5_ngmres_nrichardson
924:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type ngmres -snes_ngmres_m 10 -npc_snes_type nrichardson -npc_snes_max_it 3
925:      output_file: output/ex5_5_ngmres_richardson.out

927:    test:
928:      suffix: 5_nrichardson
929:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type nrichardson

931:    test:
932:      suffix: 5_qn
933:      args: -da_grid_x 81 -da_grid_y 81 -snes_monitor_short -snes_max_it 50 -par 6.0 -snes_type qn -snes_linesearch_type cp -snes_qn_m 10

935:    test:
936:      suffix: 6
937:      nsize: 4
938:      args: -snes_converged_reason -ksp_converged_reason -da_grid_x 129 -da_grid_y 129 -pc_type mg -pc_mg_levels 8 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.5,0,1.1 -mg_levels_ksp_max_it 2

940: TEST*/