1 static char help[] = "Test DMStag default MG components, on a Stokes-like system.\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmstag.h>
5 #include <petscksp.h>
6
7 PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b);
8
main(int argc,char ** argv)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;
17
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));
37
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 }
54
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, KSP_DMACTIVE_ALL, PETSC_FALSE));
61 PetscCall(KSPSetFromOptions(ksp));
62
63 PetscCall(VecDuplicate(b, &x));
64 PetscCall(KSPSolve(ksp, b, x));
65
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 }
74
75 /* Note: unlike in the 2D case, this does not include reasonable scaling and so will not work well */
CreateSystem1d(DM dm,Mat * A,Vec * b)76 PetscErrorCode CreateSystem1d(DM dm, Mat *A, Vec *b)
77 {
78 PetscInt dim;
79
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;
97
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 }
114
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 }
138
139 for (e = start; e < start + n; ++e) {
140 if (isFirstRank && e == start) {
141 DMStagStencil row;
142 PetscScalar val;
143
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];
152
153 row.i = e;
154 row.loc = DMSTAG_LEFT;
155 row.c = 0;
156
157 col[0].i = e;
158 col[0].loc = DMSTAG_ELEMENT;
159 col[0].c = 0;
160
161 col[1].i = e - 1;
162 col[1].loc = DMSTAG_ELEMENT;
163 col[1].c = 0;
164
165 val[0] = -1.0 / h;
166 val[1] = 1.0 / h;
167
168 col[2].i = e;
169 col[2].loc = DMSTAG_LEFT;
170 col[2].c = 0;
171 val[2] = -2.0 / (h * h);
172
173 col[3].i = e - 1;
174 col[3].loc = DMSTAG_LEFT;
175 col[3].c = 0;
176 val[3] = 1.0 / (h * h);
177
178 col[4].i = e + 1;
179 col[4].loc = DMSTAG_LEFT;
180 col[4].c = 0;
181 val[4] = 1.0 / (h * h);
182
183 PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 5, col, val, INSERT_VALUES));
184 }
185
186 /* Additional velocity point (BC) on the right */
187 if (isLastRank && e == start + n - 1) {
188 DMStagStencil row;
189 PetscScalar val;
190
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 }
197
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;
203
204 row.i = e;
205 row.loc = DMSTAG_ELEMENT;
206 row.c = 0;
207 val = 1.0;
208
209 PetscCall(DMStagMatSetValuesStencil(dm, *A, 1, &row, 1, &row, &val, INSERT_VALUES));
210 } else {
211 DMStagStencil row, col[2];
212 PetscScalar val[2];
213
214 row.i = e;
215 row.loc = DMSTAG_ELEMENT;
216 row.c = 0;
217
218 col[0].i = e;
219 col[0].loc = DMSTAG_LEFT;
220 col[0].c = 0;
221
222 col[1].i = e;
223 col[1].loc = DMSTAG_RIGHT;
224 col[1].c = 0;
225
226 val[0] = -1.0 / h;
227 val[1] = 1.0 / h;
228
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 }
239
CreateSystem2d(DM dm,Mat * A,Vec * b)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;
246
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;
257
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 }
286
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;
304
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 }
418
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;
440
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 }
554
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;
572
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 }
612
CreateSystem(DM dm,Mat * A,Vec * b)613 PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b)
614 {
615 PetscInt dim;
616
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 }
626
627 /*TEST
628
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
634
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
640
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
646
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
652
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
658
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
664
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
670
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
676
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
682
683 TEST*/
684