Actual source code: ex40.c

  1: static char help[] = "Test coloring for finite difference Jacobians with DMStag\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmstag.h>
  5: #include <petscsnes.h>

  7: /*
  8:    Note that DMStagGetValuesStencil and DMStagSetValuesStencil are inefficient,
  9:    compared to DMStagVecGetArray() and friends, and only used here for testing
 10:    purposes, as they allow the code for the Jacobian and residual functions to
 11:    be more similar. In the intended application, where users are not writing
 12:    their own Jacobian assembly routines, one should use the faster, array-based
 13:    approach.
 14: */

 16: /* A "diagonal" objective function which only couples dof living at the same "point" */
 17: PetscErrorCode FormFunction1DNoCoupling(SNES snes, Vec x, Vec f, void *ctx)
 18: {
 19:   PetscInt start, n, n_extra, N, dof[2];
 20:   Vec      x_local;
 21:   DM       dm;

 23:   (void)ctx;
 24:   SNESGetDM(snes, &dm);
 25:   DMGetLocalVector(dm, &x_local);
 26:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
 27:   DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL);
 28:   DMStagGetGlobalSizes(dm, &N, NULL, NULL);
 29:   DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL);
 30:   for (PetscInt e = start; e < start + n + n_extra; ++e) {
 31:     for (PetscInt c = 0; c < dof[0]; ++c) {
 32:       DMStagStencil row;
 33:       PetscScalar   x_val, val;

 35:       row.i   = e;
 36:       row.loc = DMSTAG_LEFT;
 37:       row.c   = c;
 38:       DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
 39:       val = (10.0 + c) * x_val * x_val * x_val; // f_i = (10 +c) * x_i^3
 40:       DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
 41:     }
 42:     if (e < N) {
 43:       for (PetscInt c = 0; c < dof[1]; ++c) {
 44:         DMStagStencil row;
 45:         PetscScalar   x_val, val;

 47:         row.i   = e;
 48:         row.loc = DMSTAG_ELEMENT;
 49:         row.c   = c;
 50:         DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
 51:         val = (20.0 + c) * x_val * x_val * x_val; // f_i = (20 + c) * x_i^3
 52:         DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
 53:       }
 54:     }
 55:   }
 56:   DMRestoreLocalVector(dm, &x_local);
 57:   VecAssemblyBegin(f);
 58:   VecAssemblyEnd(f);
 59:   return 0;
 60: }

 62: PetscErrorCode FormJacobian1DNoCoupling(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx)
 63: {
 64:   PetscInt start, n, n_extra, N, dof[2];
 65:   Vec      x_local;
 66:   DM       dm;

 68:   (void)ctx;
 69:   SNESGetDM(snes, &dm);
 70:   DMGetLocalVector(dm, &x_local);
 71:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
 72:   DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL);
 73:   DMStagGetGlobalSizes(dm, &N, NULL, NULL);
 74:   DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL);
 75:   for (PetscInt e = start; e < start + n + n_extra; ++e) {
 76:     for (PetscInt c = 0; c < dof[0]; ++c) {
 77:       DMStagStencil row_vertex;
 78:       PetscScalar   x_val, val;

 80:       row_vertex.i   = e;
 81:       row_vertex.loc = DMSTAG_LEFT;
 82:       row_vertex.c   = c;
 83:       DMStagVecGetValuesStencil(dm, x_local, 1, &row_vertex, &x_val);
 84:       val = 3.0 * (10.0 + c) * x_val * x_val;
 85:       DMStagMatSetValuesStencil(dm, Amat, 1, &row_vertex, 1, &row_vertex, &val, INSERT_VALUES);
 86:     }
 87:     if (e < N) {
 88:       for (PetscInt c = 0; c < dof[1]; ++c) {
 89:         DMStagStencil row_element;
 90:         PetscScalar   x_val, val;

 92:         row_element.i   = e;
 93:         row_element.loc = DMSTAG_ELEMENT;
 94:         row_element.c   = c;
 95:         DMStagVecGetValuesStencil(dm, x_local, 1, &row_element, &x_val);
 96:         val = 3.0 * (20.0 + c) * x_val * x_val;
 97:         DMStagMatSetValuesStencil(dm, Amat, 1, &row_element, 1, &row_element, &val, INSERT_VALUES);
 98:       }
 99:     }
100:   }
101:   DMRestoreLocalVector(dm, &x_local);

103:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
104:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
106:   return 0;
107: }

109: /* Objective functions which use the DM's stencil width. */
110: PetscErrorCode FormFunction1D(SNES snes, Vec x, Vec f, void *ctx)
111: {
112:   Vec               x_local;
113:   PetscInt          dim, stencil_width, start, n, n_extra, N, dof[2];
114:   DMStagStencilType stencil_type;
115:   DM                dm;

117:   (void)ctx;
118:   SNESGetDM(snes, &dm);
119:   DMGetDimension(dm, &dim);
121:   DMStagGetStencilType(dm, &stencil_type);
123:   DMStagGetStencilWidth(dm, &stencil_width);

125:   DMGetLocalVector(dm, &x_local);
126:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
127:   DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL);
128:   DMStagGetGlobalSizes(dm, &N, NULL, NULL);
129:   DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL);

131:   VecZeroEntries(f);

133:   for (PetscInt e = start; e < start + n + n_extra; ++e) {
134:     DMStagStencil row_vertex, row_element;

136:     row_vertex.i   = e;
137:     row_vertex.loc = DMSTAG_LEFT;

139:     row_element.i   = e;
140:     row_element.loc = DMSTAG_ELEMENT;

142:     for (PetscInt offset = -stencil_width; offset <= stencil_width; ++offset) {
143:       const PetscInt e_offset = e + offset;

145:       // vertex --> vertex
146:       if (e_offset >= 0 && e_offset < N + 1) { // Does not fully wrap in the periodic case
147:         DMStagStencil col;
148:         PetscScalar   x_val, val;

150:         for (PetscInt c_row = 0; c_row < dof[0]; ++c_row) {
151:           row_vertex.c = c_row;
152:           for (PetscInt c_col = 0; c_col < dof[0]; ++c_col) {
153:             col.c   = c_col;
154:             col.i   = e_offset;
155:             col.loc = DMSTAG_LEFT;
156:             DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
157:             val = (10.0 + offset) * x_val * x_val * x_val;
158:             DMStagVecSetValuesStencil(dm, f, 1, &row_vertex, &val, ADD_VALUES);
159:           }
160:         }
161:       }

163:       // element --> vertex
164:       if (e_offset >= 0 && e_offset < N) { // Does not fully wrap in the periodic case
165:         DMStagStencil col;
166:         PetscScalar   x_val, val;

168:         for (PetscInt c_row = 0; c_row < dof[0]; ++c_row) {
169:           row_vertex.c = c_row;
170:           for (PetscInt c_col = 0; c_col < dof[1]; ++c_col) {
171:             col.c   = c_col;
172:             col.i   = e_offset;
173:             col.loc = DMSTAG_ELEMENT;
174:             DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
175:             val = (15.0 + offset) * x_val * x_val * x_val;
176:             DMStagVecSetValuesStencil(dm, f, 1, &row_vertex, &val, ADD_VALUES);
177:           }
178:         }
179:       }

181:       if (e < N) {
182:         // vertex --> element
183:         if (e_offset >= 0 && e_offset < N + 1) { // Does not fully wrap in the periodic case
184:           DMStagStencil col;
185:           PetscScalar   x_val, val;

187:           for (PetscInt c_row = 0; c_row < dof[1]; ++c_row) {
188:             row_element.c = c_row;
189:             for (PetscInt c_col = 0; c_col < dof[0]; ++c_col) {
190:               col.c   = c_col;
191:               col.i   = e_offset;
192:               col.loc = DMSTAG_LEFT;
193:               DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
194:               val = (25.0 + offset) * x_val * x_val * x_val;
195:               DMStagVecSetValuesStencil(dm, f, 1, &row_element, &val, ADD_VALUES);
196:             }
197:           }
198:         }

200:         // element --> element
201:         if (e_offset >= 0 && e_offset < N) { // Does not fully wrap in the periodic case
202:           DMStagStencil col;
203:           PetscScalar   x_val, val;

205:           for (PetscInt c_row = 0; c_row < dof[1]; ++c_row) {
206:             row_element.c = c_row;
207:             for (PetscInt c_col = 0; c_col < dof[1]; ++c_col) {
208:               col.c   = c_col;
209:               col.i   = e_offset;
210:               col.loc = DMSTAG_ELEMENT;
211:               DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
212:               val = (20.0 + offset) * x_val * x_val * x_val;
213:               DMStagVecSetValuesStencil(dm, f, 1, &row_element, &val, ADD_VALUES);
214:             }
215:           }
216:         }
217:       }
218:     }
219:   }
220:   DMRestoreLocalVector(dm, &x_local);
221:   VecAssemblyBegin(f);
222:   VecAssemblyEnd(f);
223:   return 0;
224: }

