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