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*/