226: PetscErrorCode FormJacobian1D(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx)
227: {
228:   Vec      x_local;
229:   PetscInt dim, stencil_width, start, n, n_extra, N, dof[2];
230:   DM       dm;

232:   (void)ctx;
233:   SNESGetDM(snes, &dm);
234:   DMGetDimension(dm, &dim);
236:   DMStagGetStencilWidth(dm, &stencil_width);

238:   DMGetLocalVector(dm, &x_local);
239:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
240:   DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL);
241:   DMStagGetGlobalSizes(dm, &N, NULL, NULL);
242:   DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL);

244:   MatZeroEntries(Amat);

246:   for (PetscInt e = start; e < start + n + n_extra; ++e) {
247:     DMStagStencil row_vertex, row_element;

249:     row_vertex.i   = e;
250:     row_vertex.loc = DMSTAG_LEFT;

252:     row_element.i   = e;
253:     row_element.loc = DMSTAG_ELEMENT;

255:     for (PetscInt offset = -stencil_width; offset <= stencil_width; ++offset) {
256:       const PetscInt e_offset = e + offset;

258:       // vertex --> vertex
259:       if (e_offset >= 0 && e_offset < N + 1) {
260:         DMStagStencil col;
261:         PetscScalar   x_val, val;

263:         for (PetscInt c_row = 0; c_row < dof[0]; ++c_row) {
264:           row_vertex.c = c_row;
265:           for (PetscInt c_col = 0; c_col < dof[0]; ++c_col) {
266:             col.c   = c_col;
267:             col.i   = e_offset;
268:             col.loc = DMSTAG_LEFT;
269:             DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
270:             val = 3.0 * (10.0 + offset) * x_val * x_val;
271:             DMStagMatSetValuesStencil(dm, Amat, 1, &row_vertex, 1, &col, &val, ADD_VALUES);
272:           }
273:         }
274:       }

276:       // element --> vertex
277:       if (e_offset >= 0 && e_offset < N) {
278:         DMStagStencil col;
279:         PetscScalar   x_val, val;

281:         for (PetscInt c_row = 0; c_row < dof[0]; ++c_row) {
282:           row_vertex.c = c_row;
283:           for (PetscInt c_col = 0; c_col < dof[1]; ++c_col) {
284:             col.c   = c_col;
285:             col.i   = e_offset;
286:             col.loc = DMSTAG_ELEMENT;
287:             DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
288:             val = 3.0 * (15.0 + offset) * x_val * x_val;
289:             DMStagMatSetValuesStencil(dm, Amat, 1, &row_vertex, 1, &col, &val, ADD_VALUES);
290:           }
291:         }
292:       }

294:       if (e < N) {
295:         // vertex --> element
296:         if (e_offset >= 0 && e_offset < N + 1) {
297:           DMStagStencil col;
298:           PetscScalar   x_val, val;

300:           for (PetscInt c_row = 0; c_row < dof[1]; ++c_row) {
301:             row_element.c = c_row;
302:             for (PetscInt c_col = 0; c_col < dof[0]; ++c_col) {
303:               col.c   = c_col;
304:               col.i   = e_offset;
305:               col.loc = DMSTAG_LEFT;
306:               DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
307:               val = 3.0 * (25.0 + offset) * x_val * x_val;
308:               DMStagMatSetValuesStencil(dm, Amat, 1, &row_element, 1, &col, &val, ADD_VALUES);
309:             }
310:           }
311:         }

313:         // element --> element
314:         if (e_offset >= 0 && e_offset < N) {
315:           DMStagStencil col;
316:           PetscScalar   x_val, val;

318:           for (PetscInt c_row = 0; c_row < dof[1]; ++c_row) {
319:             row_element.c = c_row;
320:             for (PetscInt c_col = 0; c_col < dof[1]; ++c_col) {
321:               col.c   = c_col;
322:               col.i   = e_offset;
323:               col.loc = DMSTAG_ELEMENT;
324:               DMStagVecGetValuesStencil(dm, x_local, 1, &col, &x_val);
325:               val = 3.0 * (20.0 + offset) * x_val * x_val;
326:               DMStagMatSetValuesStencil(dm, Amat, 1, &row_element, 1, &col, &val, ADD_VALUES);
327:             }
328:           }
329:         }
330:       }
331:     }
332:   }
333:   DMRestoreLocalVector(dm, &x_local);
334:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
335:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
337:   return 0;
338: }

340: PetscErrorCode FormFunction2DNoCoupling(SNES snes, Vec x, Vec f, void *ctx)
341: {
342:   PetscInt start[2], n[2], n_extra[2], N[2], dof[3];
343:   Vec      x_local;
344:   DM       dm;

346:   (void)ctx;
347:   SNESGetDM(snes, &dm);
348:   DMGetLocalVector(dm, &x_local);
349:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
350:   DMStagGetCorners(dm, &start[0], &start[1], NULL, &n[0], &n[1], NULL, &n_extra[0], &n_extra[1], NULL);
351:   DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL);
352:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL);
353:   for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
354:     for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
355:       for (PetscInt c = 0; c < dof[0]; ++c) {
356:         DMStagStencil row;
357:         PetscScalar   x_val, val;

359:         row.i   = ex;
360:         row.j   = ey;
361:         row.loc = DMSTAG_DOWN_LEFT;
362:         row.c   = c;
363:         DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
364:         val = (5.0 + c) * x_val * x_val * x_val;
365:         DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
366:       }
367:       if (ex < N[0]) {
368:         for (PetscInt c = 0; c < dof[1]; ++c) {
369:           DMStagStencil row;
370:           PetscScalar   x_val, val;

372:           row.i   = ex;
373:           row.j   = ey;
374:           row.loc = DMSTAG_DOWN;
375:           row.c   = c;
376:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
377:           val = (10.0 + c) * x_val * x_val * x_val;
378:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
379:         }
380:       }
381:       if (ey < N[1]) {
382:         for (PetscInt c = 0; c < dof[1]; ++c) {
383:           DMStagStencil row;
384:           PetscScalar   x_val, val;

386:           row.i   = ex;
387:           row.j   = ey;
388:           row.loc = DMSTAG_LEFT;
389:           row.c   = c;
390:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
391:           val = (15.0 + c) * x_val * x_val * x_val;
392:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
393:         }
394:       }
395:       if (ex < N[0] && ey < N[1]) {
396:         for (PetscInt c = 0; c < dof[2]; ++c) {
397:           DMStagStencil row;
398:           PetscScalar   x_val, val;

400:           row.i   = ex;
401:           row.j   = ey;
402:           row.loc = DMSTAG_ELEMENT;
403:           row.c   = c;
404:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
405:           val = (20.0 + c) * x_val * x_val * x_val;
406:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
407:         }
408:       }
409:     }
410:   }
411:   DMRestoreLocalVector(dm, &x_local);
412:   VecAssemblyBegin(f);
413:   VecAssemblyEnd(f);
414:   return 0;
415: }

