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