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