417: PetscErrorCode FormJacobian2DNoCoupling(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx)
418: {
419:   PetscInt start[2], n[2], n_extra[2], N[2], dof[3];
420:   Vec      x_local;
421:   DM       dm;

423:   (void)ctx;
424:   SNESGetDM(snes, &dm);
425:   DMGetLocalVector(dm, &x_local);
426:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
427:   DMStagGetCorners(dm, &start[0], &start[1], NULL, &n[0], &n[1], NULL, &n_extra[0], &n_extra[1], NULL);
428:   DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL);
429:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL);
430:   for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
431:     for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
432:       for (PetscInt c = 0; c < dof[0]; ++c) {
433:         DMStagStencil row;
434:         PetscScalar   x_val, val;

436:         row.i   = ex;
437:         row.j   = ey;
438:         row.loc = DMSTAG_DOWN_LEFT;
439:         row.c   = c;
440:         DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
441:         val = 3.0 * (5.0 + c) * x_val * x_val;
442:         DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
443:       }
444:       if (ex < N[0]) {
445:         for (PetscInt c = 0; c < dof[1]; ++c) {
446:           DMStagStencil row;
447:           PetscScalar   x_val, val;

449:           row.i   = ex;
450:           row.j   = ey;
451:           row.loc = DMSTAG_DOWN;
452:           row.c   = c;
453:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
454:           val = 3.0 * (10.0 + c) * x_val * x_val;
455:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
456:         }
457:       }
458:       if (ey < N[1]) {
459:         for (PetscInt c = 0; c < dof[1]; ++c) {
460:           DMStagStencil row;
461:           PetscScalar   x_val, val;

463:           row.i   = ex;
464:           row.j   = ey;
465:           row.loc = DMSTAG_LEFT;
466:           row.c   = c;
467:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
468:           val = 3.0 * (15.0 + c) * x_val * x_val;
469:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
470:         }
471:       }
472:       if (ex < N[0] && ey < N[1]) {
473:         for (PetscInt c = 0; c < dof[2]; ++c) {
474:           DMStagStencil row;
475:           PetscScalar   x_val, val;

477:           row.i   = ex;
478:           row.j   = ey;
479:           row.loc = DMSTAG_ELEMENT;
480:           row.c   = c;
481:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
482:           val = 3.0 * (20.0 + c) * x_val * x_val;
483:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
484:         }
485:       }
486:     }
487:   }
488:   DMRestoreLocalVector(dm, &x_local);

490:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
491:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
493:   return 0;
494: }

496: PetscErrorCode FormFunction2D(SNES snes, Vec x, Vec f, void *ctx)
497: {
498:   PetscInt start[2], n[2], n_extra[2], N[2], dof[3];
499:   Vec      x_local;
500:   DM       dm;

502:   (void)ctx;
503:   SNESGetDM(snes, &dm);
504:   DMGetLocalVector(dm, &x_local);
505:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
506:   DMStagGetCorners(dm, &start[0], &start[1], NULL, &n[0], &n[1], NULL, &n_extra[0], &n_extra[1], NULL);
507:   DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL);
508:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL);

510:   VecZeroEntries(f);

512:   /* First, as in the simple case above */
513:   for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
514:     for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
515:       for (PetscInt c = 0; c < dof[0]; ++c) {
516:         DMStagStencil row;
517:         PetscScalar   x_val, val;

519:         row.i   = ex;
520:         row.j   = ey;
521:         row.loc = DMSTAG_DOWN_LEFT;
522:         row.c   = c;
523:         DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
524:         val = (5.0 + c) * x_val * x_val * x_val;
525:         DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
526:       }
527:       if (ex < N[0]) {
528:         for (PetscInt c = 0; c < dof[1]; ++c) {
529:           DMStagStencil row;
530:           PetscScalar   x_val, val;

532:           row.i   = ex;
533:           row.j   = ey;
534:           row.loc = DMSTAG_DOWN;
535:           row.c   = c;
536:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
537:           val = (10.0 + c) * x_val * x_val * x_val;
538:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
539:         }
540:       }
541:       if (ey < N[1]) {
542:         for (PetscInt c = 0; c < dof[1]; ++c) {
543:           DMStagStencil row;
544:           PetscScalar   x_val, val;

546:           row.i   = ex;
547:           row.j   = ey;
548:           row.loc = DMSTAG_LEFT;
549:           row.c   = c;
550:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
551:           val = (15.0 + c) * x_val * x_val * x_val;
552:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
553:         }
554:       }
555:       if (ex < N[0] && ey < N[1]) {
556:         for (PetscInt c = 0; c < dof[2]; ++c) {
557:           DMStagStencil row;
558:           PetscScalar   x_val, val;

560:           row.i   = ex;
561:           row.j   = ey;
562:           row.loc = DMSTAG_ELEMENT;
563:           row.c   = c;
564:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
565:           val = (20.0 + c) * x_val * x_val * x_val;
566:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
567:         }
568:       }
569:     }
570:   }

572:   /* Add additional terms fully coupling one interior element to another */
573:   {
574:     PetscMPIInt rank;

576:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
577:     if (rank == 0) {
578:       PetscInt       epe;
579:       DMStagStencil *row, *col;

581:       DMStagGetEntriesPerElement(dm, &epe);
582:       PetscMalloc1(epe, &row);
583:       PetscMalloc1(epe, &col);
584:       for (PetscInt i = 0; i < epe; ++i) {
585:         row[i].i = 0;
586:         row[i].j = 0;
587:         col[i].i = 0;
588:         col[i].j = 1;
589:       }
590:       {
591:         PetscInt nrows = 0;

593:         for (PetscInt c = 0; c < dof[0]; ++c) {
594:           row[nrows].c   = c;
595:           row[nrows].loc = DMSTAG_DOWN_LEFT;
596:           ++nrows;
597:         }
598:         for (PetscInt c = 0; c < dof[1]; ++c) {
599:           row[nrows].c   = c;
600:           row[nrows].loc = DMSTAG_LEFT;
601:           ++nrows;
602:         }
603:         for (PetscInt c = 0; c < dof[1]; ++c) {
604:           row[nrows].c   = c;
605:           row[nrows].loc = DMSTAG_DOWN;
606:           ++nrows;
607:         }
608:         for (PetscInt c = 0; c < dof[2]; ++c) {
609:           row[nrows].c   = c;
610:           row[nrows].loc = DMSTAG_ELEMENT;
611:           ++nrows;
612:         }
613:       }

615:       {
616:         PetscInt ncols = 0;

618:         for (PetscInt c = 0; c < dof[0]; ++c) {
619:           col[ncols].c   = c;
620:           col[ncols].loc = DMSTAG_DOWN_LEFT;
621:           ++ncols;
622:         }
623:         for (PetscInt c = 0; c < dof[1]; ++c) {
624:           col[ncols].c   = c;
625:           col[ncols].loc = DMSTAG_LEFT;
626:           ++ncols;
627:         }
628:         for (PetscInt c = 0; c < dof[1]; ++c) {
629:           col[ncols].c   = c;
630:           col[ncols].loc = DMSTAG_DOWN;
631:           ++ncols;
632:         }
633:         for (PetscInt c = 0; c < dof[2]; ++c) {
634:           col[ncols].c   = c;
635:           col[ncols].loc = DMSTAG_ELEMENT;
636:           ++ncols;
637:         }
638:       }

640:       for (PetscInt i = 0; i < epe; ++i) {
641:         for (PetscInt j = 0; j < epe; ++j) {
642:           PetscScalar x_val, val;

644:           DMStagVecGetValuesStencil(dm, x_local, 1, &col[j], &x_val);
645:           val = (10 * i + j) * x_val * x_val * x_val;
646:           DMStagVecSetValuesStencil(dm, f, 1, &row[i], &val, ADD_VALUES);
647:         }
648:       }
649:       PetscFree(row);
650:       PetscFree(col);
651:     }
652:   }
653:   DMRestoreLocalVector(dm, &x_local);
654:   VecAssemblyBegin(f);
655:   VecAssemblyEnd(f);
656:   return 0;
657: }

