Actual source code: ex15.c
1: static char help[] = "Test DMStag default MG components, on a Stokes-like system.\n\n";
3: #include <petscdm.h>
4: #include <petscdmstag.h>
5: #include <petscksp.h>
7: PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b);
9: int main(int argc, char **argv)
10: {
11: DM dm;
12: PetscInt dim;
13: PetscBool flg;
14: KSP ksp;
15: Mat A;
16: Vec x, b;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
20: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, &flg));
21: if (!flg) {
22: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Supply -dim option\n"));
23: return 1;
24: }
25: if (dim == 1) {
26: PetscCall(DMStagCreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 4, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, &dm));
27: } else if (dim == 2) {
28: PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 0, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, &dm));
29: } else if (dim == 3) {
30: PetscCall(DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 4, 4, 4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 0, 0, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, NULL, &dm));
31: } else {
32: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Supply -dim option with value 1, 2, or 3\n"));
33: return 1;
34: }
35: PetscCall(DMSetFromOptions(dm));
36: PetscCall(DMSetUp(dm));
38: /* Directly create a coarsened DM and transfer operators */
39: {
40: DM dmCoarse;
41: PetscCall(DMCoarsen(dm, MPI_COMM_NULL, &dmCoarse));
42: {
43: Mat Ai;
44: PetscCall(DMCreateInterpolation(dmCoarse, dm, &Ai, NULL));
45: PetscCall(MatDestroy(&Ai));
46: }
47: {
48: Mat Ar;
49: PetscCall(DMCreateRestriction(dmCoarse, dm, &Ar));
50: PetscCall(MatDestroy(&Ar));
51: }
52: PetscCall(DMDestroy(&dmCoarse));
53: }
55: /* Create and solve a system */
56: PetscCall(CreateSystem(dm, &A, &b));
57: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
58: PetscCall(KSPSetOperators(ksp, A, A));
59: PetscCall(KSPSetDM(ksp, dm));
60: PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
61: PetscCall(KSPSetFromOptions(ksp));
63: PetscCall(VecDuplicate(b, &x));
64: PetscCall(KSPSolve(ksp, b, x));
66: PetscCall(VecDestroy(&x));
67: PetscCall(VecDestroy(&b));
68: PetscCall(KSPDestroy(&ksp));
69: PetscCall(MatDestroy(&A));
70: PetscCall(DMDestroy(&dm));
71: PetscCall(PetscFinalize());
72: return 0;
73: }
75: /* Note: unlike in the 2D case, this does not include reasonable scaling and so will not work well */
76: PetscErrorCode CreateSystem1d(DM dm, Mat *A, Vec *b)
77: {
78: PetscInt dim;
80: PetscFunctionBeginUser;
81: PetscCall(DMGetDimension(dm, &dim));
82: PetscCall(DMCreateMatrix(dm, A));
83: PetscCall(DMCreateGlobalVector(dm, b));
84: if (dim == 1) {
85: PetscInt e, start, n, N;
86: PetscBool isFirstRank, isLastRank;
87: PetscScalar h;
88: PetscCall(DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, NULL, NULL, NULL));
89: PetscCall(DMStagGetGlobalSizes(dm, &N, NULL, NULL));
90: h = 1.0 / N;
91: PetscCall(DMStagGetIsFirstRank(dm, &isFirstRank, NULL, NULL));
92: PetscCall(DMStagGetIsLastRank(dm, &isLastRank, NULL, NULL));
93: for (e = start; e < start + n; ++e) {
94: DMStagStencil pos[3];
95: PetscScalar val[3];
96: PetscInt idxLoc;
98: if (isFirstRank && e == start) {
99: /* Fix first pressure node to eliminate nullspace */
100: idxLoc = 0;
101: pos[idxLoc].i = e;
102: pos[idxLoc].loc = DMSTAG_ELEMENT;
103: pos[idxLoc].c = 0;
104: val[idxLoc] = 1.0; /* 0 pressure forcing term (physical) */
105: ++idxLoc;
106: } else {
107: idxLoc = 0;
108: pos[idxLoc].i = e;
109: pos[idxLoc].loc = DMSTAG_ELEMENT;
110: pos[idxLoc].c = 0;
111: val[idxLoc] = 0.0; /* 0 pressure forcing term (physical) */
112: ++idxLoc;
113: }
115: if (isFirstRank && e == start) {
116: pos[idxLoc].i = e;
117: pos[idxLoc].loc = DMSTAG_LEFT;
118: pos[idxLoc].c = 0;
119: val[idxLoc] = 3.0; /* fixed left BC */
120: ++idxLoc;
121: } else {
122: pos[idxLoc].i = e;
123: pos[idxLoc].loc = DMSTAG_LEFT;
124: pos[idxLoc].c = 0;
125: val[idxLoc] = 1.0; /* constant forcing */
126: ++idxLoc;
127: }
128: if (isLastRank && e == start + n - 1) {
129: /* Special case on right boundary (in addition to usual case) */
130: pos[idxLoc].i = e; /* This element in the 1d ordering */
131: pos[idxLoc].loc = DMSTAG_RIGHT;
132: pos[idxLoc].c = 0;
133: val[idxLoc] = 3.0; /* fixed right BC */
134: ++idxLoc;
135: }
136: PetscCall(DMStagVecSetValuesStencil(dm, *b, idxLoc, pos, val, INSERT_VALUES));
137: }
139: for (e = start; e < start + n; ++e) {
140: if (isFirstRank && e == start) {
141: DMStagStencil row;
142: PetscScalar val;
144: row.i = e;
145: row.loc = DMSTAG_LEFT;
146: row.c = 0;
147: val = 1.0;
148: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
149: } else {
150: DMStagStencil row, col[5];
151: PetscScalar val[5];
153: row.i = e;
154: row.loc = DMSTAG_LEFT;
155: row.c = 0;
157: col[0].i = e;
158: col[0].loc = DMSTAG_ELEMENT;
159: col[0].c = 0;
161: col[1].i = e - 1;
162: col[1].loc = DMSTAG_ELEMENT;
163: col[1].c = 0;
165: val[0] = -1.0 / h;
166: val[1] = 1.0 / h;
168: col[2].i = e;
169: col[2].loc = DMSTAG_LEFT;
170: col[2].c = 0;
171: val[2] = -2.0 / (h * h);
173: col[3].i = e - 1;
174: col[3].loc = DMSTAG_LEFT;
175: col[3].c = 0;
176: val[3] = 1.0 / (h * h);
178: col[4].i = e + 1;
179: col[4].loc = DMSTAG_LEFT;
180: col[4].c = 0;
181: val[4] = 1.0 / (h * h);
183: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 5, col, val, INSERT_VALUES));
184: }
186: /* Additional velocity point (BC) on the right */
187: if (isLastRank && e == start + n - 1) {
188: DMStagStencil row;
189: PetscScalar val;
191: row.i = e;
192: row.loc = DMSTAG_RIGHT;
193: row.c = 0;
194: val = 1.0;
195: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
196: }
198: /* Equation on pressure (element) variables */
199: if (isFirstRank && e == 0) {
200: /* Fix first pressure node to eliminate nullspace */
201: DMStagStencil row;
202: PetscScalar val;
204: row.i = e;
205: row.loc = DMSTAG_ELEMENT;
206: row.c = 0;
207: val = 1.0;
209: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
210: } else {
211: DMStagStencil row, col[2];
212: PetscScalar val[2];
214: row.i = e;
215: row.loc = DMSTAG_ELEMENT;
216: row.c = 0;
218: col[0].i = e;
219: col[0].loc = DMSTAG_LEFT;
220: col[0].c = 0;
222: col[1].i = e;
223: col[1].loc = DMSTAG_RIGHT;
224: col[1].c = 0;
226: val[0] = -1.0 / h;
227: val[1] = 1.0 / h;
229: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 2, col, val, INSERT_VALUES));
230: }
231: }
232: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
233: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
234: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
235: PetscCall(VecAssemblyBegin(*b));
236: PetscCall(VecAssemblyEnd(*b));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PetscErrorCode CreateSystem2d(DM dm, Mat *A, Vec *b)
241: {
242: PetscInt N[2];
243: PetscBool isLastRankx, isLastRanky, isFirstRankx, isFirstRanky;
244: PetscInt ex, ey, startx, starty, nx, ny;
245: PetscReal hx, hy, dv;
247: PetscFunctionBeginUser;
248: PetscCall(DMCreateMatrix(dm, A));
249: PetscCall(DMCreateGlobalVector(dm, b));
250: PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
251: PetscCall(DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL));
252: PetscCall(DMStagGetIsLastRank(dm, &isLastRankx, &isLastRanky, NULL));
253: PetscCall(DMStagGetIsFirstRank(dm, &isFirstRankx, &isFirstRanky, NULL));
254: hx = 1.0 / N[0];
255: hy = 1.0 / N[1];
256: dv = hx * hy;
258: for (ey = starty; ey < starty + ny; ++ey) {
259: for (ex = startx; ex < startx + nx; ++ex) {
260: if (ex == N[0] - 1) {
261: /* Right Boundary velocity Dirichlet */
262: DMStagStencil row;
263: PetscScalar valRhs;
264: const PetscScalar valA = 1.0;
265: row.i = ex;
266: row.j = ey;
267: row.loc = DMSTAG_RIGHT;
268: row.c = 0;
269: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
270: valRhs = 0.0; /* zero Dirichlet */
271: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
272: }
273: if (ey == N[1] - 1) {
274: /* Top boundary velocity Dirichlet */
275: DMStagStencil row;
276: PetscScalar valRhs;
277: const PetscScalar valA = 1.0;
278: row.i = ex;
279: row.j = ey;
280: row.loc = DMSTAG_UP;
281: row.c = 0;
282: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
283: valRhs = 0.0; /* zero Diri */
284: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
285: }
287: if (ey == 0) {
288: /* Bottom boundary velocity Dirichlet */
289: DMStagStencil row;
290: PetscScalar valRhs;
291: const PetscScalar valA = 1.0;
292: row.i = ex;
293: row.j = ey;
294: row.loc = DMSTAG_DOWN;
295: row.c = 0;
296: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
297: valRhs = 0.0; /* zero Diri */
298: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
299: } else {
300: /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
301: DMStagStencil row, col[7];
302: PetscScalar valA[7], valRhs;
303: PetscInt nEntries;
305: row.i = ex;
306: row.j = ey;
307: row.loc = DMSTAG_DOWN;
308: row.c = 0;
309: if (ex == 0) {
310: nEntries = 6;
311: col[0].i = ex;
312: col[0].j = ey;
313: col[0].loc = DMSTAG_DOWN;
314: col[0].c = 0;
315: valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
316: col[1].i = ex;
317: col[1].j = ey - 1;
318: col[1].loc = DMSTAG_DOWN;
319: col[1].c = 0;
320: valA[1] = dv * 1.0 / (hy * hy);
321: col[2].i = ex;
322: col[2].j = ey + 1;
323: col[2].loc = DMSTAG_DOWN;
324: col[2].c = 0;
325: valA[2] = dv * 1.0 / (hy * hy);
326: /* Missing left element */
327: col[3].i = ex + 1;
328: col[3].j = ey;
329: col[3].loc = DMSTAG_DOWN;
330: col[3].c = 0;
331: valA[3] = dv * 1.0 / (hx * hx);
332: col[4].i = ex;
333: col[4].j = ey - 1;
334: col[4].loc = DMSTAG_ELEMENT;
335: col[4].c = 0;
336: valA[4] = dv * 1.0 / hy;
337: col[5].i = ex;
338: col[5].j = ey;
339: col[5].loc = DMSTAG_ELEMENT;
340: col[5].c = 0;
341: valA[5] = -dv * 1.0 / hy;
342: } else if (ex == N[0] - 1) {
343: /* Right boundary y velocity stencil */
344: nEntries = 6;
345: col[0].i = ex;
346: col[0].j = ey;
347: col[0].loc = DMSTAG_DOWN;
348: col[0].c = 0;
349: valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
350: col[1].i = ex;
351: col[1].j = ey - 1;
352: col[1].loc = DMSTAG_DOWN;
353: col[1].c = 0;
354: valA[1] = dv * 1.0 / (hy * hy);
355: col[2].i = ex;
356: col[2].j = ey + 1;
357: col[2].loc = DMSTAG_DOWN;
358: col[2].c = 0;
359: valA[2] = dv * 1.0 / (hy * hy);
360: col[3].i = ex - 1;
361: col[3].j = ey;
362: col[3].loc = DMSTAG_DOWN;
363: col[3].c = 0;
364: valA[3] = dv * 1.0 / (hx * hx);
365: /* Missing right element */
366: col[4].i = ex;
367: col[4].j = ey - 1;
368: col[4].loc = DMSTAG_ELEMENT;
369: col[4].c = 0;
370: valA[4] = dv * 1.0 / hy;
371: col[5].i = ex;
372: col[5].j = ey;
373: col[5].loc = DMSTAG_ELEMENT;
374: col[5].c = 0;
375: valA[5] = -dv * 1.0 / hy;
376: } else {
377: nEntries = 7;
378: col[0].i = ex;
379: col[0].j = ey;
380: col[0].loc = DMSTAG_DOWN;
381: col[0].c = 0;
382: valA[0] = -dv * 2.0 / (hx * hx) - dv * 2.0 / (hy * hy);
383: col[1].i = ex;
384: col[1].j = ey - 1;
385: col[1].loc = DMSTAG_DOWN;
386: col[1].c = 0;
387: valA[1] = dv * 1.0 / (hy * hy);
388: col[2].i = ex;
389: col[2].j = ey + 1;
390: col[2].loc = DMSTAG_DOWN;
391: col[2].c = 0;
392: valA[2] = dv * 1.0 / (hy * hy);
393: col[3].i = ex - 1;
394: col[3].j = ey;
395: col[3].loc = DMSTAG_DOWN;
396: col[3].c = 0;
397: valA[3] = dv * 1.0 / (hx * hx);
398: col[4].i = ex + 1;
399: col[4].j = ey;
400: col[4].loc = DMSTAG_DOWN;
401: col[4].c = 0;
402: valA[4] = dv * 1.0 / (hx * hx);
403: col[5].i = ex;
404: col[5].j = ey - 1;
405: col[5].loc = DMSTAG_ELEMENT;
406: col[5].c = 0;
407: valA[5] = dv * 1.0 / hy;
408: col[6].i = ex;
409: col[6].j = ey;
410: col[6].loc = DMSTAG_ELEMENT;
411: col[6].c = 0;
412: valA[6] = -dv * 1.0 / hy;
413: }
414: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, nEntries, col, valA, INSERT_VALUES));
415: valRhs = dv * 1.0; /* non-zero */
416: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
417: }
419: if (ex == 0) {
420: /* Left velocity Dirichlet */
421: DMStagStencil row;
422: PetscScalar valRhs;
423: const PetscScalar valA = 1.0;
424: row.i = ex;
425: row.j = ey;
426: row.loc = DMSTAG_LEFT;
427: row.c = 0;
428: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
429: valRhs = 0.0; /* zero Diri */
430: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
431: } else {
432: /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
433: DMStagStencil row, col[7];
434: PetscScalar valA[7], valRhs;
435: PetscInt nEntries;
436: row.i = ex;
437: row.j = ey;
438: row.loc = DMSTAG_LEFT;
439: row.c = 0;
441: if (ey == 0) {
442: nEntries = 6;
443: col[0].i = ex;
444: col[0].j = ey;
445: col[0].loc = DMSTAG_LEFT;
446: col[0].c = 0;
447: valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
448: /* missing term from element below */
449: col[1].i = ex;
450: col[1].j = ey + 1;
451: col[1].loc = DMSTAG_LEFT;
452: col[1].c = 0;
453: valA[1] = dv * 1.0 / (hy * hy);
454: col[2].i = ex - 1;
455: col[2].j = ey;
456: col[2].loc = DMSTAG_LEFT;
457: col[2].c = 0;
458: valA[2] = dv * 1.0 / (hx * hx);
459: col[3].i = ex + 1;
460: col[3].j = ey;
461: col[3].loc = DMSTAG_LEFT;
462: col[3].c = 0;
463: valA[3] = dv * 1.0 / (hx * hx);
464: col[4].i = ex - 1;
465: col[4].j = ey;
466: col[4].loc = DMSTAG_ELEMENT;
467: col[4].c = 0;
468: valA[4] = dv * 1.0 / hx;
469: col[5].i = ex;
470: col[5].j = ey;
471: col[5].loc = DMSTAG_ELEMENT;
472: col[5].c = 0;
473: valA[5] = -dv * 1.0 / hx;
474: } else if (ey == N[1] - 1) {
475: /* Top boundary x velocity stencil */
476: nEntries = 6;
477: row.i = ex;
478: row.j = ey;
479: row.loc = DMSTAG_LEFT;
480: row.c = 0;
481: col[0].i = ex;
482: col[0].j = ey;
483: col[0].loc = DMSTAG_LEFT;
484: col[0].c = 0;
485: valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
486: col[1].i = ex;
487: col[1].j = ey - 1;
488: col[1].loc = DMSTAG_LEFT;
489: col[1].c = 0;
490: valA[1] = dv * 1.0 / (hy * hy);
491: /* Missing element above term */
492: col[2].i = ex - 1;
493: col[2].j = ey;
494: col[2].loc = DMSTAG_LEFT;
495: col[2].c = 0;
496: valA[2] = dv * 1.0 / (hx * hx);
497: col[3].i = ex + 1;
498: col[3].j = ey;
499: col[3].loc = DMSTAG_LEFT;
500: col[3].c = 0;
501: valA[3] = dv * 1.0 / (hx * hx);
502: col[4].i = ex - 1;
503: col[4].j = ey;
504: col[4].loc = DMSTAG_ELEMENT;
505: col[4].c = 0;
506: valA[4] = dv * 1.0 / hx;
507: col[5].i = ex;
508: col[5].j = ey;
509: col[5].loc = DMSTAG_ELEMENT;
510: col[5].c = 0;
511: valA[5] = -dv * 1.0 / hx;
512: } else {
513: nEntries = 7;
514: col[0].i = ex;
515: col[0].j = ey;
516: col[0].loc = DMSTAG_LEFT;
517: col[0].c = 0;
518: valA[0] = -dv * 2.0 / (hx * hx) + -dv * 2.0 / (hy * hy);
519: col[1].i = ex;
520: col[1].j = ey - 1;
521: col[1].loc = DMSTAG_LEFT;
522: col[1].c = 0;
523: valA[1] = dv * 1.0 / (hy * hy);
524: col[2].i = ex;
525: col[2].j = ey + 1;
526: col[2].loc = DMSTAG_LEFT;
527: col[2].c = 0;
528: valA[2] = dv * 1.0 / (hy * hy);
529: col[3].i = ex - 1;
530: col[3].j = ey;
531: col[3].loc = DMSTAG_LEFT;
532: col[3].c = 0;
533: valA[3] = dv * 1.0 / (hx * hx);
534: col[4].i = ex + 1;
535: col[4].j = ey;
536: col[4].loc = DMSTAG_LEFT;
537: col[4].c = 0;
538: valA[4] = dv * 1.0 / (hx * hx);
539: col[5].i = ex - 1;
540: col[5].j = ey;
541: col[5].loc = DMSTAG_ELEMENT;
542: col[5].c = 0;
543: valA[5] = dv * 1.0 / hx;
544: col[6].i = ex;
545: col[6].j = ey;
546: col[6].loc = DMSTAG_ELEMENT;
547: col[6].c = 0;
548: valA[6] = -dv * 1.0 / hx;
549: }
550: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, nEntries, col, valA, INSERT_VALUES));
551: valRhs = dv * 0.0; /* zero */
552: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
553: }
555: /* P equation : u_x + v_y = g
556: Note that this includes an explicit zero on the diagonal. This is only needed for
557: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
558: if (ex == 0 && ey == 0) { /* Pin the first pressure node */
559: DMStagStencil row;
560: PetscScalar valA, valRhs;
561: row.i = ex;
562: row.j = ey;
563: row.loc = DMSTAG_ELEMENT;
564: row.c = 0;
565: valA = 1.0;
566: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &valA, INSERT_VALUES));
567: valRhs = 0.0; /* zero pinned pressure */
568: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
569: } else {
570: DMStagStencil row, col[5];
571: PetscScalar valA[5], valRhs;
573: /* Note: the scaling by dv here is not principled and likely suboptimal */
574: row.i = ex;
575: row.j = ey;
576: row.loc = DMSTAG_ELEMENT;
577: row.c = 0;
578: col[0].i = ex;
579: col[0].j = ey;
580: col[0].loc = DMSTAG_LEFT;
581: col[0].c = 0;
582: valA[0] = -dv * 1.0 / hx;
583: col[1].i = ex;
584: col[1].j = ey;
585: col[1].loc = DMSTAG_RIGHT;
586: col[1].c = 0;
587: valA[1] = dv * 1.0 / hx;
588: col[2].i = ex;
589: col[2].j = ey;
590: col[2].loc = DMSTAG_DOWN;
591: col[2].c = 0;
592: valA[2] = -dv * 1.0 / hy;
593: col[3].i = ex;
594: col[3].j = ey;
595: col[3].loc = DMSTAG_UP;
596: col[3].c = 0;
597: valA[3] = dv * 1.0 / hy;
598: col[4] = row;
599: valA[4] = 0.0;
600: PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 5, col, valA, INSERT_VALUES));
601: valRhs = dv * 0.0; /* zero pressure forcing */
602: PetscCall(DMStagVecSetValuesStencil(dm, *b, 1, &row, &valRhs, INSERT_VALUES));
603: }
604: }
605: }
606: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
607: PetscCall(VecAssemblyBegin(*b));
608: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
609: PetscCall(VecAssemblyEnd(*b));
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b)
614: {
615: PetscInt dim;
617: PetscFunctionBeginUser;
618: PetscCall(DMGetDimension(dm, &dim));
619: if (dim == 1) {
620: PetscCall(CreateSystem1d(dm, A, b));
621: } else if (dim == 2) {
622: PetscCall(CreateSystem2d(dm, A, b));
623: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: /*TEST
629: test:
630: suffix: 1d_directsmooth_seq
631: nsize: 1
632: requires: suitesparse
633: args: -dim 1 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type umfpack -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason
635: test:
636: suffix: 1d_directsmooth_par
637: nsize: 4
638: requires: mumps
639: args: -dim 1 /ex15 -dim 1 -stag_grid_x 16 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type mumps -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -ksp_converged_reason
641: test:
642: suffix: 1d_fssmooth_seq
643: nsize: 1
644: requires: suitesparse
645: args: -dim 1 -stag_grid_x 256 -ksp_converged_reason -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point
647: test:
648: suffix: 1d_fssmooth_par
649: nsize: 1
650: requires: mumps !single
651: args: -dim 1 -stag_grid_x 256 -ksp_converged_reason -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point
653: test:
654: suffix: 2d_directsmooth_seq
655: nsize: 1
656: requires: suitesparse
657: args: -dim 2 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type umfpack -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason
659: test:
660: suffix: 2d_directsmooth_par
661: nsize: 4
662: requires: mumps
663: args: -dim 1 /ex15 -dim 1 -stag_grid_x 16 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type lu -mg_levels_pc_factor_mat_solver_type mumps -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -ksp_converged_reason
665: test:
666: suffix: 2d_fssmooth_seq
667: nsize: 1
668: requires: suitesparse
669: args: -dim 2 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -ksp_converged_reason
671: test:
672: suffix: 2d_fssmooth_par
673: nsize: 1
674: requires: mumps
675: args: -dim 2 -ksp_type fgmres -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type mumps -ksp_converged_reason
677: test:
678: suffix: 2d_mono_mg
679: nsize: 1
680: requires: suitesparse
681: args: -dim 2 -pc_type mg -pc_mg_galerkin -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_type SCHUR -ksp_monitor_true_residual -ksp_converged_reason -mg_levels_fieldsplit_element_pc_type jacobi -mg_levels_fieldsplit_face_pc_type jacobi -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_type umfpack -pc_mg_levels 3 -mg_levels_pc_fieldsplit_schur_fact_type UPPER -mg_levels_fieldsplit_face_pc_type sor -mg_levels_fieldsplit_face_ksp_type richardson -mg_levels_fieldsplit_face_pc_sor_symmetric -mg_levels_fieldsplit_element_ksp_type richardson -mg_levels_fieldsplit_element_pc_type none -ksp_type fgmres -stag_grid_x 48 -stag_grid_y 48 -mg_levels_fieldsplit_face_ksp_max_it 1 -mg_levels_ksp_max_it 4 -mg_levels_fieldsplit_element_ksp_type preonly -mg_levels_fieldsplit_element_pc_type none -pc_mg_cycle_type w -mg_levels_ksp_max_it 4 -mg_levels_2_ksp_max_it 8
683: TEST*/