Actual source code: ex62.c
1: static char help[] = "Stokes Problem discretized with finite elements,\n\
2: using a parallel unstructured mesh (DMPLEX) to represent the domain.\n\n\n";
4: /*
5: For the isoviscous Stokes problem, which we discretize using the finite
6: element method on an unstructured mesh, the weak form equations are
8: < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > - < v, f > = 0
9: < q, -\nabla\cdot u > = 0
11: Viewing:
13: To produce nice output, use
15: -dm_refine 3 -dm_view hdf5:sol1.h5 -error_vec_view hdf5:sol1.h5::append -snes_view_solution hdf5:sol1.h5::append -exact_vec_view hdf5:sol1.h5::append
17: You can get a LaTeX view of the mesh, with point numbering using
19: -dm_view :mesh.tex:ascii_latex -dm_plex_view_scale 8.0
21: The data layout can be viewed using
23: -dm_petscsection_view
25: Lots of information about the FEM assembly can be printed using
27: -dm_plex_print_fem 3
28: */
30: #include <petscdmplex.h>
31: #include <petscpc.h>
32: #include <petscsnes.h>
33: #include <petscds.h>
34: #include <petscbag.h>
36: // TODO: Plot residual by fields after each smoother iterate
38: typedef enum {
39: SOL_QUADRATIC,
40: SOL_TRIG,
41: SOL_UNKNOWN
42: } SolType;
43: const char *SolTypes[] = {"quadratic", "trig", "unknown", "SolType", "SOL_", NULL};
45: typedef enum {
46: BC_ESSENTIAL,
47: BC_NITSCHE,
48: BC_UNKNOWN
49: } BCType;
50: const char *BCTypes[] = {"essential", "nitsche", "unknown", "BCType", "BC_", NULL};
52: typedef struct {
53: PetscScalar mu; /* dynamic shear viscosity */
54: PetscScalar eta; /* Nitsche penalty parameter (dimensionless) */
55: } Parameter;
57: typedef struct {
58: PetscBag bag; /* Problem parameters */
59: SolType sol; /* MMS solution */
60: BCType bc; /* Boundary condition type */
61: } AppCtx;
63: typedef struct {
64: PetscInt numPatchPoints;
65: PetscInt *patchPoints;
66: PetscInt dofsPerCell;
67: PetscInt numInteriorFacetCalls;
68: PetscInt numExteriorFacetCalls;
69: } PatchFacetTestCtx;
71: static PetscErrorCode TestPatchConstruct(PC pc, PetscInt *npatch, IS *patches[], IS *patchIterationSet, PetscCtx ctx)
72: {
73: PatchFacetTestCtx *test = (PatchFacetTestCtx *)ctx;
75: PetscFunctionBeginUser;
76: *npatch = 1;
77: PetscCall(PetscMalloc1(*npatch, patches));
78: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, test->numPatchPoints, test->patchPoints, PETSC_COPY_VALUES, *patches));
79: PetscCall(ISCreateStride(PETSC_COMM_SELF, *npatch, 0, 1, patchIterationSet));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode TestPatchFacetCallback(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt dofsArray[], const PetscInt dofsArrayWithAll[], PetscCtx ctx, PetscInt cellsPerFacet, PetscInt *numCalls)
84: {
85: PatchFacetTestCtx *test = (PatchFacetTestCtx *)ctx;
86: PetscInt numFacets;
88: PetscFunctionBeginUser;
89: PetscCall(ISGetLocalSize(facetIS, &numFacets));
90: PetscCheck(n == cellsPerFacet * numFacets * test->dofsPerCell, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Expected %" PetscInt_FMT " facet dofs, got %" PetscInt_FMT, cellsPerFacet * numFacets * test->dofsPerCell, n);
91: for (PetscInt i = 0; i < n; ++i) PetscCheck(dofsArray[i] == dofsArrayWithAll[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Facet dof maps differ at entry %" PetscInt_FMT, i);
92: ++(*numCalls);
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode TestPatchComputeFunctionInteriorFacets(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt dofsArray[], const PetscInt dofsArrayWithAll[], PetscCtx ctx)
97: {
98: PatchFacetTestCtx *test = (PatchFacetTestCtx *)ctx;
100: PetscFunctionBeginUser;
101: PetscCall(TestPatchFacetCallback(pc, point, x, f, facetIS, n, dofsArray, dofsArrayWithAll, ctx, 2, &test->numInteriorFacetCalls));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode TestPatchComputeFunctionExteriorFacets(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt dofsArray[], const PetscInt dofsArrayWithAll[], PetscCtx ctx)
106: {
107: PatchFacetTestCtx *test = (PatchFacetTestCtx *)ctx;
109: PetscFunctionBeginUser;
110: PetscCall(TestPatchFacetCallback(pc, point, x, f, facetIS, n, dofsArray, dofsArrayWithAll, ctx, 1, &test->numExteriorFacetCalls));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode TestPatchComputeFunction(PC pc, PetscInt point, Vec x, Vec f, IS cellIS, PetscInt n, const PetscInt dofsArray[], const PetscInt dofsArrayWithAll[], PetscCtx ctx)
115: {
116: PetscFunctionBeginUser;
117: PetscCall(PCPatchSetComputeFunctionInteriorFacets(pc, TestPatchComputeFunctionInteriorFacets, ctx));
118: PetscCall(PCPatchSetComputeFunctionExteriorFacets(pc, TestPatchComputeFunctionExteriorFacets, ctx));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode TestPatchComputeOperator(PC pc, PetscInt point, Vec x, Mat mat, IS cellIS, PetscInt n, const PetscInt dofsArray[], const PetscInt dofsArrayWithAll[], PetscCtx ctx)
123: {
124: PetscFunctionBeginUser;
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode TestPatchOuterFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
129: {
130: PetscFunctionBeginUser;
131: PetscCall(VecSet(f, 1.0));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: static PetscErrorCode TestPatchFacetResidual(void)
136: {
137: const PetscInt faces[2] = {2, 1};
138: const PetscInt nodesPerCellValue = 4;
139: DM dm;
140: PetscSection section;
141: PetscInt cStart, cEnd, pStart, pEnd, vStart, vEnd, cell, closureSize, numDofs, numCells;
142: PetscInt *cellNodeMap = NULL, *closure = NULL;
143: const PetscInt *cellNodeMaps[1];
144: PetscInt nodesPerCell[1] = {nodesPerCellValue};
145: DM dms[1];
146: IS isets[1];
147: SNES snes;
148: Vec x, f, rhs;
149: PatchFacetTestCtx test = {0, NULL, 0, 0, 0};
151: PetscFunctionBeginUser;
152: PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm));
153: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
154: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
155: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
156: PetscCheck(cStart == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected cell points to start at zero, got %" PetscInt_FMT, cStart);
157: numCells = cEnd - cStart;
158: test.dofsPerCell = nodesPerCellValue;
159: test.numPatchPoints = vEnd - vStart;
160: PetscCall(PetscMalloc1(test.numPatchPoints, &test.patchPoints));
161: for (PetscInt vertex = vStart; vertex < vEnd; ++vertex) test.patchPoints[vertex - vStart] = vertex;
163: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
164: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
165: for (PetscInt vertex = vStart; vertex < vEnd; ++vertex) PetscCall(PetscSectionSetDof(section, vertex, 1));
166: PetscCall(PetscSectionSetUp(section));
167: PetscCall(DMSetLocalSection(dm, section));
168: PetscCall(PetscSectionGetStorageSize(section, &numDofs));
169: PetscCall(ISCreateStride(PETSC_COMM_SELF, numDofs, 0, 1, &isets[0]));
171: PetscCall(PetscMalloc1(numCells * nodesPerCellValue, &cellNodeMap));
172: for (cell = cStart; cell < cEnd; ++cell) {
173: PetscInt numVertices = 0;
175: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
176: for (PetscInt i = 0; i < closureSize * 2; i += 2) {
177: const PetscInt point = closure[i];
179: if (point >= vStart && point < vEnd) {
180: PetscInt offset;
182: PetscCall(PetscSectionGetOffset(section, point, &offset));
183: PetscCheck(numVertices < nodesPerCellValue, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Found too many vertices in cell");
184: cellNodeMap[cell * nodesPerCellValue + numVertices++] = offset;
185: }
186: }
187: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
188: PetscCheck(numVertices == nodesPerCellValue, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected %" PetscInt_FMT " vertices in cell, got %" PetscInt_FMT, nodesPerCellValue, numVertices);
189: }
191: dms[0] = dm;
192: cellNodeMaps[0] = cellNodeMap;
193: PetscCall(VecCreateSeq(PETSC_COMM_SELF, numDofs, &x));
194: PetscCall(VecDuplicate(x, &f));
195: PetscCall(VecDuplicate(x, &rhs));
196: PetscCall(VecSet(rhs, 0.0));
197: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
198: PetscCall(SNESSetType(snes, SNESPATCH));
199: PetscCall(SNESSetDM(snes, dm));
200: PetscCall(SNESSetFunction(snes, f, TestPatchOuterFunction, NULL));
201: PetscCall(SNESPatchSetConstructType(snes, PC_PATCH_USER, TestPatchConstruct, &test));
202: PetscCall(SNESPatchSetDiscretisationInfo(snes, 1, dms, isets, nodesPerCell, cellNodeMaps, 0, NULL, 0, NULL));
203: PetscCall(SNESPatchSetComputeFunction(snes, TestPatchComputeFunction, &test));
204: PetscCall(SNESPatchSetComputeOperator(snes, TestPatchComputeOperator, NULL));
205: PetscCall(SNESSetTolerances(snes, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1, PETSC_CURRENT));
206: PetscCall(PetscOptionsSetValue(NULL, "-sub_snes_max_it", "1"));
207: PetscCall(SNESSetFromOptions(snes));
208: PetscCall(SNESSolve(snes, rhs, x));
209: PetscCheck(test.numInteriorFacetCalls > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior facet callback was not called");
210: PetscCheck(test.numExteriorFacetCalls > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Exterior facet callback was not called");
211: PetscCall(SNESDestroy(&snes));
212: PetscCall(VecDestroy(&rhs));
213: PetscCall(VecDestroy(&f));
214: PetscCall(VecDestroy(&x));
215: PetscCall(PetscFree(cellNodeMap));
216: PetscCall(ISDestroy(&isets[0]));
217: PetscCall(PetscSectionDestroy(§ion));
218: PetscCall(DMDestroy(&dm));
219: PetscCall(PetscFree(test.patchPoints));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
224: {
225: const PetscReal mu = PetscRealPart(constants[0]);
226: const PetscInt Nc = uOff[1] - uOff[0];
227: PetscInt c, d;
229: for (c = 0; c < Nc; ++c) {
230: for (d = 0; d < dim; ++d) f1[c * dim + d] = mu * (u_x[c * dim + d] + u_x[d * dim + c]);
231: f1[c * dim + c] -= u[uOff[1]];
232: }
233: }
235: static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
236: {
237: PetscInt d;
238: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] -= u_x[d * dim + d];
239: }
241: static void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
242: {
243: PetscInt d;
244: for (d = 0; d < dim; ++d) g1[d * dim + d] = -1.0; /* < q, -\nabla\cdot u > */
245: }
247: static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
248: {
249: PetscInt d;
250: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* -< \nabla\cdot v, p > */
251: }
253: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
254: {
255: const PetscReal mu = PetscRealPart(constants[0]);
256: const PetscInt Nc = uOff[1] - uOff[0];
257: PetscInt c, d;
259: for (c = 0; c < Nc; ++c) {
260: for (d = 0; d < dim; ++d) {
261: g3[((c * Nc + c) * dim + d) * dim + d] += mu; /* < \nabla v, \nabla u > */
262: g3[((c * Nc + d) * dim + d) * dim + c] += mu; /* < \nabla v, {\nabla u}^T > */
263: }
264: }
265: }
267: static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
268: {
269: const PetscReal mu = PetscRealPart(constants[0]);
271: g0[0] = 1.0 / mu;
272: }
274: /* Quadratic MMS Solution
275: 2D:
277: u = x^2 + y^2
278: v = 2 x^2 - 2xy
279: p = x + y - 1
280: f = <1 - 4 mu, 1 - 4 mu>
282: so that
284: e(u) = (grad u + grad u^T) = / 4x 4x \
285: \ 4x -4x /
286: div mu e(u) - \nabla p + f = mu <4, 4> - <1, 1> + <1 - 4 mu, 1 - 4 mu> = 0
287: \nabla \cdot u = 2x - 2x = 0
289: 3D:
291: u = 2 x^2 + y^2 + z^2
292: v = 2 x^2 - 2xy
293: w = 2 x^2 - 2xz
294: p = x + y + z - 3/2
295: f = <1 - 8 mu, 1 - 4 mu, 1 - 4 mu>
297: so that
299: e(u) = (grad u + grad u^T) = / 8x 4x 4x \
300: | 4x -4x 0 |
301: \ 4x 0 -4x /
302: div mu e(u) - \nabla p + f = mu <8, 4, 4> - <1, 1, 1> + <1 - 8 mu, 1 - 4 mu, 1 - 4 mu> = 0
303: \nabla \cdot u = 4x - 2x - 2x = 0
304: */
305: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
306: {
307: PetscInt c;
309: u[0] = (dim - 1) * PetscSqr(x[0]);
310: for (c = 1; c < Nc; ++c) {
311: u[0] += PetscSqr(x[c]);
312: u[c] = 2.0 * PetscSqr(x[0]) - 2.0 * x[0] * x[c];
313: }
314: return PETSC_SUCCESS;
315: }
317: static PetscErrorCode quadratic_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
318: {
319: PetscInt d;
321: u[0] = -0.5 * dim;
322: for (d = 0; d < dim; ++d) u[0] += x[d];
323: return PETSC_SUCCESS;
324: }
326: static void f0_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
327: {
328: const PetscReal mu = PetscRealPart(constants[0]);
329: PetscInt d;
331: f0[0] = (dim - 1) * 4.0 * mu - 1.0;
332: for (d = 1; d < dim; ++d) f0[d] = 4.0 * mu - 1.0;
333: }
335: /* Trigonometric MMS Solution
336: 2D:
338: u = sin(pi x) + sin(pi y)
339: v = -pi cos(pi x) y
340: p = sin(2 pi x) + sin(2 pi y)
341: f = <2pi cos(2 pi x) + mu pi^2 sin(pi x) + mu pi^2 sin(pi y), 2pi cos(2 pi y) - mu pi^3 cos(pi x) y>
343: so that
345: e(u) = (grad u + grad u^T) = / 2pi cos(pi x) pi cos(pi y) + pi^2 sin(pi x) y \
346: \ pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) /
347: div mu e(u) - \nabla p + f = mu <-pi^2 sin(pi x) - pi^2 sin(pi y), pi^3 cos(pi x) y> - <2pi cos(2 pi x), 2pi cos(2 pi y)> + <f_x, f_y> = 0
348: \nabla \cdot u = pi cos(pi x) - pi cos(pi x) = 0
350: 3D:
352: u = 2 sin(pi x) + sin(pi y) + sin(pi z)
353: v = -pi cos(pi x) y
354: w = -pi cos(pi x) z
355: p = sin(2 pi x) + sin(2 pi y) + sin(2 pi z)
356: f = <2pi cos(2 pi x) + mu 2pi^2 sin(pi x) + mu pi^2 sin(pi y) + mu pi^2 sin(pi z), 2pi cos(2 pi y) - mu pi^3 cos(pi x) y, 2pi cos(2 pi z) - mu pi^3 cos(pi x) z>
358: so that
360: e(u) = (grad u + grad u^T) = / 4pi cos(pi x) pi cos(pi y) + pi^2 sin(pi x) y pi cos(pi z) + pi^2 sin(pi x) z \
361: | pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) 0 |
362: \ pi cos(pi z) + pi^2 sin(pi x) z 0 -2pi cos(pi x) /
363: div mu e(u) - \nabla p + f = mu <-2pi^2 sin(pi x) - pi^2 sin(pi y) - pi^2 sin(pi z), pi^3 cos(pi x) y, pi^3 cos(pi x) z> - <2pi cos(2 pi x), 2pi cos(2 pi y), 2pi cos(2 pi z)> + <f_x, f_y, f_z> = 0
364: \nabla \cdot u = 2 pi cos(pi x) - pi cos(pi x) - pi cos(pi x) = 0
365: */
366: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
367: {
368: PetscInt c;
370: u[0] = (dim - 1) * PetscSinReal(PETSC_PI * x[0]);
371: for (c = 1; c < Nc; ++c) {
372: u[0] += PetscSinReal(PETSC_PI * x[c]);
373: u[c] = -PETSC_PI * PetscCosReal(PETSC_PI * x[0]) * x[c];
374: }
375: return PETSC_SUCCESS;
376: }
378: static PetscErrorCode trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
379: {
380: PetscInt d;
382: for (d = 0, u[0] = 0.0; d < dim; ++d) u[0] += PetscSinReal(2.0 * PETSC_PI * x[d]);
383: return PETSC_SUCCESS;
384: }
386: static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
387: {
388: const PetscReal mu = PetscRealPart(constants[0]);
389: PetscInt d;
391: f0[0] = -2.0 * PETSC_PI * PetscCosReal(2.0 * PETSC_PI * x[0]) - (dim - 1) * mu * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x[0]);
392: for (d = 1; d < dim; ++d) {
393: f0[0] -= mu * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x[d]);
394: f0[d] = -2.0 * PETSC_PI * PetscCosReal(2.0 * PETSC_PI * x[d]) + mu * PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * x[0]) * x[d];
395: }
396: }
398: /* Inline helpers for computing exact velocity in void boundary kernels */
399: static inline void ExactVelocityQuadratic(PetscInt dim, const PetscReal x[], PetscScalar g[])
400: {
401: PetscInt c;
403: g[0] = (dim - 1) * x[0] * x[0];
404: for (c = 1; c < dim; ++c) {
405: g[0] += x[c] * x[c];
406: g[c] = 2.0 * x[0] * x[0] - 2.0 * x[0] * x[c];
407: }
408: }
410: static inline void ExactVelocityTrig(PetscInt dim, const PetscReal x[], PetscScalar g[])
411: {
412: PetscInt c;
414: g[0] = (dim - 1) * PetscSinReal(PETSC_PI * x[0]);
415: for (c = 1; c < dim; ++c) {
416: g[0] += PetscSinReal(PETSC_PI * x[c]);
417: g[c] = -PETSC_PI * PetscCosReal(PETSC_PI * x[0]) * x[c];
418: }
419: }
421: /* Nitsche boundary residual kernels for velocity (field 0)
422: f0_bd_u[c] = -mu * sum_d (u_x[c*dim+d] + u_x[d*dim+c]) * n[d] (consistency: stress flux from IBP)
423: + p * n[c] (pressure flux from IBP)
424: + penalty * (u[c] - g[c]) (penalty) */
425: static void f0_bd_nitsche_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
426: {
427: const PetscReal mu = PetscRealPart(constants[0]);
428: const PetscReal penalty = PetscRealPart(constants[1]);
429: PetscScalar g[3];
430: PetscInt c, d;
432: ExactVelocityQuadratic(dim, x, g);
433: for (c = 0; c < dim; ++c) {
434: f0[c] = penalty * (u[c] - g[c]) + u[uOff[1]] * n[c];
435: for (d = 0; d < dim; ++d) f0[c] -= mu * (u_x[c * dim + d] + u_x[d * dim + c]) * n[d];
436: }
437: }
439: static void f0_bd_nitsche_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
440: {
441: const PetscReal mu = PetscRealPart(constants[0]);
442: const PetscReal penalty = PetscRealPart(constants[1]);
443: PetscScalar g[3];
444: PetscInt c, d;
446: ExactVelocityTrig(dim, x, g);
447: for (c = 0; c < dim; ++c) {
448: f0[c] = penalty * (u[c] - g[c]) + u[uOff[1]] * n[c];
449: for (d = 0; d < dim; ++d) f0[c] -= mu * (u_x[c * dim + d] + u_x[d * dim + c]) * n[d];
450: }
451: }
453: /* f1_bd_u[c*dim+d] = -mu * (n[d]*(u[c]-g[c]) + n[c]*(u[d]-g[d])) (symmetry / adjoint consistency) */
454: static void f1_bd_nitsche_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
455: {
456: const PetscReal mu = PetscRealPart(constants[0]);
457: PetscScalar g[3];
458: PetscInt c, d;
460: ExactVelocityQuadratic(dim, x, g);
461: for (c = 0; c < dim; ++c)
462: for (d = 0; d < dim; ++d) f1[c * dim + d] = -mu * (n[d] * (u[c] - g[c]) + n[c] * (u[d] - g[d]));
463: }
465: static void f1_bd_nitsche_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
466: {
467: const PetscReal mu = PetscRealPart(constants[0]);
468: PetscScalar g[3];
469: PetscInt c, d;
471: ExactVelocityTrig(dim, x, g);
472: for (c = 0; c < dim; ++c)
473: for (d = 0; d < dim; ++d) f1[c * dim + d] = -mu * (n[d] * (u[c] - g[c]) + n[c] * (u[d] - g[d]));
474: }
476: /* Nitsche boundary residual kernels for pressure (field 1)
477: f0_bd_p = sum_d n[d] * (u[d] - g[d]) (continuity equation boundary correction) */
478: static void f0_bd_nitsche_quadratic_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
479: {
480: PetscScalar g[3];
481: PetscInt d;
483: ExactVelocityQuadratic(dim, x, g);
484: f0[0] = 0.0;
485: for (d = 0; d < dim; ++d) f0[0] += n[d] * (u[d] - g[d]);
486: }
488: static void f0_bd_nitsche_trig_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
489: {
490: PetscScalar g[3];
491: PetscInt d;
493: ExactVelocityTrig(dim, x, g);
494: f0[0] = 0.0;
495: for (d = 0; d < dim; ++d) f0[0] += n[d] * (u[d] - g[d]);
496: }
498: /* Nitsche boundary Jacobian kernels (solution-independent)
499: g0_bd_uu[c*Nc+d] = delta(c,d) * penalty (penalty Jacobian) */
500: static void g0_bd_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
501: {
502: const PetscReal penalty = PetscRealPart(constants[1]);
503: const PetscInt Nc = uOff[1] - uOff[0];
504: PetscInt c;
506: for (c = 0; c < Nc; ++c) g0[c * Nc + c] = penalty;
507: }
509: /* g1_bd_uu[(c*Nc+d)*dim+e] = -mu * (delta(c,d)*n[e] + delta(c,e)*n[d]) (consistency Jacobian: df0/du_x) */
510: static void g1_bd_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
511: {
512: const PetscReal mu = PetscRealPart(constants[0]);
513: const PetscInt Nc = uOff[1] - uOff[0];
514: PetscInt c, d, e;
516: for (c = 0; c < Nc; ++c)
517: for (d = 0; d < Nc; ++d)
518: for (e = 0; e < dim; ++e) g1[(c * Nc + d) * dim + e] = -mu * ((c == d ? 1.0 : 0.0) * n[e] + (c == e ? 1.0 : 0.0) * n[d]);
519: }
521: /* g2_bd_uu[(c*Nc+d)*dim+e] = -mu * (n[e]*delta(c,d) + n[c]*delta(e,d)) (symmetry Jacobian: df1/du) */
522: static void g2_bd_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
523: {
524: const PetscReal mu = PetscRealPart(constants[0]);
525: const PetscInt Nc = uOff[1] - uOff[0];
526: PetscInt c, d, e;
528: for (c = 0; c < Nc; ++c)
529: for (d = 0; d < Nc; ++d)
530: for (e = 0; e < dim; ++e) g2[(c * Nc + d) * dim + e] = -mu * (n[e] * (c == d ? 1.0 : 0.0) + n[c] * (e == d ? 1.0 : 0.0));
531: }
533: /* g0_bd_up[c*1+0] = n[c] (velocity-pressure coupling: df0_u/dp) */
534: static void g0_bd_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
535: {
536: for (PetscInt c = 0; c < dim; ++c) g0[c] = n[c];
537: }
539: /* g0_bd_pu[0*Nc+d] = n[d] (pressure-velocity coupling: df0_p/du) */
540: static void g0_bd_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
541: {
542: for (PetscInt d = 0; d < dim; ++d) g0[d] = n[d];
543: }
545: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
546: {
547: PetscInt sol, bc;
549: PetscFunctionBeginUser;
550: options->sol = SOL_QUADRATIC;
551: options->bc = BC_ESSENTIAL;
552: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
553: sol = options->sol;
554: PetscCall(PetscOptionsEList("-sol", "The MMS solution", "ex62.c", SolTypes, PETSC_STATIC_ARRAY_LENGTH(SolTypes) - 3, SolTypes[options->sol], &sol, NULL));
555: options->sol = (SolType)sol;
556: bc = options->bc;
557: PetscCall(PetscOptionsEList("-bc", "The boundary condition type", "ex62.c", BCTypes, PETSC_STATIC_ARRAY_LENGTH(BCTypes) - 3, BCTypes[options->bc], &bc, NULL));
558: options->bc = (BCType)bc;
559: PetscOptionsEnd();
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
564: {
565: PetscFunctionBeginUser;
566: PetscCall(DMCreate(comm, dm));
567: PetscCall(DMSetType(*dm, DMPLEX));
568: PetscCall(DMSetFromOptions(*dm));
569: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
574: {
575: Parameter *p;
577: PetscFunctionBeginUser;
578: /* setup PETSc parameter bag */
579: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx->bag));
580: PetscCall(PetscBagGetData(ctx->bag, &p));
581: PetscCall(PetscBagSetName(ctx->bag, "par", "Stokes Parameters"));
582: PetscCall(PetscBagRegisterScalar(ctx->bag, &p->mu, 1.0, "mu", "Dynamic Shear Viscosity, Pa s"));
583: PetscCall(PetscBagRegisterScalar(ctx->bag, &p->eta, 100.0, "eta", "Nitsche penalty parameter (dimensionless)"));
584: PetscCall(PetscBagSetFromOptions(ctx->bag));
585: {
586: PetscViewer viewer;
587: PetscViewerFormat format;
588: PetscBool flg;
590: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
591: if (flg) {
592: PetscCall(PetscViewerPushFormat(viewer, format));
593: PetscCall(PetscBagView(ctx->bag, viewer));
594: PetscCall(PetscViewerFlush(viewer));
595: PetscCall(PetscViewerPopFormat(viewer));
596: PetscCall(PetscViewerDestroy(&viewer));
597: }
598: }
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode SetupEqn(DM dm, AppCtx *user)
603: {
604: PetscErrorCode (*exactFuncs[2])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
605: void (*f0_bd_u)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
606: void (*f1_bd_u)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
607: void (*f0_bd_p)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
608: PetscDS ds;
609: DMLabel label;
610: const PetscInt id = 1;
612: PetscFunctionBeginUser;
613: PetscCall(DMGetDS(dm, &ds));
614: switch (user->sol) {
615: case SOL_QUADRATIC:
616: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_u, f1_u));
617: exactFuncs[0] = quadratic_u;
618: exactFuncs[1] = quadratic_p;
619: f0_bd_u = f0_bd_nitsche_quadratic_u;
620: f1_bd_u = f1_bd_nitsche_quadratic_u;
621: f0_bd_p = f0_bd_nitsche_quadratic_p;
622: break;
623: case SOL_TRIG:
624: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
625: exactFuncs[0] = trig_u;
626: exactFuncs[1] = trig_p;
627: f0_bd_u = f0_bd_nitsche_trig_u;
628: f1_bd_u = f1_bd_nitsche_trig_u;
629: f0_bd_p = f0_bd_nitsche_trig_p;
630: break;
631: default:
632: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", SolTypes[PetscMin(user->sol, SOL_UNKNOWN)], user->sol);
633: }
634: PetscCall(PetscDSSetResidual(ds, 1, f0_p, NULL));
635: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
636: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL));
637: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL));
638: PetscCall(PetscDSSetJacobianPreconditioner(ds, 0, 0, NULL, NULL, NULL, g3_uu));
639: PetscCall(PetscDSSetJacobianPreconditioner(ds, 1, 1, g0_pp, NULL, NULL, NULL));
641: PetscCall(PetscDSSetExactSolution(ds, 0, exactFuncs[0], user));
642: PetscCall(PetscDSSetExactSolution(ds, 1, exactFuncs[1], user));
644: PetscCall(DMGetLabel(dm, "marker", &label));
645: switch (user->bc) {
646: case BC_ESSENTIAL:
647: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, user, NULL));
648: break;
649: case BC_NITSCHE: {
650: PetscWeakForm wf;
651: DMLabel faceSetsLabel;
652: IS valueIS;
653: const PetscInt *faceSetValues;
654: PetscInt numValues, bd, i;
656: PetscCall(DMGetLabel(dm, "Face Sets", &faceSetsLabel));
657: PetscCall(DMLabelGetNumValues(faceSetsLabel, &numValues));
658: PetscCall(DMLabelGetValueIS(faceSetsLabel, &valueIS));
659: PetscCall(ISGetIndices(valueIS, &faceSetValues));
661: /* Velocity boundary: natural BC with Nitsche terms on all boundary faces */
662: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "wall", faceSetsLabel, numValues, faceSetValues, 0, 0, NULL, NULL, NULL, user, &bd));
663: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
664: for (i = 0; i < numValues; ++i) {
665: /* Velocity residual (field 0): f0 and f1 */
666: PetscCall(PetscWeakFormSetIndexBdResidual(wf, faceSetsLabel, faceSetValues[i], 0, 0, 0, f0_bd_u, 0, f1_bd_u));
667: /* Velocity-velocity Jacobian (field 0, field 0): g0 (penalty), g1 (consistency), g2 (symmetry) */
668: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, faceSetsLabel, faceSetValues[i], 0, 0, 0, 0, g0_bd_uu, 0, g1_bd_uu, 0, g2_bd_uu, 0, NULL));
669: /* Velocity-pressure Jacobian (field 0, field 1): g0 (pressure coupling) */
670: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, faceSetsLabel, faceSetValues[i], 0, 1, 0, 0, g0_bd_up, 0, NULL, 0, NULL, 0, NULL));
671: }
673: /* Pressure boundary: natural BC for continuity equation correction */
674: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "wall_pres", faceSetsLabel, numValues, faceSetValues, 1, 0, NULL, NULL, NULL, user, &bd));
675: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
676: for (i = 0; i < numValues; ++i) {
677: /* Pressure residual (field 1): f0 */
678: PetscCall(PetscWeakFormSetIndexBdResidual(wf, faceSetsLabel, faceSetValues[i], 1, 0, 0, f0_bd_p, 0, NULL));
679: /* Pressure-velocity Jacobian (field 1, field 0): g0 */
680: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, faceSetsLabel, faceSetValues[i], 1, 0, 0, 0, g0_bd_pu, 0, NULL, 0, NULL, 0, NULL));
681: }
682: PetscCall(ISRestoreIndices(valueIS, &faceSetValues));
683: PetscCall(ISDestroy(&valueIS));
684: } break;
685: default:
686: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unsupported BC type: %s (%d)", BCTypes[PetscMin(user->bc, BC_UNKNOWN)], user->bc);
687: }
689: /* Make constant values available to pointwise functions */
690: {
691: Parameter *param;
692: PetscScalar constants[2];
694: PetscCall(PetscBagGetData(user->bag, ¶m));
695: constants[0] = param->mu; /* dynamic shear viscosity, Pa s */
696: constants[1] = 0.0; /* Nitsche penalty (set below if needed) */
697: if (user->bc == BC_NITSCHE) {
698: /* Compute cell size h from mesh */
699: PetscInt dim, cStart;
700: PetscReal vol, h;
702: PetscCall(DMGetDimension(dm, &dim));
703: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
704: PetscCall(DMPlexComputeCellGeometryFVM(dm, cStart, &vol, NULL, NULL));
705: h = PetscPowReal(vol, 1.0 / dim);
706: constants[1] = PetscRealPart(param->eta) * PetscRealPart(param->mu) / h;
707: }
708: PetscCall(PetscDSSetConstants(ds, 2, constants));
709: }
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
714: {
715: for (PetscInt c = 0; c < Nc; ++c) u[c] = 0.0;
716: return PETSC_SUCCESS;
717: }
718: static PetscErrorCode one(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
719: {
720: for (PetscInt c = 0; c < Nc; ++c) u[c] = 1.0;
721: return PETSC_SUCCESS;
722: }
724: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
725: {
726: Vec vec;
727: PetscErrorCode (*funcs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx) = {zero, one};
729: PetscFunctionBeginUser;
730: PetscCheck(origField == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Field %" PetscInt_FMT " should be 1 for pressure", origField);
731: funcs[field] = one;
732: {
733: PetscDS ds;
734: PetscCall(DMGetDS(dm, &ds));
735: PetscCall(PetscObjectViewFromOptions((PetscObject)ds, NULL, "-ds_view"));
736: }
737: PetscCall(DMCreateGlobalVector(dm, &vec));
738: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
739: PetscCall(VecNormalize(vec, NULL));
740: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullspace));
741: PetscCall(VecDestroy(&vec));
742: /* New style for field null spaces */
743: {
744: PetscObject pressure;
745: MatNullSpace nullspacePres;
747: PetscCall(DMGetField(dm, field, NULL, &pressure));
748: PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
749: PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
750: PetscCall(MatNullSpaceDestroy(&nullspacePres));
751: }
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: static PetscErrorCode SetupProblem(DM dm, PetscErrorCode (*setupEqn)(DM, AppCtx *), AppCtx *user)
756: {
757: DM cdm = dm;
758: PetscQuadrature q = NULL;
759: PetscBool simplex;
760: PetscInt dim, Nf = 2, f, Nc[2];
761: const char *name[2] = {"velocity", "pressure"};
762: const char *prefix[2] = {"vel_", "pres_"};
764: PetscFunctionBegin;
765: PetscCall(DMGetDimension(dm, &dim));
766: PetscCall(DMPlexIsSimplex(dm, &simplex));
767: Nc[0] = dim;
768: Nc[1] = 1;
769: for (f = 0; f < Nf; ++f) {
770: PetscFE fe;
772: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, prefix[f], -1, &fe));
773: PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
774: if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
775: PetscCall(PetscFESetQuadrature(fe, q));
776: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
777: PetscCall(PetscFEDestroy(&fe));
778: }
779: PetscCall(DMCreateDS(dm));
780: PetscCall((*setupEqn)(dm, user));
781: while (cdm) {
782: PetscCall(DMCopyDisc(dm, cdm));
783: PetscCall(DMSetNullSpaceConstructor(cdm, 1, CreatePressureNullSpace));
784: PetscCall(DMGetCoarseDM(cdm, &cdm));
785: }
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: int main(int argc, char **argv)
790: {
791: SNES snes;
792: DM dm;
793: Vec u;
794: AppCtx user;
795: PetscBool testPatchFacetResidual = PETSC_FALSE;
797: PetscFunctionBeginUser;
798: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
799: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_patch_facet_residual", &testPatchFacetResidual, NULL));
800: if (testPatchFacetResidual) PetscCall(TestPatchFacetResidual());
801: else {
802: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
803: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
804: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
805: PetscCall(SNESSetDM(snes, dm));
806: PetscCall(DMSetApplicationContext(dm, &user));
808: PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
809: PetscCall(SetupProblem(dm, SetupEqn, &user));
810: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
812: PetscCall(DMCreateGlobalVector(dm, &u));
813: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
814: PetscCall(SNESSetFromOptions(snes));
815: PetscCall(DMSNESCheckFromOptions(snes, u));
816: PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
817: {
818: Mat J;
819: MatNullSpace sp;
821: PetscCall(SNESSetUp(snes));
822: PetscCall(CreatePressureNullSpace(dm, 1, 1, &sp));
823: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
824: PetscCall(MatSetNullSpace(J, sp));
825: PetscCall(MatNullSpaceDestroy(&sp));
826: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
827: PetscCall(MatViewFromOptions(J, NULL, "-J_view"));
828: }
829: PetscCall(SNESSolve(snes, NULL, u));
831: PetscCall(VecDestroy(&u));
832: PetscCall(SNESDestroy(&snes));
833: PetscCall(DMDestroy(&dm));
834: PetscCall(PetscBagDestroy(&user.bag));
835: }
836: PetscCall(PetscFinalize());
837: return 0;
838: }
839: /*TEST
841: test:
842: suffix: 2d_p2_p1_check
843: requires: triangle
844: args: -sol quadratic -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
846: test:
847: suffix: 2d_p2_p1_check_parallel
848: nsize: {{2 3 5}}
849: requires: triangle
850: args: -sol quadratic -dm_refine 2 -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
852: test:
853: suffix: 3d_p2_p1_check
854: requires: ctetgen
855: args: -sol quadratic -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
857: test:
858: suffix: 3d_p2_p1_check_parallel
859: nsize: {{2 3 5}}
860: requires: ctetgen
861: args: -sol quadratic -dm_refine 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
863: test:
864: suffix: 2d_p2_p1_conv
865: requires: triangle
866: # Using -dm_refine 3 gives L_2 convergence rate: [3.0, 2.1]
867: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
868: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
869: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
870: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
872: test:
873: suffix: 2d_p2_p1_conv_gamg
874: requires: triangle
875: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 \
876: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
877: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_explicit_operator_mat_type aij -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_coarse_pc_type svd
879: test:
880: suffix: 3d_p2_p1_conv
881: requires: ctetgen !single
882: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.8]
883: args: -sol trig -dm_plex_dim 3 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 \
884: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
885: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
886: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
888: test:
889: suffix: 2d_q2_q1_check
890: args: -sol quadratic -dm_plex_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
892: test:
893: suffix: 3d_q2_q1_check
894: args: -sol quadratic -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
896: test:
897: suffix: 2d_q2_q1_conv
898: # Using -dm_refine 3 -convest_num_refine 1 gives L_2 convergence rate: [3.0, 2.1]
899: args: -sol trig -dm_plex_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -ksp_error_if_not_converged \
900: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
901: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
902: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
904: test:
905: suffix: 3d_q2_q1_conv
906: requires: !single
907: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.4]
908: args: -sol trig -dm_plex_simplex 0 -dm_plex_dim 3 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 \
909: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
910: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
911: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
913: test:
914: suffix: 2d_p3_p2_check
915: requires: triangle
916: args: -sol quadratic -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
918: test:
919: suffix: 3d_p3_p2_check
920: requires: ctetgen !single
921: args: -sol quadratic -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
923: test:
924: suffix: 2d_p3_p2_conv
925: requires: triangle
926: # Using -dm_refine 2 gives L_2 convergence rate: [3.8, 3.0]
927: args: -sol trig -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
928: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
929: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
930: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
932: test:
933: suffix: 3d_p3_p2_conv
934: requires: ctetgen long_runtime
935: # Using -dm_refine 1 -convest_num_refine 2 gives L_2 convergence rate: [3.6, 3.9]
936: args: -sol trig -dm_plex_dim 3 -dm_refine 1 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 \
937: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
938: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
939: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
941: test:
942: suffix: 2d_q1_p0_conv
943: requires: !single
944: # Using -dm_refine 3 gives L_2 convergence rate: [1.9, 1.0]
945: args: -sol quadratic -dm_plex_simplex 0 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 2 \
946: -ksp_atol 1e-10 -petscds_jac_pre 0 \
947: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
948: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_explicit_operator_mat_type aij -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_levels_pc_type jacobi -fieldsplit_pressure_mg_coarse_pc_type svd -fieldsplit_pressure_pc_gamg_aggressive_coarsening 0
950: test:
951: suffix: 3d_q1_p0_conv
952: requires: !single
953: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [1.7, 1.0]
954: args: -sol quadratic -dm_plex_simplex 0 -dm_plex_dim 3 -dm_refine 1 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 1 \
955: -ksp_atol 1e-10 -petscds_jac_pre 0 \
956: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
957: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_explicit_operator_mat_type aij -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_mg_levels_pc_type jacobi -fieldsplit_pressure_mg_coarse_pc_type svd -fieldsplit_pressure_pc_gamg_aggressive_coarsening 0
959: # Stokes preconditioners
960: # Block diagonal \begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix}
961: test:
962: suffix: 2d_p2_p1_block_diagonal
963: requires: triangle
964: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
965: -snes_error_if_not_converged \
966: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-4 -ksp_error_if_not_converged \
967: -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
968: output_file: output/empty.out
969: # Block triangular \begin{pmatrix} A & B \\ 0 & I \end{pmatrix}
970: test:
971: suffix: 2d_p2_p1_block_triangular
972: requires: triangle
973: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
974: -snes_error_if_not_converged \
975: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
976: -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
977: output_file: output/empty.out
978: # Diagonal Schur complement \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix}
979: test:
980: suffix: 2d_p2_p1_schur_diagonal
981: requires: triangle
982: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
983: -snes_error_if_not_converged \
984: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
985: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -pc_fieldsplit_off_diag_use_amat \
986: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
987: output_file: output/empty.out
988: # Upper triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
989: test:
990: suffix: 2d_p2_p1_schur_upper
991: requires: triangle
992: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001 \
993: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
994: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -pc_fieldsplit_off_diag_use_amat \
995: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
996: # Lower triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
997: test:
998: suffix: 2d_p2_p1_schur_lower
999: requires: triangle
1000: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1001: -snes_error_if_not_converged \
1002: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
1003: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -pc_fieldsplit_off_diag_use_amat \
1004: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1005: output_file: output/empty.out
1006: # Full Schur complement \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix} \begin{pmatrix} I & A^{-1} B \\ 0 & I \end{pmatrix}
1007: test:
1008: suffix: 2d_p2_p1_schur_full
1009: requires: triangle
1010: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1011: -snes_error_if_not_converged \
1012: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
1013: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_off_diag_use_amat \
1014: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1015: output_file: output/empty.out
1016: # Full Schur + Velocity GMG
1017: test:
1018: suffix: 2d_p2_p1_gmg_vcycle
1019: TODO: broken (requires subDMs hooks)
1020: requires: triangle
1021: args: -sol quadratic -dm_refine_hierarchy 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1022: -ksp_type fgmres -ksp_atol 1e-9 -snes_error_if_not_converged -pc_use_amat \
1023: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_off_diag_use_amat \
1024: -fieldsplit_velocity_pc_type mg -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type gamg -fieldsplit_pressure_pc_gamg_esteig_ksp_max_it 10 -fieldsplit_pressure_mg_levels_pc_type sor -fieldsplit_pressure_mg_coarse_pc_type svd
1025: # SIMPLE \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & B^T diag(A)^{-1} B \end{pmatrix} \begin{pmatrix} I & diag(A)^{-1} B \\ 0 & I \end{pmatrix}
1026: test:
1027: suffix: 2d_p2_p1_simple
1028: requires: triangle
1029: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1030: -snes_error_if_not_converged \
1031: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
1032: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
1033: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
1034: -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi
1035: output_file: output/empty.out
1036: # FETI-DP solvers (these solvers are quite inefficient, they are here to exercise the code)
1037: test:
1038: suffix: 2d_p2_p1_fetidp
1039: requires: triangle mumps
1040: nsize: 5
1041: args: -sol quadratic -dm_refine 2 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1042: -snes_error_if_not_converged \
1043: -ksp_type fetidp -ksp_rtol 1.0e-8 \
1044: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1045: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 200 -fetidp_fieldsplit_p_pc_type none \
1046: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type mumps -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type mumps -fetidp_fieldsplit_lag_ksp_type preonly
1047: output_file: output/empty.out
1048: test:
1049: suffix: 2d_q2_q1_fetidp
1050: requires: mumps
1051: nsize: 5
1052: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1053: -ksp_type fetidp -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
1054: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1055: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 200 -fetidp_fieldsplit_p_pc_type none \
1056: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type mumps -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type mumps -fetidp_fieldsplit_lag_ksp_type preonly
1057: output_file: output/empty.out
1058: test:
1059: suffix: 3d_p2_p1_fetidp
1060: requires: ctetgen mumps suitesparse
1061: nsize: 5
1062: args: -sol quadratic -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1063: -snes_error_if_not_converged \
1064: -ksp_type fetidp -ksp_rtol 1.0e-9 \
1065: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1066: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 1000 -fetidp_fieldsplit_p_pc_type none \
1067: -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_benign_trick -fetidp_bddc_pc_bddc_deluxe_singlemat \
1068: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
1069: -fetidp_bddelta_pc_factor_mat_solver_type umfpack -fetidp_fieldsplit_lag_ksp_type preonly -fetidp_bddc_sub_schurs_mat_solver_type mumps -fetidp_bddc_sub_schurs_mat_mumps_icntl_14 100000 \
1070: -fetidp_bddelta_pc_factor_mat_ordering_type external \
1071: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
1072: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
1073: output_file: output/empty.out
1074: test:
1075: suffix: 3d_q2_q1_fetidp
1076: requires: suitesparse
1077: nsize: 5
1078: args: -sol quadratic -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1079: -snes_error_if_not_converged \
1080: -ksp_type fetidp -ksp_rtol 1.0e-8 \
1081: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1082: -fetidp_fieldsplit_p_ksp_max_it 1 -fetidp_fieldsplit_p_ksp_type richardson -fetidp_fieldsplit_p_ksp_richardson_scale 2000 -fetidp_fieldsplit_p_pc_type none \
1083: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
1084: -fetidp_bddc_pc_bddc_symmetric -fetidp_fieldsplit_lag_ksp_type preonly \
1085: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
1086: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
1087: output_file: output/empty.out
1088: # BDDC solvers (these solvers are quite inefficient, they are here to exercise the code)
1089: test:
1090: suffix: 2d_p2_p1_bddc
1091: nsize: 2
1092: requires: triangle !single
1093: args: -sol quadratic -dm_plex_box_faces 2,2,2 -dm_refine 1 -dm_mat_type is -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1094: -snes_error_if_not_converged \
1095: -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
1096: -pc_type bddc -pc_bddc_corner_selection -pc_bddc_dirichlet_pc_type svd -pc_bddc_neumann_pc_type svd -pc_bddc_coarse_redundant_pc_type svd
1097: output_file: output/empty.out
1098: # Vanka
1099: test:
1100: suffix: patch_facet_residual
1101: args: -test_patch_facet_residual
1102: output_file: output/empty.out
1103: test:
1104: suffix: 2d_q1_p0_vanka
1105: output_file: output/empty.out
1106: requires: double !complex
1107: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1108: -snes_rtol 1.0e-4 \
1109: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1110: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
1111: -sub_ksp_type preonly -sub_pc_type lu
1112: test:
1113: suffix: 2d_q1_p0_vanka_denseinv
1114: output_file: output/empty.out
1115: requires: double !complex
1116: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1117: -snes_rtol 1.0e-4 \
1118: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1119: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
1120: -pc_patch_dense_inverse -pc_patch_sub_mat_type seqdense
1121: # Vanka smoother
1122: test:
1123: suffix: 2d_q1_p0_gmg_vanka
1124: output_file: output/empty.out
1125: requires: double !complex
1126: args: -sol quadratic -dm_plex_simplex 0 -dm_refine_hierarchy 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1127: -snes_rtol 1.0e-4 \
1128: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1129: -pc_type mg \
1130: -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 30 \
1131: -mg_levels_pc_type patch -mg_levels_pc_patch_partition_of_unity 0 -mg_levels_pc_patch_construct_codim 0 -mg_levels_pc_patch_construct_type vanka \
1132: -mg_levels_sub_ksp_type preonly -mg_levels_sub_pc_type lu \
1133: -mg_coarse_pc_type svd
1134: # Nitsche BC consistency check
1135: test:
1136: suffix: 2d_q2_q1_nitsche_check
1137: requires: double !complex
1138: args: -sol quadratic -bc nitsche -dm_plex_simplex 0 -dm_refine 1 \
1139: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
1140: # Nitsche BC + Vanka
1141: test:
1142: suffix: 2d_q1_p0_nitsche_vanka
1143: output_file: output/empty.out
1144: requires: double !complex
1145: args: -sol quadratic -bc nitsche -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1146: -snes_rtol 1.0e-4 \
1147: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1148: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
1149: -sub_ksp_type preonly -sub_pc_type lu
1150: # Nitsche BC + GMRES (sanity check that Nitsche formulation solves correctly)
1151: test:
1152: suffix: 2d_q2_q1_nitsche_gmres
1153: output_file: output/empty.out
1154: requires: double !complex
1155: args: -sol quadratic -bc nitsche -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1156: -snes_error_if_not_converged \
1157: -ksp_type gmres -ksp_rtol 1e-12 -pc_type jacobi
1159: TEST*/