659: PetscErrorCode FormJacobian2D(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx)
660: {
661:   PetscInt start[2], n[2], n_extra[2], N[2], dof[3];
662:   Vec      x_local;
663:   DM       dm;

665:   (void)ctx;
666:   SNESGetDM(snes, &dm);
667:   DMGetLocalVector(dm, &x_local);
668:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
669:   DMStagGetCorners(dm, &start[0], &start[1], NULL, &n[0], &n[1], NULL, &n_extra[0], &n_extra[1], NULL);
670:   DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL);
671:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL);
672:   MatZeroEntries(Amat);
673:   for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
674:     for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
675:       for (PetscInt c = 0; c < dof[0]; ++c) {
676:         DMStagStencil row;
677:         PetscScalar   x_val, val;

679:         row.i   = ex;
680:         row.j   = ey;
681:         row.loc = DMSTAG_DOWN_LEFT;
682:         row.c   = c;
683:         DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
684:         val = 3.0 * (5.0 + c) * x_val * x_val;
685:         DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
686:       }
687:       if (ex < N[0]) {
688:         for (PetscInt c = 0; c < dof[1]; ++c) {
689:           DMStagStencil row;
690:           PetscScalar   x_val, val;

692:           row.i   = ex;
693:           row.j   = ey;
694:           row.loc = DMSTAG_DOWN;
695:           row.c   = c;
696:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
697:           val = 3.0 * (10.0 + c) * x_val * x_val;
698:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
699:         }
700:       }
701:       if (ey < N[1]) {
702:         for (PetscInt c = 0; c < dof[1]; ++c) {
703:           DMStagStencil row;
704:           PetscScalar   x_val, val;

706:           row.i   = ex;
707:           row.j   = ey;
708:           row.loc = DMSTAG_LEFT;
709:           row.c   = c;
710:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
711:           val = 3.0 * (15.0 + c) * x_val * x_val;
712:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
713:         }
714:       }
715:       if (ex < N[0] && ey < N[1]) {
716:         for (PetscInt c = 0; c < dof[2]; ++c) {
717:           DMStagStencil row;
718:           PetscScalar   x_val, val;

720:           row.i   = ex;
721:           row.j   = ey;
722:           row.loc = DMSTAG_ELEMENT;
723:           row.c   = c;
724:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
725:           val = 3.0 * (20.0 + c) * x_val * x_val;
726:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
727:         }
728:       }
729:     }
730:   }

732:   /* Add additional terms fully coupling one interior element to another */
733:   {
734:     PetscMPIInt rank;

736:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
737:     if (rank == 0) {
738:       PetscInt       epe;
739:       DMStagStencil *row, *col;

741:       DMStagGetEntriesPerElement(dm, &epe);
742:       PetscMalloc1(epe, &row);
743:       PetscMalloc1(epe, &col);
744:       for (PetscInt i = 0; i < epe; ++i) {
745:         row[i].i = 0;
746:         row[i].j = 0;
747:         col[i].i = 0;
748:         col[i].j = 1;
749:       }
750:       {
751:         PetscInt nrows = 0;

753:         for (PetscInt c = 0; c < dof[0]; ++c) {
754:           row[nrows].c   = c;
755:           row[nrows].loc = DMSTAG_DOWN_LEFT;
756:           ++nrows;
757:         }
758:         for (PetscInt c = 0; c < dof[1]; ++c) {
759:           row[nrows].c   = c;
760:           row[nrows].loc = DMSTAG_LEFT;
761:           ++nrows;
762:         }
763:         for (PetscInt c = 0; c < dof[1]; ++c) {
764:           row[nrows].c   = c;
765:           row[nrows].loc = DMSTAG_DOWN;
766:           ++nrows;
767:         }
768:         for (PetscInt c = 0; c < dof[2]; ++c) {
769:           row[nrows].c   = c;
770:           row[nrows].loc = DMSTAG_ELEMENT;
771:           ++nrows;
772:         }
773:       }

775:       {
776:         PetscInt ncols = 0;

778:         for (PetscInt c = 0; c < dof[0]; ++c) {
779:           col[ncols].c   = c;
780:           col[ncols].loc = DMSTAG_DOWN_LEFT;
781:           ++ncols;
782:         }
783:         for (PetscInt c = 0; c < dof[1]; ++c) {
784:           col[ncols].c   = c;
785:           col[ncols].loc = DMSTAG_LEFT;
786:           ++ncols;
787:         }
788:         for (PetscInt c = 0; c < dof[1]; ++c) {
789:           col[ncols].c   = c;
790:           col[ncols].loc = DMSTAG_DOWN;
791:           ++ncols;
792:         }
793:         for (PetscInt c = 0; c < dof[2]; ++c) {
794:           col[ncols].c   = c;
795:           col[ncols].loc = DMSTAG_ELEMENT;
796:           ++ncols;
797:         }
798:       }

800:       for (PetscInt i = 0; i < epe; ++i) {
801:         for (PetscInt j = 0; j < epe; ++j) {
802:           PetscScalar x_val, val;

804:           DMStagVecGetValuesStencil(dm, x_local, 1, &col[j], &x_val);
805:           val = 3.0 * (10 * i + j) * x_val * x_val;
806:           DMStagMatSetValuesStencil(dm, Amat, 1, &row[i], 1, &col[j], &val, ADD_VALUES);
807:         }
808:       }
809:       PetscFree(row);
810:       PetscFree(col);
811:     }
812:   }
813:   DMRestoreLocalVector(dm, &x_local);
814:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
815:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
817:   return 0;
818: }

820: PetscErrorCode FormFunction3DNoCoupling(SNES snes, Vec x, Vec f, void *ctx)
821: {
822:   PetscInt start[3], n[3], n_extra[3], N[3], dof[4];
823:   Vec      x_local;
824:   DM       dm;

826:   (void)ctx;
827:   SNESGetDM(snes, &dm);
828:   DMGetLocalVector(dm, &x_local);
829:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
830:   DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &n_extra[0], &n_extra[1], &n_extra[2]);
831:   DMStagGetGlobalSizes(dm, &N[0], &N[1], &N[2]);
832:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], &dof[3]);
833:   for (PetscInt ez = start[2]; ez < start[2] + n[2] + n_extra[2]; ++ez) {
834:     for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
835:       for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
836:         for (PetscInt c = 0; c < dof[0]; ++c) {
837:           DMStagStencil row;
838:           PetscScalar   x_val, val;

840:           row.i   = ex;
841:           row.j   = ey;
842:           row.k   = ez;
843:           row.loc = DMSTAG_BACK_DOWN_LEFT;
844:           row.c   = c;
845:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
846:           val = (5.0 + c) * x_val * x_val * x_val;
847:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
848:         }
849:         if (ez < N[2]) {
850:           for (PetscInt c = 0; c < dof[1]; ++c) {
851:             DMStagStencil row;
852:             PetscScalar   x_val, val;

854:             row.i   = ex;
855:             row.j   = ey;
856:             row.k   = ez;
857:             row.loc = DMSTAG_DOWN_LEFT;
858:             row.c   = c;
859:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
860:             val = (50.0 + c) * x_val * x_val * x_val;
861:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
862:           }
863:         }
864:         if (ey < N[1]) {
865:           for (PetscInt c = 0; c < dof[1]; ++c) {
866:             DMStagStencil row;
867:             PetscScalar   x_val, val;

869:             row.i   = ex;
870:             row.j   = ey;
871:             row.k   = ez;
872:             row.loc = DMSTAG_BACK_LEFT;
873:             row.c   = c;
874:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
875:             val = (55.0 + c) * x_val * x_val * x_val;
876:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
877:           }
878:         }
879:         if (ex < N[0]) {
880:           for (PetscInt c = 0; c < dof[1]; ++c) {
881:             DMStagStencil row;
882:             PetscScalar   x_val, val;

884:             row.i   = ex;
885:             row.j   = ey;
886:             row.k   = ez;
887:             row.loc = DMSTAG_BACK_DOWN;
888:             row.c   = c;
889:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
890:             val = (60.0 + c) * x_val * x_val * x_val;
891:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
892:           }
893:         }
894:         if (ex < N[0] && ez < N[2]) {
895:           for (PetscInt c = 0; c < dof[2]; ++c) {
896:             DMStagStencil row;
897:             PetscScalar   x_val, val;

899:             row.i   = ex;
900:             row.j   = ey;
901:             row.k   = ez;
902:             row.loc = DMSTAG_DOWN;
903:             row.c   = c;
904:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
905:             val = (10.0 + c) * x_val * x_val * x_val;
906:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
907:           }
908:         }
909:         if (ey < N[1] && ez < N[2]) {
910:           for (PetscInt c = 0; c < dof[2]; ++c) {
911:             DMStagStencil row;
912:             PetscScalar   x_val, val;

914:             row.i   = ex;
915:             row.j   = ey;
916:             row.k   = ez;
917:             row.loc = DMSTAG_LEFT;
918:             row.c   = c;
919:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
920:             val = (15.0 + c) * x_val * x_val * x_val;
921:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
922:           }
923:         }
924:         if (ex < N[0] && ey < N[1]) {
925:           for (PetscInt c = 0; c < dof[2]; ++c) {
926:             DMStagStencil row;
927:             PetscScalar   x_val, val;

929:             row.i   = ex;
930:             row.j   = ey;
931:             row.k   = ez;
932:             row.loc = DMSTAG_BACK;
933:             row.c   = c;
934:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
935:             val = (15.0 + c) * x_val * x_val * x_val;
936:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
937:           }
938:         }
939:         if (ex < N[0] && ey < N[1] && ez < N[2]) {
940:           for (PetscInt c = 0; c < dof[3]; ++c) {
941:             DMStagStencil row;
942:             PetscScalar   x_val, val;

944:             row.i   = ex;
945:             row.j   = ey;
946:             row.k   = ez;
947:             row.loc = DMSTAG_ELEMENT;
948:             row.c   = c;
949:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
950:             val = (20.0 + c) * x_val * x_val * x_val;
951:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, INSERT_VALUES);
952:           }
953:         }
954:       }
955:     }
956:   }
957:   DMRestoreLocalVector(dm, &x_local);
958:   VecAssemblyBegin(f);
959:   VecAssemblyEnd(f);
960:   return 0;
961: }

