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_", 0};
45: typedef enum {
46: BC_ESSENTIAL,
47: BC_NITSCHE,
48: BC_UNKNOWN
49: } BCType;
50: const char *BCTypes[] = {"essential", "nitsche", "unknown", "BCType", "BC_", 0};
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 bs[1] = {1}, nodesPerCell[1] = {nodesPerCellValue}, subspaceOffsets[2] = {0, 0};
145: DM dms[1];
146: SNES snes;
147: Vec x, f, rhs;
148: PatchFacetTestCtx test = {0};
150: PetscFunctionBeginUser;
151: PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm));
152: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
153: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
154: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
155: PetscCheck(cStart == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected cell points to start at zero, got %" PetscInt_FMT, cStart);
156: numCells = cEnd - cStart;
157: test.dofsPerCell = nodesPerCellValue;
158: test.numPatchPoints = vEnd - vStart;
159: PetscCall(PetscMalloc1(test.numPatchPoints, &test.patchPoints));
160: for (PetscInt vertex = vStart; vertex < vEnd; ++vertex) test.patchPoints[vertex - vStart] = vertex;
162: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, §ion));
163: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
164: for (PetscInt vertex = vStart; vertex < vEnd; ++vertex) PetscCall(PetscSectionSetDof(section, vertex, 1));
165: PetscCall(PetscSectionSetUp(section));
166: PetscCall(DMSetLocalSection(dm, section));
167: PetscCall(PetscSectionGetStorageSize(section, &numDofs));
168: subspaceOffsets[1] = numDofs;
170: PetscCall(PetscMalloc1(numCells * nodesPerCellValue, &cellNodeMap));
171: for (cell = cStart; cell < cEnd; ++cell) {
172: PetscInt numVertices = 0;
174: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
175: for (PetscInt i = 0; i < closureSize * 2; i += 2) {
176: const PetscInt point = closure[i];
178: if (point >= vStart && point < vEnd) {
179: PetscInt offset;
181: PetscCall(PetscSectionGetOffset(section, point, &offset));
182: PetscCheck(numVertices < nodesPerCellValue, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Found too many vertices in cell");
183: cellNodeMap[cell * nodesPerCellValue + numVertices++] = offset;
184: }
185: }
186: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
187: PetscCheck(numVertices == nodesPerCellValue, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected %" PetscInt_FMT " vertices in cell, got %" PetscInt_FMT, nodesPerCellValue, numVertices);
188: }
190: dms[0] = dm;
191: cellNodeMaps[0] = cellNodeMap;
192: PetscCall(VecCreateSeq(PETSC_COMM_SELF, numDofs, &x));
193: PetscCall(VecDuplicate(x, &f));
194: PetscCall(VecDuplicate(x, &rhs));
195: PetscCall(VecSet(rhs, 0.0));
196: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
197: PetscCall(SNESSetType(snes, SNESPATCH));
198: PetscCall(SNESSetDM(snes, dm));
199: PetscCall(SNESSetFunction(snes, f, TestPatchOuterFunction, NULL));
200: PetscCall(SNESPatchSetConstructType(snes, PC_PATCH_USER, TestPatchConstruct, &test));
201: PetscCall(SNESPatchSetDiscretisationInfo(snes, 1, dms, bs, nodesPerCell, cellNodeMaps, subspaceOffsets, 0, NULL, 0, NULL));
202: PetscCall(SNESPatchSetComputeFunction(snes, TestPatchComputeFunction, &test));
203: PetscCall(SNESPatchSetComputeOperator(snes, TestPatchComputeOperator, NULL));
204: PetscCall(SNESSetTolerances(snes, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1, PETSC_CURRENT));
205: PetscCall(PetscOptionsSetValue(NULL, "-sub_snes_max_it", "1"));
206: PetscCall(SNESSetFromOptions(snes));
207: PetscCall(SNESSolve(snes, rhs, x));
208: PetscCheck(test.numInteriorFacetCalls > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior facet callback was not called");
209: PetscCheck(test.numExteriorFacetCalls > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Exterior facet callback was not called");
210: PetscCall(SNESDestroy(&snes));
211: PetscCall(VecDestroy(&rhs));
212: PetscCall(VecDestroy(&f));
213: PetscCall(VecDestroy(&x));
214: PetscCall(PetscFree(cellNodeMap));
215: PetscCall(PetscSectionDestroy(§ion));
216: PetscCall(DMDestroy(&dm));
217: PetscCall(PetscFree(test.patchPoints));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: 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[])
222: {
223: const PetscReal mu = PetscRealPart(constants[0]);
224: const PetscInt Nc = uOff[1] - uOff[0];
225: PetscInt c, d;
227: for (c = 0; c < Nc; ++c) {
228: for (d = 0; d < dim; ++d) f1[c * dim + d] = mu * (u_x[c * dim + d] + u_x[d * dim + c]);
229: f1[c * dim + c] -= u[uOff[1]];
230: }
231: }
233: 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[])
234: {
235: PetscInt d;
236: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] -= u_x[d * dim + d];
237: }
239: 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[])
240: {
241: PetscInt d;
242: for (d = 0; d < dim; ++d) g1[d * dim + d] = -1.0; /* < q, -\nabla\cdot u > */
243: }
245: 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[])
246: {
247: PetscInt d;
248: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* -< \nabla\cdot v, p > */
249: }
251: 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[])
252: {
253: const PetscReal mu = PetscRealPart(constants[0]);
254: const PetscInt Nc = uOff[1] - uOff[0];
255: PetscInt c, d;
257: for (c = 0; c < Nc; ++c) {
258: for (d = 0; d < dim; ++d) {
259: g3[((c * Nc + c) * dim + d) * dim + d] += mu; /* < \nabla v, \nabla u > */
260: g3[((c * Nc + d) * dim + d) * dim + c] += mu; /* < \nabla v, {\nabla u}^T > */
261: }
262: }
263: }
265: 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[])
266: {
267: const PetscReal mu = PetscRealPart(constants[0]);
269: g0[0] = 1.0 / mu;
270: }
272: /* Quadratic MMS Solution
273: 2D:
275: u = x^2 + y^2
276: v = 2 x^2 - 2xy
277: p = x + y - 1
278: f = <1 - 4 mu, 1 - 4 mu>
280: so that
282: e(u) = (grad u + grad u^T) = / 4x 4x \
283: \ 4x -4x /
284: div mu e(u) - \nabla p + f = mu <4, 4> - <1, 1> + <1 - 4 mu, 1 - 4 mu> = 0
285: \nabla \cdot u = 2x - 2x = 0
287: 3D:
289: u = 2 x^2 + y^2 + z^2
290: v = 2 x^2 - 2xy
291: w = 2 x^2 - 2xz
292: p = x + y + z - 3/2
293: f = <1 - 8 mu, 1 - 4 mu, 1 - 4 mu>
295: so that
297: e(u) = (grad u + grad u^T) = / 8x 4x 4x \
298: | 4x -4x 0 |
299: \ 4x 0 -4x /
300: div mu e(u) - \nabla p + f = mu <8, 4, 4> - <1, 1, 1> + <1 - 8 mu, 1 - 4 mu, 1 - 4 mu> = 0
301: \nabla \cdot u = 4x - 2x - 2x = 0
302: */
303: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
304: {
305: PetscInt c;
307: u[0] = (dim - 1) * PetscSqr(x[0]);
308: for (c = 1; c < Nc; ++c) {
309: u[0] += PetscSqr(x[c]);
310: u[c] = 2.0 * PetscSqr(x[0]) - 2.0 * x[0] * x[c];
311: }
312: return PETSC_SUCCESS;
313: }
315: static PetscErrorCode quadratic_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
316: {
317: PetscInt d;
319: u[0] = -0.5 * dim;
320: for (d = 0; d < dim; ++d) u[0] += x[d];
321: return PETSC_SUCCESS;
322: }
324: 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[])
325: {
326: const PetscReal mu = PetscRealPart(constants[0]);
327: PetscInt d;
329: f0[0] = (dim - 1) * 4.0 * mu - 1.0;
330: for (d = 1; d < dim; ++d) f0[d] = 4.0 * mu - 1.0;
331: }
333: /* Trigonometric MMS Solution
334: 2D:
336: u = sin(pi x) + sin(pi y)
337: v = -pi cos(pi x) y
338: p = sin(2 pi x) + sin(2 pi y)
339: 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>
341: so that
343: e(u) = (grad u + grad u^T) = / 2pi cos(pi x) pi cos(pi y) + pi^2 sin(pi x) y \
344: \ pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) /
345: 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
346: \nabla \cdot u = pi cos(pi x) - pi cos(pi x) = 0
348: 3D:
350: u = 2 sin(pi x) + sin(pi y) + sin(pi z)
351: v = -pi cos(pi x) y
352: w = -pi cos(pi x) z
353: p = sin(2 pi x) + sin(2 pi y) + sin(2 pi z)
354: 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>
356: so that
358: 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 \
359: | pi cos(pi y) + pi^2 sin(pi x) y -2pi cos(pi x) 0 |
360: \ pi cos(pi z) + pi^2 sin(pi x) z 0 -2pi cos(pi x) /
361: 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
362: \nabla \cdot u = 2 pi cos(pi x) - pi cos(pi x) - pi cos(pi x) = 0
363: */
364: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
365: {
366: PetscInt c;
368: u[0] = (dim - 1) * PetscSinReal(PETSC_PI * x[0]);
369: for (c = 1; c < Nc; ++c) {
370: u[0] += PetscSinReal(PETSC_PI * x[c]);
371: u[c] = -PETSC_PI * PetscCosReal(PETSC_PI * x[0]) * x[c];
372: }
373: return PETSC_SUCCESS;
374: }
376: static PetscErrorCode trig_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
377: {
378: PetscInt d;
380: for (d = 0, u[0] = 0.0; d < dim; ++d) u[0] += PetscSinReal(2.0 * PETSC_PI * x[d]);
381: return PETSC_SUCCESS;
382: }
384: 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[])
385: {
386: const PetscReal mu = PetscRealPart(constants[0]);
387: PetscInt d;
389: f0[0] = -2.0 * PETSC_PI * PetscCosReal(2.0 * PETSC_PI * x[0]) - (dim - 1) * mu * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x[0]);
390: for (d = 1; d < dim; ++d) {
391: f0[0] -= mu * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * x[d]);
392: 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];
393: }
394: }
396: /* Inline helpers for computing exact velocity in void boundary kernels */
397: static inline void ExactVelocityQuadratic(PetscInt dim, const PetscReal x[], PetscScalar g[])
398: {
399: PetscInt c;
401: g[0] = (dim - 1) * x[0] * x[0];
402: for (c = 1; c < dim; ++c) {
403: g[0] += x[c] * x[c];
404: g[c] = 2.0 * x[0] * x[0] - 2.0 * x[0] * x[c];
405: }
406: }
408: static inline void ExactVelocityTrig(PetscInt dim, const PetscReal x[], PetscScalar g[])
409: {
410: PetscInt c;
412: g[0] = (dim - 1) * PetscSinReal(PETSC_PI * x[0]);
413: for (c = 1; c < dim; ++c) {
414: g[0] += PetscSinReal(PETSC_PI * x[c]);
415: g[c] = -PETSC_PI * PetscCosReal(PETSC_PI * x[0]) * x[c];
416: }
417: }
419: /* Nitsche boundary residual kernels for velocity (field 0)
420: f0_bd_u[c] = -mu * sum_d (u_x[c*dim+d] + u_x[d*dim+c]) * n[d] (consistency: stress flux from IBP)
421: + p * n[c] (pressure flux from IBP)
422: + penalty * (u[c] - g[c]) (penalty) */
423: 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[])
424: {
425: const PetscReal mu = PetscRealPart(constants[0]);
426: const PetscReal penalty = PetscRealPart(constants[1]);
427: PetscScalar g[3];
428: PetscInt c, d;
430: ExactVelocityQuadratic(dim, x, g);
431: for (c = 0; c < dim; ++c) {
432: f0[c] = penalty * (u[c] - g[c]) + u[uOff[1]] * n[c];
433: for (d = 0; d < dim; ++d) f0[c] -= mu * (u_x[c * dim + d] + u_x[d * dim + c]) * n[d];
434: }
435: }
437: 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[])
438: {
439: const PetscReal mu = PetscRealPart(constants[0]);
440: const PetscReal penalty = PetscRealPart(constants[1]);
441: PetscScalar g[3];
442: PetscInt c, d;
444: ExactVelocityTrig(dim, x, g);
445: for (c = 0; c < dim; ++c) {
446: f0[c] = penalty * (u[c] - g[c]) + u[uOff[1]] * n[c];
447: for (d = 0; d < dim; ++d) f0[c] -= mu * (u_x[c * dim + d] + u_x[d * dim + c]) * n[d];
448: }
449: }
451: /* f1_bd_u[c*dim+d] = -mu * (n[d]*(u[c]-g[c]) + n[c]*(u[d]-g[d])) (symmetry / adjoint consistency) */
452: 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[])
453: {
454: const PetscReal mu = PetscRealPart(constants[0]);
455: PetscScalar g[3];
456: PetscInt c, d;
458: ExactVelocityQuadratic(dim, x, g);
459: for (c = 0; c < dim; ++c)
460: for (d = 0; d < dim; ++d) f1[c * dim + d] = -mu * (n[d] * (u[c] - g[c]) + n[c] * (u[d] - g[d]));
461: }
463: 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[])
464: {
465: const PetscReal mu = PetscRealPart(constants[0]);
466: PetscScalar g[3];
467: PetscInt c, d;
469: ExactVelocityTrig(dim, x, g);
470: for (c = 0; c < dim; ++c)
471: for (d = 0; d < dim; ++d) f1[c * dim + d] = -mu * (n[d] * (u[c] - g[c]) + n[c] * (u[d] - g[d]));
472: }
474: /* Nitsche boundary residual kernels for pressure (field 1)
475: f0_bd_p = sum_d n[d] * (u[d] - g[d]) (continuity equation boundary correction) */
476: 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[])
477: {
478: PetscScalar g[3];
479: PetscInt d;
481: ExactVelocityQuadratic(dim, x, g);
482: f0[0] = 0.0;
483: for (d = 0; d < dim; ++d) f0[0] += n[d] * (u[d] - g[d]);
484: }
486: 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[])
487: {
488: PetscScalar g[3];
489: PetscInt d;
491: ExactVelocityTrig(dim, x, g);
492: f0[0] = 0.0;
493: for (d = 0; d < dim; ++d) f0[0] += n[d] * (u[d] - g[d]);
494: }
496: /* Nitsche boundary Jacobian kernels (solution-independent)
497: g0_bd_uu[c*Nc+d] = delta(c,d) * penalty (penalty Jacobian) */
498: 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[])
499: {
500: const PetscReal penalty = PetscRealPart(constants[1]);
501: const PetscInt Nc = uOff[1] - uOff[0];
502: PetscInt c;
504: for (c = 0; c < Nc; ++c) g0[c * Nc + c] = penalty;
505: }
507: /* g1_bd_uu[(c*Nc+d)*dim+e] = -mu * (delta(c,d)*n[e] + delta(c,e)*n[d]) (consistency Jacobian: df0/du_x) */
508: 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[])
509: {
510: const PetscReal mu = PetscRealPart(constants[0]);
511: const PetscInt Nc = uOff[1] - uOff[0];
512: PetscInt c, d, e;
514: for (c = 0; c < Nc; ++c)
515: for (d = 0; d < Nc; ++d)
516: 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]);
517: }
519: /* g2_bd_uu[(c*Nc+d)*dim+e] = -mu * (n[e]*delta(c,d) + n[c]*delta(e,d)) (symmetry Jacobian: df1/du) */
520: 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[])
521: {
522: const PetscReal mu = PetscRealPart(constants[0]);
523: const PetscInt Nc = uOff[1] - uOff[0];
524: PetscInt c, d, e;
526: for (c = 0; c < Nc; ++c)
527: for (d = 0; d < Nc; ++d)
528: 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));
529: }
531: /* g0_bd_up[c*1+0] = n[c] (velocity-pressure coupling: df0_u/dp) */
532: 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[])
533: {
534: for (PetscInt c = 0; c < dim; ++c) g0[c] = n[c];
535: }
537: /* g0_bd_pu[0*Nc+d] = n[d] (pressure-velocity coupling: df0_p/du) */
538: 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[])
539: {
540: for (PetscInt d = 0; d < dim; ++d) g0[d] = n[d];
541: }
543: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
544: {
545: PetscInt sol, bc;
547: PetscFunctionBeginUser;
548: options->sol = SOL_QUADRATIC;
549: options->bc = BC_ESSENTIAL;
550: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
551: sol = options->sol;
552: PetscCall(PetscOptionsEList("-sol", "The MMS solution", "ex62.c", SolTypes, PETSC_STATIC_ARRAY_LENGTH(SolTypes) - 3, SolTypes[options->sol], &sol, NULL));
553: options->sol = (SolType)sol;
554: bc = options->bc;
555: PetscCall(PetscOptionsEList("-bc", "The boundary condition type", "ex62.c", BCTypes, PETSC_STATIC_ARRAY_LENGTH(BCTypes) - 3, BCTypes[options->bc], &bc, NULL));
556: options->bc = (BCType)bc;
557: PetscOptionsEnd();
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
562: {
563: PetscFunctionBeginUser;
564: PetscCall(DMCreate(comm, dm));
565: PetscCall(DMSetType(*dm, DMPLEX));
566: PetscCall(DMSetFromOptions(*dm));
567: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
572: {
573: Parameter *p;
575: PetscFunctionBeginUser;
576: /* setup PETSc parameter bag */
577: PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx->bag));
578: PetscCall(PetscBagGetData(ctx->bag, &p));
579: PetscCall(PetscBagSetName(ctx->bag, "par", "Stokes Parameters"));
580: PetscCall(PetscBagRegisterScalar(ctx->bag, &p->mu, 1.0, "mu", "Dynamic Shear Viscosity, Pa s"));
581: PetscCall(PetscBagRegisterScalar(ctx->bag, &p->eta, 100.0, "eta", "Nitsche penalty parameter (dimensionless)"));
582: PetscCall(PetscBagSetFromOptions(ctx->bag));
583: {
584: PetscViewer viewer;
585: PetscViewerFormat format;
586: PetscBool flg;
588: PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
589: if (flg) {
590: PetscCall(PetscViewerPushFormat(viewer, format));
591: PetscCall(PetscBagView(ctx->bag, viewer));
592: PetscCall(PetscViewerFlush(viewer));
593: PetscCall(PetscViewerPopFormat(viewer));
594: PetscCall(PetscViewerDestroy(&viewer));
595: }
596: }
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: static PetscErrorCode SetupEqn(DM dm, AppCtx *user)
601: {
602: PetscErrorCode (*exactFuncs[2])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
603: 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[]);
604: 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[]);
605: 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[]);
606: PetscDS ds;
607: DMLabel label;
608: const PetscInt id = 1;
610: PetscFunctionBeginUser;
611: PetscCall(DMGetDS(dm, &ds));
612: switch (user->sol) {
613: case SOL_QUADRATIC:
614: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_u, f1_u));
615: exactFuncs[0] = quadratic_u;
616: exactFuncs[1] = quadratic_p;
617: f0_bd_u = f0_bd_nitsche_quadratic_u;
618: f1_bd_u = f1_bd_nitsche_quadratic_u;
619: f0_bd_p = f0_bd_nitsche_quadratic_p;
620: break;
621: case SOL_TRIG:
622: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
623: exactFuncs[0] = trig_u;
624: exactFuncs[1] = trig_p;
625: f0_bd_u = f0_bd_nitsche_trig_u;
626: f1_bd_u = f1_bd_nitsche_trig_u;
627: f0_bd_p = f0_bd_nitsche_trig_p;
628: break;
629: default:
630: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", SolTypes[PetscMin(user->sol, SOL_UNKNOWN)], user->sol);
631: }
632: PetscCall(PetscDSSetResidual(ds, 1, f0_p, NULL));
633: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
634: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL));
635: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL));
636: PetscCall(PetscDSSetJacobianPreconditioner(ds, 0, 0, NULL, NULL, NULL, g3_uu));
637: PetscCall(PetscDSSetJacobianPreconditioner(ds, 1, 1, g0_pp, NULL, NULL, NULL));
639: PetscCall(PetscDSSetExactSolution(ds, 0, exactFuncs[0], user));
640: PetscCall(PetscDSSetExactSolution(ds, 1, exactFuncs[1], user));
642: PetscCall(DMGetLabel(dm, "marker", &label));
643: switch (user->bc) {
644: case BC_ESSENTIAL:
645: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, user, NULL));
646: break;
647: case BC_NITSCHE: {
648: PetscWeakForm wf;
649: DMLabel faceSetsLabel;
650: IS valueIS;
651: const PetscInt *faceSetValues;
652: PetscInt numValues, bd, i;
654: PetscCall(DMGetLabel(dm, "Face Sets", &faceSetsLabel));
655: PetscCall(DMLabelGetNumValues(faceSetsLabel, &numValues));
656: PetscCall(DMLabelGetValueIS(faceSetsLabel, &valueIS));
657: PetscCall(ISGetIndices(valueIS, &faceSetValues));
659: /* Velocity boundary: natural BC with Nitsche terms on all boundary faces */
660: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "wall", faceSetsLabel, numValues, faceSetValues, 0, 0, NULL, NULL, NULL, user, &bd));
661: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
662: for (i = 0; i < numValues; ++i) {
663: /* Velocity residual (field 0): f0 and f1 */
664: PetscCall(PetscWeakFormSetIndexBdResidual(wf, faceSetsLabel, faceSetValues[i], 0, 0, 0, f0_bd_u, 0, f1_bd_u));
665: /* Velocity-velocity Jacobian (field 0, field 0): g0 (penalty), g1 (consistency), g2 (symmetry) */
666: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, faceSetsLabel, faceSetValues[i], 0, 0, 0, 0, g0_bd_uu, 0, g1_bd_uu, 0, g2_bd_uu, 0, NULL));
667: /* Velocity-pressure Jacobian (field 0, field 1): g0 (pressure coupling) */
668: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, faceSetsLabel, faceSetValues[i], 0, 1, 0, 0, g0_bd_up, 0, NULL, 0, NULL, 0, NULL));
669: }
671: /* Pressure boundary: natural BC for continuity equation correction */
672: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "wall_pres", faceSetsLabel, numValues, faceSetValues, 1, 0, NULL, NULL, NULL, user, &bd));
673: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
674: for (i = 0; i < numValues; ++i) {
675: /* Pressure residual (field 1): f0 */
676: PetscCall(PetscWeakFormSetIndexBdResidual(wf, faceSetsLabel, faceSetValues[i], 1, 0, 0, f0_bd_p, 0, NULL));
677: /* Pressure-velocity Jacobian (field 1, field 0): g0 */
678: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, faceSetsLabel, faceSetValues[i], 1, 0, 0, 0, g0_bd_pu, 0, NULL, 0, NULL, 0, NULL));
679: }
680: PetscCall(ISRestoreIndices(valueIS, &faceSetValues));
681: PetscCall(ISDestroy(&valueIS));
682: } break;
683: default:
684: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unsupported BC type: %s (%d)", BCTypes[PetscMin(user->bc, BC_UNKNOWN)], user->bc);
685: }
687: /* Make constant values available to pointwise functions */
688: {
689: Parameter *param;
690: PetscScalar constants[2];
692: PetscCall(PetscBagGetData(user->bag, ¶m));
693: constants[0] = param->mu; /* dynamic shear viscosity, Pa s */
694: constants[1] = 0.0; /* Nitsche penalty (set below if needed) */
695: if (user->bc == BC_NITSCHE) {
696: /* Compute cell size h from mesh */
697: PetscInt dim, cStart;
698: PetscReal vol, h;
700: PetscCall(DMGetDimension(dm, &dim));
701: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
702: PetscCall(DMPlexComputeCellGeometryFVM(dm, cStart, &vol, NULL, NULL));
703: h = PetscPowReal(vol, 1.0 / dim);
704: constants[1] = PetscRealPart(param->eta) * PetscRealPart(param->mu) / h;
705: }
706: PetscCall(PetscDSSetConstants(ds, 2, constants));
707: }
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
712: {
713: for (PetscInt c = 0; c < Nc; ++c) u[c] = 0.0;
714: return PETSC_SUCCESS;
715: }
716: static PetscErrorCode one(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
717: {
718: for (PetscInt c = 0; c < Nc; ++c) u[c] = 1.0;
719: return PETSC_SUCCESS;
720: }
722: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
723: {
724: Vec vec;
725: PetscErrorCode (*funcs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx) = {zero, one};
727: PetscFunctionBeginUser;
728: PetscCheck(origField == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Field %" PetscInt_FMT " should be 1 for pressure", origField);
729: funcs[field] = one;
730: {
731: PetscDS ds;
732: PetscCall(DMGetDS(dm, &ds));
733: PetscCall(PetscObjectViewFromOptions((PetscObject)ds, NULL, "-ds_view"));
734: }
735: PetscCall(DMCreateGlobalVector(dm, &vec));
736: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
737: PetscCall(VecNormalize(vec, NULL));
738: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullspace));
739: PetscCall(VecDestroy(&vec));
740: /* New style for field null spaces */
741: {
742: PetscObject pressure;
743: MatNullSpace nullspacePres;
745: PetscCall(DMGetField(dm, field, NULL, &pressure));
746: PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
747: PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
748: PetscCall(MatNullSpaceDestroy(&nullspacePres));
749: }
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: static PetscErrorCode SetupProblem(DM dm, PetscErrorCode (*setupEqn)(DM, AppCtx *), AppCtx *user)
754: {
755: DM cdm = dm;
756: PetscQuadrature q = NULL;
757: PetscBool simplex;
758: PetscInt dim, Nf = 2, f, Nc[2];
759: const char *name[2] = {"velocity", "pressure"};
760: const char *prefix[2] = {"vel_", "pres_"};
762: PetscFunctionBegin;
763: PetscCall(DMGetDimension(dm, &dim));
764: PetscCall(DMPlexIsSimplex(dm, &simplex));
765: Nc[0] = dim;
766: Nc[1] = 1;
767: for (f = 0; f < Nf; ++f) {
768: PetscFE fe;
770: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, prefix[f], -1, &fe));
771: PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
772: if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
773: PetscCall(PetscFESetQuadrature(fe, q));
774: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
775: PetscCall(PetscFEDestroy(&fe));
776: }
777: PetscCall(DMCreateDS(dm));
778: PetscCall((*setupEqn)(dm, user));
779: while (cdm) {
780: PetscCall(DMCopyDisc(dm, cdm));
781: PetscCall(DMSetNullSpaceConstructor(cdm, 1, CreatePressureNullSpace));
782: PetscCall(DMGetCoarseDM(cdm, &cdm));
783: }
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: int main(int argc, char **argv)
788: {
789: SNES snes;
790: DM dm;
791: Vec u;
792: AppCtx user;
793: PetscBool testPatchFacetResidual = PETSC_FALSE;
795: PetscFunctionBeginUser;
796: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
797: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_patch_facet_residual", &testPatchFacetResidual, NULL));
798: if (testPatchFacetResidual) PetscCall(TestPatchFacetResidual());
799: else {
800: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
801: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
802: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
803: PetscCall(SNESSetDM(snes, dm));
804: PetscCall(DMSetApplicationContext(dm, &user));
806: PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
807: PetscCall(SetupProblem(dm, SetupEqn, &user));
808: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
810: PetscCall(DMCreateGlobalVector(dm, &u));
811: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
812: PetscCall(SNESSetFromOptions(snes));
813: PetscCall(DMSNESCheckFromOptions(snes, u));
814: PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
815: {
816: Mat J;
817: MatNullSpace sp;
819: PetscCall(SNESSetUp(snes));
820: PetscCall(CreatePressureNullSpace(dm, 1, 1, &sp));
821: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
822: PetscCall(MatSetNullSpace(J, sp));
823: PetscCall(MatNullSpaceDestroy(&sp));
824: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
825: PetscCall(MatViewFromOptions(J, NULL, "-J_view"));
826: }
827: PetscCall(SNESSolve(snes, NULL, u));
829: PetscCall(VecDestroy(&u));
830: PetscCall(SNESDestroy(&snes));
831: PetscCall(DMDestroy(&dm));
832: PetscCall(PetscBagDestroy(&user.bag));
833: }
834: PetscCall(PetscFinalize());
835: return 0;
836: }
837: /*TEST
839: test:
840: suffix: 2d_p2_p1_check
841: requires: triangle
842: args: -sol quadratic -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
844: test:
845: suffix: 2d_p2_p1_check_parallel
846: nsize: {{2 3 5}}
847: requires: triangle
848: args: -sol quadratic -dm_refine 2 -petscpartitioner_type simple -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
850: test:
851: suffix: 3d_p2_p1_check
852: requires: ctetgen
853: 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
855: test:
856: suffix: 3d_p2_p1_check_parallel
857: nsize: {{2 3 5}}
858: requires: ctetgen
859: 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
861: test:
862: suffix: 2d_p2_p1_conv
863: requires: triangle
864: # Using -dm_refine 3 gives L_2 convergence rate: [3.0, 2.1]
865: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
866: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
867: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
868: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
870: test:
871: suffix: 2d_p2_p1_conv_gamg
872: requires: triangle
873: args: -sol trig -vel_petscspace_degree 2 -pres_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 \
874: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
875: -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
877: test:
878: suffix: 3d_p2_p1_conv
879: requires: ctetgen !single
880: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.8]
881: 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 \
882: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
883: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
884: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
886: test:
887: suffix: 2d_q2_q1_check
888: args: -sol quadratic -dm_plex_simplex 0 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
890: test:
891: suffix: 3d_q2_q1_check
892: 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
894: test:
895: suffix: 2d_q2_q1_conv
896: # Using -dm_refine 3 -convest_num_refine 1 gives L_2 convergence rate: [3.0, 2.1]
897: 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 \
898: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
899: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
900: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
902: test:
903: suffix: 3d_q2_q1_conv
904: requires: !single
905: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [2.8, 2.4]
906: 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 \
907: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
908: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
909: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
911: test:
912: suffix: 2d_p3_p2_check
913: requires: triangle
914: args: -sol quadratic -vel_petscspace_degree 3 -pres_petscspace_degree 2 -dmsnes_check 0.0001
916: test:
917: suffix: 3d_p3_p2_check
918: requires: ctetgen !single
919: 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
921: test:
922: suffix: 2d_p3_p2_conv
923: requires: triangle
924: # Using -dm_refine 2 gives L_2 convergence rate: [3.8, 3.0]
925: args: -sol trig -vel_petscspace_degree 3 -pres_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 -ksp_error_if_not_converged \
926: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
927: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
928: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
930: test:
931: suffix: 3d_p3_p2_conv
932: requires: ctetgen long_runtime
933: # Using -dm_refine 1 -convest_num_refine 2 gives L_2 convergence rate: [3.6, 3.9]
934: 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 \
935: -ksp_atol 1e-10 -ksp_error_if_not_converged -pc_use_amat \
936: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition a11 -pc_fieldsplit_off_diag_use_amat \
937: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type lu
939: test:
940: suffix: 2d_q1_p0_conv
941: requires: !single
942: # Using -dm_refine 3 gives L_2 convergence rate: [1.9, 1.0]
943: args: -sol quadratic -dm_plex_simplex 0 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -snes_convergence_estimate -convest_num_refine 2 \
944: -ksp_atol 1e-10 -petscds_jac_pre 0 \
945: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
946: -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
948: test:
949: suffix: 3d_q1_p0_conv
950: requires: !single
951: # Using -dm_refine 2 -convest_num_refine 2 gives L_2 convergence rate: [1.7, 1.0]
952: 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 \
953: -ksp_atol 1e-10 -petscds_jac_pre 0 \
954: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full \
955: -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
957: # Stokes preconditioners
958: # Block diagonal \begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix}
959: test:
960: suffix: 2d_p2_p1_block_diagonal
961: requires: triangle
962: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
963: -snes_error_if_not_converged \
964: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-4 -ksp_error_if_not_converged \
965: -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
966: output_file: output/empty.out
967: # Block triangular \begin{pmatrix} A & B \\ 0 & I \end{pmatrix}
968: test:
969: suffix: 2d_p2_p1_block_triangular
970: requires: triangle
971: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
972: -snes_error_if_not_converged \
973: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
974: -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi
975: output_file: output/empty.out
976: # Diagonal Schur complement \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix}
977: test:
978: suffix: 2d_p2_p1_schur_diagonal
979: requires: triangle
980: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
981: -snes_error_if_not_converged \
982: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
983: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -pc_fieldsplit_off_diag_use_amat \
984: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
985: output_file: output/empty.out
986: # Upper triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
987: test:
988: suffix: 2d_p2_p1_schur_upper
989: requires: triangle
990: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001 \
991: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
992: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -pc_fieldsplit_off_diag_use_amat \
993: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
994: # Lower triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
995: test:
996: suffix: 2d_p2_p1_schur_lower
997: requires: triangle
998: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
999: -snes_error_if_not_converged \
1000: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
1001: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -pc_fieldsplit_off_diag_use_amat \
1002: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1003: output_file: output/empty.out
1004: # 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}
1005: test:
1006: suffix: 2d_p2_p1_schur_full
1007: requires: triangle
1008: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1009: -snes_error_if_not_converged \
1010: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged -pc_use_amat \
1011: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_off_diag_use_amat \
1012: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1013: output_file: output/empty.out
1014: # Full Schur + Velocity GMG
1015: test:
1016: suffix: 2d_p2_p1_gmg_vcycle
1017: TODO: broken (requires subDMs hooks)
1018: requires: triangle
1019: args: -sol quadratic -dm_refine_hierarchy 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1020: -ksp_type fgmres -ksp_atol 1e-9 -snes_error_if_not_converged -pc_use_amat \
1021: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_off_diag_use_amat \
1022: -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
1023: # 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}
1024: test:
1025: suffix: 2d_p2_p1_simple
1026: requires: triangle
1027: args: -sol quadratic -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -petscds_jac_pre 0 \
1028: -snes_error_if_not_converged \
1029: -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
1030: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
1031: -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
1032: -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi
1033: output_file: output/empty.out
1034: # FETI-DP solvers (these solvers are quite inefficient, they are here to exercise the code)
1035: test:
1036: suffix: 2d_p2_p1_fetidp
1037: requires: triangle mumps
1038: nsize: 5
1039: 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 \
1040: -snes_error_if_not_converged \
1041: -ksp_type fetidp -ksp_rtol 1.0e-8 \
1042: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1043: -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 \
1044: -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
1045: output_file: output/empty.out
1046: test:
1047: suffix: 2d_q2_q1_fetidp
1048: requires: mumps
1049: nsize: 5
1050: 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 \
1051: -ksp_type fetidp -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
1052: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1053: -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 \
1054: -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
1055: output_file: output/empty.out
1056: test:
1057: suffix: 3d_p2_p1_fetidp
1058: requires: ctetgen mumps suitesparse
1059: nsize: 5
1060: 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 \
1061: -snes_error_if_not_converged \
1062: -ksp_type fetidp -ksp_rtol 1.0e-9 \
1063: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1064: -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 \
1065: -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_benign_trick -fetidp_bddc_pc_bddc_deluxe_singlemat \
1066: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
1067: -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 \
1068: -fetidp_bddelta_pc_factor_mat_ordering_type external \
1069: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
1070: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
1071: output_file: output/empty.out
1072: test:
1073: suffix: 3d_q2_q1_fetidp
1074: requires: suitesparse
1075: nsize: 5
1076: 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 \
1077: -snes_error_if_not_converged \
1078: -ksp_type fetidp -ksp_rtol 1.0e-8 \
1079: -ksp_fetidp_saddlepoint -fetidp_ksp_type cg \
1080: -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 \
1081: -fetidp_pc_discrete_harmonic -fetidp_harmonic_pc_factor_mat_solver_type petsc -fetidp_harmonic_pc_type cholesky \
1082: -fetidp_bddc_pc_bddc_symmetric -fetidp_fieldsplit_lag_ksp_type preonly \
1083: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_solver_type umfpack -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_solver_type umfpack \
1084: -fetidp_bddc_pc_bddc_dirichlet_pc_factor_mat_ordering_type external -fetidp_bddc_pc_bddc_neumann_pc_factor_mat_ordering_type external
1085: output_file: output/empty.out
1086: # BDDC solvers (these solvers are quite inefficient, they are here to exercise the code)
1087: test:
1088: suffix: 2d_p2_p1_bddc
1089: nsize: 2
1090: requires: triangle !single
1091: 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 \
1092: -snes_error_if_not_converged \
1093: -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-8 -ksp_error_if_not_converged \
1094: -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
1095: output_file: output/empty.out
1096: # Vanka
1097: test:
1098: suffix: patch_facet_residual
1099: args: -test_patch_facet_residual
1100: output_file: output/empty.out
1101: test:
1102: suffix: 2d_q1_p0_vanka
1103: output_file: output/empty.out
1104: requires: double !complex
1105: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1106: -snes_rtol 1.0e-4 \
1107: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1108: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
1109: -sub_ksp_type preonly -sub_pc_type lu
1110: test:
1111: suffix: 2d_q1_p0_vanka_denseinv
1112: output_file: output/empty.out
1113: requires: double !complex
1114: args: -sol quadratic -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1115: -snes_rtol 1.0e-4 \
1116: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1117: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
1118: -pc_patch_dense_inverse -pc_patch_sub_mat_type seqdense
1119: # Vanka smoother
1120: test:
1121: suffix: 2d_q1_p0_gmg_vanka
1122: output_file: output/empty.out
1123: requires: double !complex
1124: args: -sol quadratic -dm_plex_simplex 0 -dm_refine_hierarchy 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1125: -snes_rtol 1.0e-4 \
1126: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1127: -pc_type mg \
1128: -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 30 \
1129: -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 \
1130: -mg_levels_sub_ksp_type preonly -mg_levels_sub_pc_type lu \
1131: -mg_coarse_pc_type svd
1132: # Nitsche BC consistency check
1133: test:
1134: suffix: 2d_q2_q1_nitsche_check
1135: requires: double !complex
1136: args: -sol quadratic -bc nitsche -dm_plex_simplex 0 -dm_refine 1 \
1137: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -dmsnes_check 0.0001
1138: # Nitsche BC + Vanka
1139: test:
1140: suffix: 2d_q1_p0_nitsche_vanka
1141: output_file: output/empty.out
1142: requires: double !complex
1143: args: -sol quadratic -bc nitsche -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 1 -pres_petscspace_degree 0 -petscds_jac_pre 0 \
1144: -snes_rtol 1.0e-4 \
1145: -ksp_type fgmres -ksp_atol 1e-5 -ksp_error_if_not_converged \
1146: -pc_type patch -pc_patch_partition_of_unity 0 -pc_patch_construct_codim 0 -pc_patch_construct_type vanka \
1147: -sub_ksp_type preonly -sub_pc_type lu
1148: # Nitsche BC + GMRES (sanity check that Nitsche formulation solves correctly)
1149: test:
1150: suffix: 2d_q2_q1_nitsche_gmres
1151: output_file: output/empty.out
1152: requires: double !complex
1153: args: -sol quadratic -bc nitsche -dm_plex_simplex 0 -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
1154: -snes_error_if_not_converged \
1155: -ksp_type gmres -ksp_rtol 1e-12 -pc_type jacobi
1157: TEST*/