xref: /petsc/src/snes/tutorials/ex14.c (revision 76ce8f1ac04011f726b87b0654d63ff8fc64561d)
1 
2 static char help[] = "Bratu nonlinear PDE in 3d.\n\
3 We solve the  Bratu (SFI - solid fuel ignition) problem in a 3D 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 
9 /*T
10    Concepts: SNES^parallel Bratu example
11    Concepts: DMDA^using distributed arrays;
12    Processors: n
13 T*/
14 
15 /* ------------------------------------------------------------------------
16 
17     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
18     the partial differential equation
19 
20             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
21 
22     with boundary conditions
23 
24              u = 0  for  x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
25 
26     A finite difference approximation with the usual 7-point stencil
27     is used to discretize the boundary value problem to obtain a nonlinear
28     system of equations.
29 
30   ------------------------------------------------------------------------- */
31 
32 /*
33    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
34    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
35    file automatically includes:
36      petscsys.h       - base PETSc routines   petscvec.h - vectors
37      petscmat.h - matrices
38      petscis.h     - index sets            petscksp.h - Krylov subspace methods
39      petscviewer.h - viewers               petscpc.h  - preconditioners
40      petscksp.h   - linear solvers
41 */
42 #include <petscdm.h>
43 #include <petscdmda.h>
44 #include <petscsnes.h>
45 
46 /*
47    User-defined application context - contains data needed by the
48    application-provided call-back routines, FormJacobian() and
49    FormFunction().
50 */
51 typedef struct {
52   PetscReal param;             /* test problem parameter */
53   DM        da;                /* distributed array data structure */
54 } AppCtx;
55 
56 /*
57    User-defined routines
58 */
59 extern PetscErrorCode FormFunctionLocal(SNES,Vec,Vec,void*);
60 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
61 extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
62 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
63 
64 int main(int argc,char **argv)
65 {
66   SNES           snes;                         /* nonlinear solver */
67   Vec            x,r;                          /* solution, residual vectors */
68   Mat            J = NULL;                            /* Jacobian matrix */
69   AppCtx         user;                         /* user-defined work context */
70   PetscInt       its;                          /* iterations for convergence */
71   MatFDColoring  matfdcoloring = NULL;
72   PetscBool      matrix_free = PETSC_FALSE,coloring = PETSC_FALSE, coloring_ds = PETSC_FALSE,local_coloring = PETSC_FALSE;
73   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm;
74 
75   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76      Initialize program
77      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78 
79   PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
80 
81   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82      Initialize problem parameters
83   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84   user.param = 6.0;
85   PetscCall(PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL));
86   PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lambda is out of range");
87 
88   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89      Create nonlinear solver context
90      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91   PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
92 
93   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94      Create distributed array (DMDA) to manage parallel grid and vectors
95   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96   PetscCall(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,4,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&user.da));
97   PetscCall(DMSetFromOptions(user.da));
98   PetscCall(DMSetUp(user.da));
99   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100      Extract global vectors from DMDA; then duplicate for remaining
101      vectors that are the same types
102    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103   PetscCall(DMCreateGlobalVector(user.da,&x));
104   PetscCall(VecDuplicate(x,&r));
105 
106   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107      Set function evaluation routine and vector
108   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109   PetscCall(SNESSetFunction(snes,r,FormFunction,(void*)&user));
110 
111   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112      Create matrix data structure; set Jacobian evaluation routine
113 
114      Set Jacobian matrix data structure and default Jacobian evaluation
115      routine. User can override with:
116      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
117                 (unless user explicitly sets preconditioner)
118      -snes_mf_operator : form preconditioning matrix as set by the user,
119                          but use matrix-free approx for Jacobian-vector
120                          products within Newton-Krylov method
121      -fdcoloring : using finite differences with coloring to compute the Jacobian
122 
123      Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option
124      below is to test the call to MatFDColoringSetType().
125      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126   PetscCall(PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL));
127   PetscCall(PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL));
128   PetscCall(PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL));
129   PetscCall(PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL));
130   if (!matrix_free) {
131     PetscCall(DMSetMatType(user.da,MATAIJ));
132     PetscCall(DMCreateMatrix(user.da,&J));
133     if (coloring) {
134       ISColoring iscoloring;
135       if (!local_coloring) {
136         PetscCall(DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring));
137         PetscCall(MatFDColoringCreate(J,iscoloring,&matfdcoloring));
138         PetscCall(MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user));
139       } else {
140         PetscCall(DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring));
141         PetscCall(MatFDColoringCreate(J,iscoloring,&matfdcoloring));
142         PetscCall(MatFDColoringUseDM(J,matfdcoloring));
143         PetscCall(MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunctionLocal,&user));
144       }
145       if (coloring_ds) {
146         PetscCall(MatFDColoringSetType(matfdcoloring,MATMFFD_DS));
147       }
148       PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
149       PetscCall(MatFDColoringSetUp(J,iscoloring,matfdcoloring));
150       PetscCall(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring));
151       PetscCall(ISColoringDestroy(&iscoloring));
152     } else {
153       PetscCall(SNESSetJacobian(snes,J,J,FormJacobian,&user));
154     }
155   }
156 
157   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158      Customize nonlinear solver; set runtime options
159    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160   PetscCall(SNESSetDM(snes,user.da));
161   PetscCall(SNESSetFromOptions(snes));
162 
163   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164      Evaluate initial guess
165      Note: The user should initialize the vector, x, with the initial guess
166      for the nonlinear solver prior to calling SNESSolve().  In particular,
167      to employ an initial guess of zero, the user should explicitly set
168      this vector to zero by calling VecSet().
169   */
170   PetscCall(FormInitialGuess(&user,x));
171 
172   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173      Solve nonlinear system
174      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175   PetscCall(SNESSolve(snes,NULL,x));
176   PetscCall(SNESGetIterationNumber(snes,&its));
177 
178   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179      Explicitly check norm of the residual of the solution
180    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181   PetscCall(FormFunction(snes,x,r,(void*)&user));
182   PetscCall(VecNorm(r,NORM_2,&fnorm));
183   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm));
184 
185   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186      Free work space.  All PETSc objects should be destroyed when they
187      are no longer needed.
188    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189 
190   PetscCall(MatDestroy(&J));
191   PetscCall(VecDestroy(&x));
192   PetscCall(VecDestroy(&r));
193   PetscCall(SNESDestroy(&snes));
194   PetscCall(DMDestroy(&user.da));
195   PetscCall(MatFDColoringDestroy(&matfdcoloring));
196   PetscCall(PetscFinalize());
197   return 0;
198 }
199 /* ------------------------------------------------------------------- */
200 /*
201    FormInitialGuess - Forms initial approximation.
202 
203    Input Parameters:
204    user - user-defined application context
205    X - vector
206 
207    Output Parameter:
208    X - vector
209  */
210 PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
211 {
212   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
213   PetscReal      lambda,temp1,hx,hy,hz,tempk,tempj;
214   PetscScalar    ***x;
215 
216   PetscFunctionBeginUser;
217   PetscCall(DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE));
218 
219   lambda = user->param;
220   hx     = 1.0/(PetscReal)(Mx-1);
221   hy     = 1.0/(PetscReal)(My-1);
222   hz     = 1.0/(PetscReal)(Mz-1);
223   temp1  = lambda/(lambda + 1.0);
224 
225   /*
226      Get a pointer to vector data.
227        - For default PETSc vectors, VecGetArray() returns a pointer to
228          the data array.  Otherwise, the routine is implementation dependent.
229        - You MUST call VecRestoreArray() when you no longer need access to
230          the array.
231   */
232   PetscCall(DMDAVecGetArray(user->da,X,&x));
233 
234   /*
235      Get local grid boundaries (for 3-dimensional DMDA):
236        xs, ys, zs   - starting grid indices (no ghost points)
237        xm, ym, zm   - widths of local grid (no ghost points)
238 
239   */
240   PetscCall(DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm));
241 
242   /*
243      Compute initial guess over the locally owned part of the grid
244   */
245   for (k=zs; k<zs+zm; k++) {
246     tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
247     for (j=ys; j<ys+ym; j++) {
248       tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
249       for (i=xs; i<xs+xm; i++) {
250         if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
251           /* boundary conditions are all zero Dirichlet */
252           x[k][j][i] = 0.0;
253         } else {
254           x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
255         }
256       }
257     }
258   }
259 
260   /*
261      Restore vector
262   */
263   PetscCall(DMDAVecRestoreArray(user->da,X,&x));
264   PetscFunctionReturn(0);
265 }
266 /* ------------------------------------------------------------------- */
267 /*
268    FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch
269 
270    Input Parameters:
271 .  snes - the SNES context
272 .  localX - input vector, this contains the ghosted region needed
273 .  ptr - optional user-defined context, as set by SNESSetFunction()
274 
275    Output Parameter:
276 .  F - function vector, this does not contain a ghosted region
277  */
278 PetscErrorCode FormFunctionLocal(SNES snes,Vec localX,Vec F,void *ptr)
279 {
280   AppCtx         *user = (AppCtx*)ptr;
281   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
282   PetscReal      two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
283   PetscScalar    u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
284   DM             da;
285 
286   PetscFunctionBeginUser;
287   PetscCall(SNESGetDM(snes,&da));
288   PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE));
289 
290   lambda  = user->param;
291   hx      = 1.0/(PetscReal)(Mx-1);
292   hy      = 1.0/(PetscReal)(My-1);
293   hz      = 1.0/(PetscReal)(Mz-1);
294   sc      = hx*hy*hz*lambda;
295   hxhzdhy = hx*hz/hy;
296   hyhzdhx = hy*hz/hx;
297   hxhydhz = hx*hy/hz;
298 
299   /*
300      Get pointers to vector data
301   */
302   PetscCall(DMDAVecGetArrayRead(da,localX,&x));
303   PetscCall(DMDAVecGetArray(da,F,&f));
304 
305   /*
306      Get local grid boundaries
307   */
308   PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm));
309 
310   /*
311      Compute function over the locally owned part of the grid
312   */
313   for (k=zs; k<zs+zm; k++) {
314     for (j=ys; j<ys+ym; j++) {
315       for (i=xs; i<xs+xm; i++) {
316         if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
317           f[k][j][i] = x[k][j][i];
318         } else {
319           u          = x[k][j][i];
320           u_east     = x[k][j][i+1];
321           u_west     = x[k][j][i-1];
322           u_north    = x[k][j+1][i];
323           u_south    = x[k][j-1][i];
324           u_up       = x[k+1][j][i];
325           u_down     = x[k-1][j][i];
326           u_xx       = (-u_east + two*u - u_west)*hyhzdhx;
327           u_yy       = (-u_north + two*u - u_south)*hxhzdhy;
328           u_zz       = (-u_up + two*u - u_down)*hxhydhz;
329           f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
330         }
331       }
332     }
333   }
334 
335   /*
336      Restore vectors
337   */
338   PetscCall(DMDAVecRestoreArrayRead(da,localX,&x));
339   PetscCall(DMDAVecRestoreArray(da,F,&f));
340   PetscCall(PetscLogFlops(11.0*ym*xm));
341   PetscFunctionReturn(0);
342 }
343 /* ------------------------------------------------------------------- */
344 /*
345    FormFunction - Evaluates nonlinear function, F(x) on the entire domain
346 
347    Input Parameters:
348 .  snes - the SNES context
349 .  X - input vector
350 .  ptr - optional user-defined context, as set by SNESSetFunction()
351 
352    Output Parameter:
353 .  F - function vector
354  */
355 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
356 {
357   Vec            localX;
358   DM             da;
359 
360   PetscFunctionBeginUser;
361   PetscCall(SNESGetDM(snes,&da));
362   PetscCall(DMGetLocalVector(da,&localX));
363 
364   /*
365      Scatter ghost points to local vector,using the 2-step process
366         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
367      By placing code between these two statements, computations can be
368      done while messages are in transition.
369   */
370   PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX));
371   PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX));
372 
373   PetscCall(FormFunctionLocal(snes,localX,F,ptr));
374   PetscCall(DMRestoreLocalVector(da,&localX));
375   PetscFunctionReturn(0);
376 }
377 /* ------------------------------------------------------------------- */
378 /*
379    FormJacobian - Evaluates Jacobian matrix.
380 
381    Input Parameters:
382 .  snes - the SNES context
383 .  x - input vector
384 .  ptr - optional user-defined context, as set by SNESSetJacobian()
385 
386    Output Parameters:
387 .  A - Jacobian matrix
388 .  B - optionally different preconditioning matrix
389 
390 */
391 PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
392 {
393   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
394   Vec            localX;
395   PetscInt       i,j,k,Mx,My,Mz;
396   MatStencil     col[7],row;
397   PetscInt       xs,ys,zs,xm,ym,zm;
398   PetscScalar    lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;
399   DM             da;
400 
401   PetscFunctionBeginUser;
402   PetscCall(SNESGetDM(snes,&da));
403   PetscCall(DMGetLocalVector(da,&localX));
404   PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE));
405 
406   lambda  = user->param;
407   hx      = 1.0/(PetscReal)(Mx-1);
408   hy      = 1.0/(PetscReal)(My-1);
409   hz      = 1.0/(PetscReal)(Mz-1);
410   sc      = hx*hy*hz*lambda;
411   hxhzdhy = hx*hz/hy;
412   hyhzdhx = hy*hz/hx;
413   hxhydhz = hx*hy/hz;
414 
415   /*
416      Scatter ghost points to local vector, using the 2-step process
417         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
418      By placing code between these two statements, computations can be
419      done while messages are in transition.
420   */
421   PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX));
422   PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX));
423 
424   /*
425      Get pointer to vector data
426   */
427   PetscCall(DMDAVecGetArrayRead(da,localX,&x));
428 
429   /*
430      Get local grid boundaries
431   */
432   PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm));
433 
434   /*
435      Compute entries for the locally owned part of the Jacobian.
436       - Currently, all PETSc parallel matrix formats are partitioned by
437         contiguous chunks of rows across the processors.
438       - Each processor needs to insert only elements that it owns
439         locally (but any non-local elements will be sent to the
440         appropriate processor during matrix assembly).
441       - Here, we set all entries for a particular row at once.
442       - We can set matrix entries either using either
443         MatSetValuesLocal() or MatSetValues(), as discussed above.
444   */
445   for (k=zs; k<zs+zm; k++) {
446     for (j=ys; j<ys+ym; j++) {
447       for (i=xs; i<xs+xm; i++) {
448         row.k = k; row.j = j; row.i = i;
449         /* boundary points */
450         if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
451           v[0] = 1.0;
452           PetscCall(MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES));
453         } else {
454           /* interior grid points */
455           v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j;  col[0].i = i;
456           v[1] = -hxhzdhy; col[1].k=k;  col[1].j=j-1;col[1].i = i;
457           v[2] = -hyhzdhx; col[2].k=k;  col[2].j=j;  col[2].i = i-1;
458           v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i;
459           v[4] = -hyhzdhx; col[4].k=k;  col[4].j=j;  col[4].i = i+1;
460           v[5] = -hxhzdhy; col[5].k=k;  col[5].j=j+1;col[5].i = i;
461           v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j;  col[6].i = i;
462           PetscCall(MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES));
463         }
464       }
465     }
466   }
467   PetscCall(DMDAVecRestoreArrayRead(da,localX,&x));
468   PetscCall(DMRestoreLocalVector(da,&localX));
469 
470   /*
471      Assemble matrix, using the 2-step process:
472        MatAssemblyBegin(), MatAssemblyEnd().
473   */
474   PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
475   PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
476 
477   /*
478      Normally since the matrix has already been assembled above; this
479      would do nothing. But in the matrix free mode -snes_mf_operator
480      this tells the "matrix-free" matrix that a new linear system solve
481      is about to be done.
482   */
483 
484   PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY));
485   PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY));
486 
487   /*
488      Tell the matrix we will never add a new nonzero location to the
489      matrix. If we do, it will generate an error.
490   */
491   PetscCall(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE));
492   PetscFunctionReturn(0);
493 }
494 
495 /*TEST
496 
497    test:
498       nsize: 4
499       args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
500 
501    test:
502       suffix: 2
503       nsize: 4
504       args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
505 
506    test:
507       suffix: 3
508       nsize: 4
509       args: -fdcoloring -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
510 
511    test:
512       suffix: 3_ds
513       nsize: 4
514       args: -fdcoloring -fdcoloring_ds -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
515 
516    test:
517       suffix: 4
518       nsize: 4
519       args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1
520       requires: !single
521 
522    test:
523       suffix: 5
524       nsize: 4
525       args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1 -snes_type newtontrdc
526       requires: !single
527 
528 TEST*/
529