963: PetscErrorCode FormJacobian3DNoCoupling(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx)
964: {
965:   PetscInt start[3], n[3], n_extra[3], N[3], dof[4];
966:   Vec      x_local;
967:   DM       dm;

969:   (void)ctx;
970:   SNESGetDM(snes, &dm);
971:   DMGetLocalVector(dm, &x_local);
972:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
973:   DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &n_extra[0], &n_extra[1], &n_extra[2]);
974:   DMStagGetGlobalSizes(dm, &N[0], &N[1], &N[2]);
975:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], &dof[3]);
976:   for (PetscInt ez = start[2]; ez < start[2] + n[2] + n_extra[2]; ++ez) {
977:     for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
978:       for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
979:         for (PetscInt c = 0; c < dof[0]; ++c) {
980:           DMStagStencil row;
981:           PetscScalar   x_val, val;

983:           row.i   = ex;
984:           row.j   = ey;
985:           row.k   = ez;
986:           row.loc = DMSTAG_BACK_DOWN_LEFT;
987:           row.c   = c;
988:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
989:           val = 3.0 * (5.0 + c) * x_val * x_val;
990:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
991:         }
992:         if (ez < N[2]) {
993:           for (PetscInt c = 0; c < dof[1]; ++c) {
994:             DMStagStencil row;
995:             PetscScalar   x_val, val;

997:             row.i   = ex;
998:             row.j   = ey;
999:             row.k   = ez;
1000:             row.loc = DMSTAG_DOWN_LEFT;
1001:             row.c   = c;
1002:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1003:             val = 3.0 * (50.0 + c) * x_val * x_val;
1004:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1005:           }
1006:         }
1007:         if (ey < N[1]) {
1008:           for (PetscInt c = 0; c < dof[1]; ++c) {
1009:             DMStagStencil row;
1010:             PetscScalar   x_val, val;

1012:             row.i   = ex;
1013:             row.j   = ey;
1014:             row.k   = ez;
1015:             row.loc = DMSTAG_BACK_LEFT;
1016:             row.c   = c;
1017:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1018:             val = 3.0 * (55.0 + c) * x_val * x_val;
1019:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1020:           }
1021:         }
1022:         if (ex < N[0]) {
1023:           for (PetscInt c = 0; c < dof[1]; ++c) {
1024:             DMStagStencil row;
1025:             PetscScalar   x_val, val;

1027:             row.i   = ex;
1028:             row.j   = ey;
1029:             row.k   = ez;
1030:             row.loc = DMSTAG_BACK_DOWN;
1031:             row.c   = c;
1032:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1033:             val = 3.0 * (60.0 + c) * x_val * x_val;
1034:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1035:           }
1036:         }
1037:         if (ex < N[0] && ez < N[2]) {
1038:           for (PetscInt c = 0; c < dof[2]; ++c) {
1039:             DMStagStencil row;
1040:             PetscScalar   x_val, val;

1042:             row.i   = ex;
1043:             row.j   = ey;
1044:             row.k   = ez;
1045:             row.loc = DMSTAG_DOWN;
1046:             row.c   = c;
1047:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1048:             val = 3.0 * (10.0 + c) * x_val * x_val;
1049:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1050:           }
1051:         }
1052:         if (ey < N[1] && ez < N[2]) {
1053:           for (PetscInt c = 0; c < dof[2]; ++c) {
1054:             DMStagStencil row;
1055:             PetscScalar   x_val, val;

1057:             row.i   = ex;
1058:             row.j   = ey;
1059:             row.k   = ez;
1060:             row.loc = DMSTAG_LEFT;
1061:             row.c   = c;
1062:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1063:             val = 3.0 * (15.0 + c) * x_val * x_val;
1064:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1065:           }
1066:         }
1067:         if (ex < N[0] && ey < N[1]) {
1068:           for (PetscInt c = 0; c < dof[2]; ++c) {
1069:             DMStagStencil row;
1070:             PetscScalar   x_val, val;

1072:             row.i   = ex;
1073:             row.j   = ey;
1074:             row.k   = ez;
1075:             row.loc = DMSTAG_BACK;
1076:             row.c   = c;
1077:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1078:             val = 3.0 * (15.0 + c) * x_val * x_val;
1079:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1080:           }
1081:         }
1082:         if (ex < N[0] && ey < N[1] && ez < N[2]) {
1083:           for (PetscInt c = 0; c < dof[3]; ++c) {
1084:             DMStagStencil row;
1085:             PetscScalar   x_val, val;

1087:             row.i   = ex;
1088:             row.j   = ey;
1089:             row.k   = ez;
1090:             row.loc = DMSTAG_ELEMENT;
1091:             row.c   = c;
1092:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1093:             val = 3.0 * (20.0 + c) * x_val * x_val;
1094:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, INSERT_VALUES);
1095:           }
1096:         }
1097:       }
1098:     }
1099:   }
1100:   DMRestoreLocalVector(dm, &x_local);
1101:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
1102:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
1104:   return 0;
1105: }

1107: PetscErrorCode FormFunction3D(SNES snes, Vec x, Vec f, void *ctx)
1108: {
1109:   PetscInt start[3], n[3], n_extra[3], N[3], dof[4];
1110:   Vec      x_local;
1111:   DM       dm;

1113:   (void)ctx;
1114:   SNESGetDM(snes, &dm);
1115:   DMGetLocalVector(dm, &x_local);
1116:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
1117:   DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &n_extra[0], &n_extra[1], &n_extra[2]);
1118:   DMStagGetGlobalSizes(dm, &N[0], &N[1], &N[2]);
1119:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], &dof[3]);
1120:   VecZeroEntries(f);
1121:   for (PetscInt ez = start[2]; ez < start[2] + n[2] + n_extra[2]; ++ez) {
1122:     for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
1123:       for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
1124:         for (PetscInt c = 0; c < dof[0]; ++c) {
1125:           DMStagStencil row;
1126:           PetscScalar   x_val, val;

1128:           row.i   = ex;
1129:           row.j   = ey;
1130:           row.k   = ez;
1131:           row.loc = DMSTAG_BACK_DOWN_LEFT;
1132:           row.c   = c;
1133:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1134:           val = (5.0 + c) * x_val * x_val * x_val;
1135:           DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1136:         }
1137:         if (ez < N[2]) {
1138:           for (PetscInt c = 0; c < dof[1]; ++c) {
1139:             DMStagStencil row;
1140:             PetscScalar   x_val, val;

1142:             row.i   = ex;
1143:             row.j   = ey;
1144:             row.k   = ez;
1145:             row.loc = DMSTAG_DOWN_LEFT;
1146:             row.c   = c;
1147:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1148:             val = (50.0 + c) * x_val * x_val * x_val;
1149:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1150:           }
1151:         }
1152:         if (ey < N[1]) {
1153:           for (PetscInt c = 0; c < dof[1]; ++c) {
1154:             DMStagStencil row;
1155:             PetscScalar   x_val, val;

1157:             row.i   = ex;
1158:             row.j   = ey;
1159:             row.k   = ez;
1160:             row.loc = DMSTAG_BACK_LEFT;
1161:             row.c   = c;
1162:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1163:             val = (55.0 + c) * x_val * x_val * x_val;
1164:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1165:           }
1166:         }
1167:         if (ex < N[0]) {
1168:           for (PetscInt c = 0; c < dof[1]; ++c) {
1169:             DMStagStencil row;
1170:             PetscScalar   x_val, val;

1172:             row.i   = ex;
1173:             row.j   = ey;
1174:             row.k   = ez;
1175:             row.loc = DMSTAG_BACK_DOWN;
1176:             row.c   = c;
1177:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1178:             val = (60.0 + c) * x_val * x_val * x_val;
1179:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1180:           }
1181:         }
1182:         if (ex < N[0] && ez < N[2]) {
1183:           for (PetscInt c = 0; c < dof[2]; ++c) {
1184:             DMStagStencil row;
1185:             PetscScalar   x_val, val;

1187:             row.i   = ex;
1188:             row.j   = ey;
1189:             row.k   = ez;
1190:             row.loc = DMSTAG_DOWN;
1191:             row.c   = c;
1192:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1193:             val = (10.0 + c) * x_val * x_val * x_val;
1194:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1195:           }
1196:         }
1197:         if (ey < N[1] && ez < N[2]) {
1198:           for (PetscInt c = 0; c < dof[2]; ++c) {
1199:             DMStagStencil row;
1200:             PetscScalar   x_val, val;

1202:             row.i   = ex;
1203:             row.j   = ey;
1204:             row.k   = ez;
1205:             row.loc = DMSTAG_LEFT;
1206:             row.c   = c;
1207:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1208:             val = (15.0 + c) * x_val * x_val * x_val;
1209:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1210:           }
1211:         }
1212:         if (ex < N[0] && ey < N[1]) {
1213:           for (PetscInt c = 0; c < dof[2]; ++c) {
1214:             DMStagStencil row;
1215:             PetscScalar   x_val, val;

1217:             row.i   = ex;
1218:             row.j   = ey;
1219:             row.k   = ez;
1220:             row.loc = DMSTAG_BACK;
1221:             row.c   = c;
1222:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1223:             val = (15.0 + c) * x_val * x_val * x_val;
1224:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1225:           }
1226:         }
1227:         if (ex < N[0] && ey < N[1] && ez < N[2]) {
1228:           for (PetscInt c = 0; c < dof[3]; ++c) {
1229:             DMStagStencil row;
1230:             PetscScalar   x_val, val;

1232:             row.i   = ex;
1233:             row.j   = ey;
1234:             row.k   = ez;
1235:             row.loc = DMSTAG_ELEMENT;
1236:             row.c   = c;
1237:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1238:             val = (20.0 + c) * x_val * x_val * x_val;
1239:             DMStagVecSetValuesStencil(dm, f, 1, &row, &val, ADD_VALUES);
1240:           }
1241:         }
1242:       }
1243:     }
1244:   }

1246:   /* Add additional terms fully coupling one interior element to another */
1247:   {
1248:     PetscMPIInt rank;

1250:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1251:     if (rank == 0) {
1252:       PetscInt       epe;
1253:       DMStagStencil *row, *col;

1255:       DMStagGetEntriesPerElement(dm, &epe);
1256:       PetscMalloc1(epe, &row);
1257:       PetscMalloc1(epe, &col);
1258:       for (PetscInt i = 0; i < epe; ++i) {
1259:         row[i].i = 0;
1260:         row[i].j = 0;
1261:         row[i].k = 0;
1262:         col[i].i = 0;
1263:         col[i].j = 0;
1264:         col[i].k = 1;
1265:       }

1267:       {
1268:         PetscInt nrows = 0;

1270:         for (PetscInt c = 0; c < dof[0]; ++c) {
1271:           row[nrows].c   = c;
1272:           row[nrows].loc = DMSTAG_BACK_DOWN_LEFT;
1273:           ++nrows;
1274:         }
1275:         for (PetscInt c = 0; c < dof[1]; ++c) {
1276:           row[nrows].c   = c;
1277:           row[nrows].loc = DMSTAG_DOWN_LEFT;
1278:           ++nrows;
1279:         }
1280:         for (PetscInt c = 0; c < dof[1]; ++c) {
1281:           row[nrows].c   = c;
1282:           row[nrows].loc = DMSTAG_BACK_LEFT;
1283:           ++nrows;
1284:         }
1285:         for (PetscInt c = 0; c < dof[1]; ++c) {
1286:           row[nrows].c   = c;
1287:           row[nrows].loc = DMSTAG_BACK_DOWN;
1288:           ++nrows;
1289:         }
1290:         for (PetscInt c = 0; c < dof[2]; ++c) {
1291:           row[nrows].c   = c;
1292:           row[nrows].loc = DMSTAG_LEFT;
1293:           ++nrows;
1294:         }
1295:         for (PetscInt c = 0; c < dof[2]; ++c) {
1296:           row[nrows].c   = c;
1297:           row[nrows].loc = DMSTAG_DOWN;
1298:           ++nrows;
1299:         }
1300:         for (PetscInt c = 0; c < dof[2]; ++c) {
1301:           row[nrows].c   = c;
1302:           row[nrows].loc = DMSTAG_BACK;
1303:           ++nrows;
1304:         }
1305:         for (PetscInt c = 0; c < dof[3]; ++c) {
1306:           row[nrows].c   = c;
1307:           row[nrows].loc = DMSTAG_ELEMENT;
1308:           ++nrows;
1309:         }
1310:       }

1312:       {
1313:         PetscInt ncols = 0;

1315:         for (PetscInt c = 0; c < dof[0]; ++c) {
1316:           col[ncols].c   = c;
1317:           col[ncols].loc = DMSTAG_BACK_DOWN_LEFT;
1318:           ++ncols;
1319:         }
1320:         for (PetscInt c = 0; c < dof[1]; ++c) {
1321:           col[ncols].c   = c;
1322:           col[ncols].loc = DMSTAG_DOWN_LEFT;
1323:           ++ncols;
1324:         }
1325:         for (PetscInt c = 0; c < dof[1]; ++c) {
1326:           col[ncols].c   = c;
1327:           col[ncols].loc = DMSTAG_BACK_LEFT;
1328:           ++ncols;
1329:         }
1330:         for (PetscInt c = 0; c < dof[1]; ++c) {
1331:           col[ncols].c   = c;
1332:           col[ncols].loc = DMSTAG_BACK_DOWN;
1333:           ++ncols;
1334:         }
1335:         for (PetscInt c = 0; c < dof[2]; ++c) {
1336:           col[ncols].c   = c;
1337:           col[ncols].loc = DMSTAG_LEFT;
1338:           ++ncols;
1339:         }
1340:         for (PetscInt c = 0; c < dof[2]; ++c) {
1341:           col[ncols].c   = c;
1342:           col[ncols].loc = DMSTAG_DOWN;
1343:           ++ncols;
1344:         }
1345:         for (PetscInt c = 0; c < dof[2]; ++c) {
1346:           col[ncols].c   = c;
1347:           col[ncols].loc = DMSTAG_BACK;
1348:           ++ncols;
1349:         }
1350:         for (PetscInt c = 0; c < dof[3]; ++c) {
1351:           col[ncols].c   = c;
1352:           col[ncols].loc = DMSTAG_ELEMENT;
1353:           ++ncols;
1354:         }
1355:       }

1357:       for (PetscInt i = 0; i < epe; ++i) {
1358:         for (PetscInt j = 0; j < epe; ++j) {
1359:           PetscScalar x_val, val;

1361:           DMStagVecGetValuesStencil(dm, x_local, 1, &col[j], &x_val);
1362:           val = (10 * i + j) * x_val * x_val * x_val;
1363:           DMStagVecSetValuesStencil(dm, f, 1, &row[i], &val, ADD_VALUES);
1364:         }
1365:       }
1366:       PetscFree(row);
1367:       PetscFree(col);
1368:     }
1369:   }
1370:   DMRestoreLocalVector(dm, &x_local);
1371:   VecAssemblyBegin(f);
1372:   VecAssemblyEnd(f);
1373:   return 0;
1374: }

1376: PetscErrorCode FormJacobian3D(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx)
1377: {
1378:   PetscInt start[3], n[3], n_extra[3], N[3], dof[4];
1379:   Vec      x_local;
1380:   DM       dm;

1382:   (void)ctx;
1383:   SNESGetDM(snes, &dm);
1384:   DMGetLocalVector(dm, &x_local);
1385:   DMGlobalToLocal(dm, x, INSERT_VALUES, x_local);
1386:   DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &n_extra[0], &n_extra[1], &n_extra[2]);
1387:   DMStagGetGlobalSizes(dm, &N[0], &N[1], &N[2]);
1388:   DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], &dof[3]);
1389:   MatZeroEntries(Amat);
1390:   for (PetscInt ez = start[2]; ez < start[2] + n[2] + n_extra[2]; ++ez) {
1391:     for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
1392:       for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
1393:         for (PetscInt c = 0; c < dof[0]; ++c) {
1394:           DMStagStencil row;
1395:           PetscScalar   x_val, val;

1397:           row.i   = ex;
1398:           row.j   = ey;
1399:           row.k   = ez;
1400:           row.loc = DMSTAG_BACK_DOWN_LEFT;
1401:           row.c   = c;
1402:           DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1403:           val = 3.0 * (5.0 + c) * x_val * x_val;
1404:           DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1405:         }
1406:         if (ez < N[2]) {
1407:           for (PetscInt c = 0; c < dof[1]; ++c) {
1408:             DMStagStencil row;
1409:             PetscScalar   x_val, val;

1411:             row.i   = ex;
1412:             row.j   = ey;
1413:             row.k   = ez;
1414:             row.loc = DMSTAG_DOWN_LEFT;
1415:             row.c   = c;
1416:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1417:             val = 3.0 * (50.0 + c) * x_val * x_val;
1418:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1419:           }
1420:         }
1421:         if (ey < N[1]) {
1422:           for (PetscInt c = 0; c < dof[1]; ++c) {
1423:             DMStagStencil row;
1424:             PetscScalar   x_val, val;

1426:             row.i   = ex;
1427:             row.j   = ey;
1428:             row.k   = ez;
1429:             row.loc = DMSTAG_BACK_LEFT;
1430:             row.c   = c;
1431:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1432:             val = 3.0 * (55.0 + c) * x_val * x_val;
1433:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1434:           }
1435:         }
1436:         if (ex < N[0]) {
1437:           for (PetscInt c = 0; c < dof[1]; ++c) {
1438:             DMStagStencil row;
1439:             PetscScalar   x_val, val;

1441:             row.i   = ex;
1442:             row.j   = ey;
1443:             row.k   = ez;
1444:             row.loc = DMSTAG_BACK_DOWN;
1445:             row.c   = c;
1446:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1447:             val = 3.0 * (60.0 + c) * x_val * x_val;
1448:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1449:           }
1450:         }
1451:         if (ex < N[0] && ez < N[2]) {
1452:           for (PetscInt c = 0; c < dof[2]; ++c) {
1453:             DMStagStencil row;
1454:             PetscScalar   x_val, val;

1456:             row.i   = ex;
1457:             row.j   = ey;
1458:             row.k   = ez;
1459:             row.loc = DMSTAG_DOWN;
1460:             row.c   = c;
1461:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1462:             val = 3.0 * (10.0 + c) * x_val * x_val;
1463:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1464:           }
1465:         }
1466:         if (ey < N[1] && ez < N[2]) {
1467:           for (PetscInt c = 0; c < dof[2]; ++c) {
1468:             DMStagStencil row;
1469:             PetscScalar   x_val, val;

1471:             row.i   = ex;
1472:             row.j   = ey;
1473:             row.k   = ez;
1474:             row.loc = DMSTAG_LEFT;
1475:             row.c   = c;
1476:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1477:             val = 3.0 * (15.0 + c) * x_val * x_val;
1478:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1479:           }
1480:         }
1481:         if (ex < N[0] && ey < N[1]) {
1482:           for (PetscInt c = 0; c < dof[2]; ++c) {
1483:             DMStagStencil row;
1484:             PetscScalar   x_val, val;

1486:             row.i   = ex;
1487:             row.j   = ey;
1488:             row.k   = ez;
1489:             row.loc = DMSTAG_BACK;
1490:             row.c   = c;
1491:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1492:             val = 3.0 * (15.0 + c) * x_val * x_val;
1493:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1494:           }
1495:         }
1496:         if (ex < N[0] && ey < N[1] && ez < N[2]) {
1497:           for (PetscInt c = 0; c < dof[3]; ++c) {
1498:             DMStagStencil row;
1499:             PetscScalar   x_val, val;

1501:             row.i   = ex;
1502:             row.j   = ey;
1503:             row.k   = ez;
1504:             row.loc = DMSTAG_ELEMENT;
1505:             row.c   = c;
1506:             DMStagVecGetValuesStencil(dm, x_local, 1, &row, &x_val);
1507:             val = 3.0 * (20.0 + c) * x_val * x_val;
1508:             DMStagMatSetValuesStencil(dm, Amat, 1, &row, 1, &row, &val, ADD_VALUES);
1509:           }
1510:         }
1511:       }
1512:     }
1513:   }

1515:   /* Add additional terms fully coupling one interior element to another */
1516:   {
1517:     PetscMPIInt rank;

1519:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1520:     if (rank == 0) {
1521:       PetscInt       epe;
1522:       DMStagStencil *row, *col;

1524:       DMStagGetEntriesPerElement(dm, &epe);
1525:       PetscMalloc1(epe, &row);
1526:       PetscMalloc1(epe, &col);
1527:       for (PetscInt i = 0; i < epe; ++i) {
1528:         row[i].i = 0;
1529:         row[i].j = 0;
1530:         row[i].k = 0;
1531:         col[i].i = 0;
1532:         col[i].j = 0;
1533:         col[i].k = 1;
1534:       }

1536:       {
1537:         PetscInt nrows = 0;

1539:         for (PetscInt c = 0; c < dof[0]; ++c) {
1540:           row[nrows].c   = c;
1541:           row[nrows].loc = DMSTAG_BACK_DOWN_LEFT;
1542:           ++nrows;
1543:         }
1544:         for (PetscInt c = 0; c < dof[1]; ++c) {
1545:           row[nrows].c   = c;
1546:           row[nrows].loc = DMSTAG_DOWN_LEFT;
1547:           ++nrows;
1548:         }
1549:         for (PetscInt c = 0; c < dof[1]; ++c) {
1550:           row[nrows].c   = c;
1551:           row[nrows].loc = DMSTAG_BACK_LEFT;
1552:           ++nrows;
1553:         }
1554:         for (PetscInt c = 0; c < dof[1]; ++c) {
1555:           row[nrows].c   = c;
1556:           row[nrows].loc = DMSTAG_BACK_DOWN;
1557:           ++nrows;
1558:         }
1559:         for (PetscInt c = 0; c < dof[2]; ++c) {
1560:           row[nrows].c   = c;
1561:           row[nrows].loc = DMSTAG_LEFT;
1562:           ++nrows;
1563:         }
1564:         for (PetscInt c = 0; c < dof[2]; ++c) {
1565:           row[nrows].c   = c;
1566:           row[nrows].loc = DMSTAG_DOWN;
1567:           ++nrows;
1568:         }
1569:         for (PetscInt c = 0; c < dof[2]; ++c) {
1570:           row[nrows].c   = c;
1571:           row[nrows].loc = DMSTAG_BACK;
1572:           ++nrows;
1573:         }
1574:         for (PetscInt c = 0; c < dof[3]; ++c) {
1575:           row[nrows].c   = c;
1576:           row[nrows].loc = DMSTAG_ELEMENT;
1577:           ++nrows;
1578:         }
1579:       }

1581:       {
1582:         PetscInt ncols = 0;

1584:         for (PetscInt c = 0; c < dof[0]; ++c) {
1585:           col[ncols].c   = c;
1586:           col[ncols].loc = DMSTAG_BACK_DOWN_LEFT;
1587:           ++ncols;
1588:         }
1589:         for (PetscInt c = 0; c < dof[1]; ++c) {
1590:           col[ncols].c   = c;
1591:           col[ncols].loc = DMSTAG_DOWN_LEFT;
1592:           ++ncols;
1593:         }
1594:         for (PetscInt c = 0; c < dof[1]; ++c) {
1595:           col[ncols].c   = c;
1596:           col[ncols].loc = DMSTAG_BACK_LEFT;
1597:           ++ncols;
1598:         }
1599:         for (PetscInt c = 0; c < dof[1]; ++c) {
1600:           col[ncols].c   = c;
1601:           col[ncols].loc = DMSTAG_BACK_DOWN;
1602:           ++ncols;
1603:         }
1604:         for (PetscInt c = 0; c < dof[2]; ++c) {
1605:           col[ncols].c   = c;
1606:           col[ncols].loc = DMSTAG_LEFT;
1607:           ++ncols;
1608:         }
1609:         for (PetscInt c = 0; c < dof[2]; ++c) {
1610:           col[ncols].c   = c;
1611:           col[ncols].loc = DMSTAG_DOWN;
1612:           ++ncols;
1613:         }
1614:         for (PetscInt c = 0; c < dof[2]; ++c) {
1615:           col[ncols].c   = c;
1616:           col[ncols].loc = DMSTAG_BACK;
1617:           ++ncols;
1618:         }
1619:         for (PetscInt c = 0; c < dof[3]; ++c) {
1620:           col[ncols].c   = c;
1621:           col[ncols].loc = DMSTAG_ELEMENT;
1622:           ++ncols;
1623:         }
1624:       }

1626:       for (PetscInt i = 0; i < epe; ++i) {
1627:         for (PetscInt j = 0; j < epe; ++j) {
1628:           PetscScalar x_val, val;

1630:           DMStagVecGetValuesStencil(dm, x_local, 1, &col[j], &x_val);
1631:           val = 3.0 * (10 * i + j) * x_val * x_val;
1632:           DMStagMatSetValuesStencil(dm, Amat, 1, &row[i], 1, &col[j], &val, ADD_VALUES);
1633:         }
1634:       }
1635:       PetscFree(row);
1636:       PetscFree(col);
1637:     }
1638:   }
1639:   DMRestoreLocalVector(dm, &x_local);
1640:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
1641:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
1643:   return 0;
1644: }

1646: int main(int argc, char **argv)
1647: {
1648:   DM        dm;
1649:   PetscInt  dim;
1650:   PetscBool no_coupling;
1651:   Vec       x, b;
1652:   SNES      snes;

1655:   PetscInitialize(&argc, &argv, (char *)0, help);
1656:   dim = 3;
1657:   PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
1658:   no_coupling = PETSC_FALSE;
1659:   PetscOptionsGetBool(NULL, NULL, "-no_coupling", &no_coupling, NULL);

1661:   switch (dim) {
1662:   case 1:
1663:     DMStagCreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, &dm);
1664:     break;
1665:   case 2:
1666:     DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 3, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, &dm);
1667:     break;
1668:   case 3:
1669:     DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, NULL, &dm);
1670:     break;
1671:   default:
1672:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
1673:   }
1674:   DMSetFromOptions(dm);
1675:   DMSetUp(dm);

1677:   SNESCreate(PETSC_COMM_WORLD, &snes);
1678:   SNESSetDM(snes, dm);
1679:   if (no_coupling) {
1680:     switch (dim) {
1681:     case 1:
1682:       SNESSetFunction(snes, NULL, FormFunction1DNoCoupling, NULL);
1683:       SNESSetJacobian(snes, NULL, NULL, FormJacobian1DNoCoupling, NULL);
1684:       break;
1685:     case 2:
1686:       SNESSetFunction(snes, NULL, FormFunction2DNoCoupling, NULL);
1687:       SNESSetJacobian(snes, NULL, NULL, FormJacobian2DNoCoupling, NULL);
1688:       break;
1689:     case 3:
1690:       SNESSetFunction(snes, NULL, FormFunction3DNoCoupling, NULL);
1691:       SNESSetJacobian(snes, NULL, NULL, FormJacobian3DNoCoupling, NULL);
1692:       break;
1693:     default:
1694:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
1695:     }
1696:   } else {
1697:     switch (dim) {
1698:     case 1:
1699:       SNESSetFunction(snes, NULL, FormFunction1D, NULL);
1700:       SNESSetJacobian(snes, NULL, NULL, FormJacobian1D, NULL);
1701:       break;
1702:     case 2:
1703:       SNESSetFunction(snes, NULL, FormFunction2D, NULL);
1704:       SNESSetJacobian(snes, NULL, NULL, FormJacobian2D, NULL);
1705:       break;
1706:     case 3:
1707:       SNESSetFunction(snes, NULL, FormFunction3D, NULL);
1708:       SNESSetJacobian(snes, NULL, NULL, FormJacobian3D, NULL);
1709:       break;
1710:     default:
1711:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
1712:     }
1713:   }
1714:   SNESSetFromOptions(snes);

1716:   DMCreateGlobalVector(dm, &x);
1717:   VecDuplicate(x, &b);
1718:   VecSet(x, 2.0); // Initial guess
1719:   VecSet(b, 0.0); // RHS
1720:   SNESSolve(snes, b, x);

1722:   SNESDestroy(&snes);
1723:   VecDestroy(&x);
1724:   VecDestroy(&b);
1725:   DMDestroy(&dm);
1726:   PetscFinalize();
1727:   return 0;
1728: }

1730: /*TEST

1732:    test:
1733:       suffix: 1d_no_coupling
1734:       nsize: {{1 2}separate output}
1735:       args: -dim 1 -no_coupling -stag_stencil_type none -pc_type jacobi -snes_converged_reason -snes_test_jacobian -stag_dof_0 {{1 2}separate output} -stag_dof_1 {{1 2}separate output} -snes_max_it 2
1736:    test:
1737:       suffix: 1d_test_jac
1738:       nsize: {{1 2}separate output}
1739:       args: -dim 1 -stag_stencil_width {{0 1}separate output} -pc_type jacobi -snes_converged_reason -snes_test_jacobian -snes_max_it 2
1740:    test:
1741:       suffix: 1d_fd_coloring
1742:       nsize: {{1 2}separate output}
1743:       args: -dim 1 -stag_stencil_width {{0 1 2}separate output} -pc_type jacobi -snes_converged_reason -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type {{natural sl}} -snes_max_it 2
1744:    test:
1745:       suffix: 1d_periodic
1746:       nsize: {{1 2}separate output}
1747:       args: -dim 1 -stag_boundary_type_x periodic -stag_stencil_width {{1 2}separate output} -pc_type jacobi -snes_converged_reason -snes_test_jacobian -snes_max_it 2
1748:    test:
1749:       suffix: 1d_multidof
1750:       nsize: 2
1751:       args: -dim 1 -stag_stencil_width 2 -stag_dof_0 2 -stag_dof_1 3 -pc_type jacobi -snes_converged_reason -snes_test_jacobian -snes_max_it 2
1752:    test:
1753:       suffix: 2d_no_coupling
1754:       nsize: {{1 4}separate output}
1755:       args: -dim 2 -no_coupling -stag_stencil_type none -pc_type jacobi -snes_test_jacobian -stag_dof_0 {{1 2}separate output} -stag_dof_1 {{1 2}separate output} -stag_dof_2 {{1 2}separate output} -snes_max_it 2
1756:    test:
1757:       suffix: 3d_no_coupling
1758:       nsize: 2
1759:       args: -dim 3 -no_coupling -stag_stencil_type none -pc_type jacobi -snes_test_jacobian -stag_dof_0 2 -stag_dof_1 2 -stag_dof_2 2 -stag_dof_3 2 -snes_max_it 2
1760:    test:
1761:       suffix: 2d_fd_coloring
1762:       nsize: {{1 2}separate output}
1763:       args: -dim 2 -stag_stencil_width {{1 2}separate output} -pc_type jacobi -snes_converged_reason -snes_fd_color -snes_fd_color_use_mat -stag_stencil_type {{star box}separate output} -snes_max_it 2
1764:    test:
1765:       suffix: 3d_fd_coloring
1766:       nsize: {{1 2}separate output}
1767:       args: -dim 3 -stag_stencil_width {{1 2}separate output} -pc_type jacobi -snes_converged_reason -snes_fd_color -snes_fd_color_use_mat -stag_stencil_type {{star box}separate output} -snes_max_it 2
1768: TEST*/