Actual source code: ex73.c
1: static char help[] = "Tests for Gauss' Law\n\n";
3: /* We want to check the weak version of Gauss' Law, namely that
5: \int_\Omega v div q - \int_\Gamma v (q \cdot n) = 0
7: */
9: #include <petscdmplex.h>
10: #include <petscsnes.h>
11: #include <petscds.h>
13: typedef struct {
14: PetscInt degree; // The degree of the discretization
15: PetscBool divFree; // True if the solution is divergence-free
16: } AppCtx;
18: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
19: {
20: u[0] = 0.0;
21: return PETSC_SUCCESS;
22: }
24: // div = 0
25: static void solenoidal_2d(PetscInt n, const PetscReal x[], PetscScalar u[])
26: {
27: u[0] = PetscPowRealInt(x[0], n) * PetscPowRealInt(x[1], n - 1);
28: u[1] = -PetscPowRealInt(x[0], n - 1) * PetscPowRealInt(x[1], n);
29: }
30: // div = 0
31: static void solenoidal_3d(PetscInt n, const PetscReal x[], PetscScalar u[])
32: {
33: u[0] = PetscPowRealInt(x[0], n) * PetscPowRealInt(x[1], n - 1) * PetscPowRealInt(x[2], n - 1);
34: u[1] = -2. * PetscPowRealInt(x[0], n - 1) * PetscPowRealInt(x[1], n) * PetscPowRealInt(x[2], n - 1);
35: u[2] = PetscPowRealInt(x[0], n - 1) * PetscPowRealInt(x[1], n - 1) * PetscPowRealInt(x[2], n);
36: }
38: static PetscErrorCode solenoidal_totaldeg_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
39: {
40: const PetscInt deg = *(PetscInt *)ctx;
41: const PetscInt n = deg / 2 + deg % 2;
43: solenoidal_2d(n, x, u);
44: return PETSC_SUCCESS;
45: }
47: static PetscErrorCode solenoidal_totaldeg_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
48: {
49: const PetscInt deg = *(PetscInt *)ctx;
50: const PetscInt n = deg / 3 + (deg % 3 ? 1 : 0);
52: solenoidal_3d(n, x, u);
53: return PETSC_SUCCESS;
54: }
56: // This is in P_n^{-}
57: static PetscErrorCode source_totaldeg(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
58: {
59: const PetscInt n = *(PetscInt *)ctx;
61: for (PetscInt d = 0; d < dim; ++d) u[d] = PetscPowRealInt(x[d], n + 1);
62: return PETSC_SUCCESS;
63: }
65: static void identity(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[])
66: {
67: const PetscInt deg = (PetscInt)PetscRealPart(constants[0]);
68: PetscScalar p[3];
70: if (dim == 2) PetscCallVoid(solenoidal_totaldeg_2d(dim, t, x, uOff[1] - uOff[0], p, (void *)°));
71: else PetscCallVoid(solenoidal_totaldeg_3d(dim, t, x, uOff[1] - uOff[0], p, (void *)°));
72: for (PetscInt c = 0; c < dim; ++c) f0[c] = -u[c] + p[c];
73: }
75: static void zero_bd(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[])
76: {
77: for (PetscInt d = 0; d < dim; ++d) f0[0] = 0.;
78: }
80: static void flux(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[])
81: {
82: for (PetscInt d = 0; d < dim; ++d) f0[0] -= u[d] * n[d];
83: }
85: static void divergence(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[])
86: {
87: for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d];
88: }
90: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
91: {
92: PetscFunctionBegin;
93: PetscCall(DMCreate(comm, dm));
94: PetscCall(DMSetType(*dm, DMPLEX));
95: PetscCall(DMSetFromOptions(*dm));
96: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
101: {
102: options->degree = -1;
104: PetscFunctionBeginUser;
105: PetscOptionsBegin(comm, "", "Gauss' Law Test Options", "DMPLEX");
106: PetscOptionsEnd();
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user)
111: {
112: PetscFE feq, fep;
113: PetscSpace sp;
114: PetscQuadrature quad, fquad;
115: PetscDS ds;
116: DMLabel label;
117: DMPolytopeType ct;
118: PetscInt dim, cStart, minDeg, maxDeg;
119: PetscBool isTrimmed, isSum;
121: PetscFunctionBeginUser;
122: PetscCall(DMGetDimension(dm, &dim));
123: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
124: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
125: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, dim, ct, "field_", -1, &feq));
126: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
127: PetscCall(PetscFEGetQuadrature(feq, &quad));
128: PetscCall(PetscFEGetFaceQuadrature(feq, &fquad));
129: PetscCall(PetscFEGetBasisSpace(feq, &sp));
130: PetscCall(PetscSpaceGetDegree(sp, &minDeg, &maxDeg));
131: PetscCall(PetscObjectTypeCompare((PetscObject)sp, PETSCSPACEPTRIMMED, &isTrimmed));
132: PetscCall(PetscObjectTypeCompare((PetscObject)sp, PETSCSPACESUM, &isSum));
133: if (isSum) {
134: PetscSpace subsp, xsp, ysp;
135: PetscInt xdeg, ydeg;
136: PetscBool isTensor;
138: PetscCall(PetscSpaceSumGetSubspace(sp, 0, &subsp));
139: PetscCall(PetscObjectTypeCompare((PetscObject)subsp, PETSCSPACETENSOR, &isTensor));
140: if (isTensor) {
141: PetscCall(PetscSpaceTensorGetSubspace(subsp, 0, &xsp));
142: PetscCall(PetscSpaceTensorGetSubspace(subsp, 1, &ysp));
143: PetscCall(PetscSpaceGetDegree(xsp, &xdeg, NULL));
144: PetscCall(PetscSpaceGetDegree(ysp, &ydeg, NULL));
145: isTrimmed = xdeg != ydeg ? PETSC_TRUE : PETSC_FALSE;
146: }
147: }
148: user->degree = minDeg;
149: if (isTrimmed) user->divFree = PETSC_FALSE;
150: else user->divFree = PETSC_TRUE;
151: PetscCheck(!user->divFree || user->degree, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Degree 0 solution not available");
152: PetscCall(PetscFEDestroy(&feq));
153: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, "pot_", -1, &fep));
154: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fep));
155: PetscCall(PetscFESetQuadrature(fep, quad));
156: PetscCall(PetscFESetFaceQuadrature(fep, fquad));
157: PetscCall(PetscFEDestroy(&fep));
158: PetscCall(DMCreateDS(dm));
160: PetscCall(DMGetDS(dm, &ds));
161: PetscCall(PetscDSSetResidual(ds, 0, identity, NULL));
162: PetscCall(PetscDSSetResidual(ds, 1, divergence, NULL));
163: if (user->divFree) {
164: if (dim == 2) PetscCall(PetscDSSetExactSolution(ds, 0, solenoidal_totaldeg_2d, &user->degree));
165: else PetscCall(PetscDSSetExactSolution(ds, 0, solenoidal_totaldeg_3d, &user->degree));
166: } else {
167: PetscCall(PetscDSSetExactSolution(ds, 0, source_totaldeg, &user->degree));
168: }
169: PetscCall(PetscDSSetExactSolution(ds, 1, zero, &user->degree));
170: PetscCall(DMGetLabel(dm, "marker", &label));
172: // TODO Can we also test the boundary residual integration?
173: //PetscWeakForm wf;
174: //PetscInt bd, id = 1;
175: //PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "boundary", label, 1, &id, 1, 0, NULL, NULL, NULL, user, &bd));
176: //PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
177: //PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 1, 0, 0, flux, 0, NULL));
179: {
180: PetscScalar constants[1];
182: constants[0] = user->degree;
183: PetscCall(PetscDSSetConstants(ds, 1, constants));
184: }
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: int main(int argc, char **argv)
189: {
190: DM dm;
191: SNES snes;
192: Vec u;
193: PetscReal error[2], residual;
194: PetscScalar source[2], outflow[2];
195: AppCtx user;
197: PetscFunctionBeginUser;
198: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
199: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
200: PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
201: PetscCall(CreateDiscretization(dm, &user));
202: PetscCall(DMGetGlobalVector(dm, &u));
203: PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
204: PetscCall(DMComputeExactSolution(dm, 0., u, NULL));
206: PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
207: PetscCall(SNESSetDM(snes, dm));
208: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
209: PetscCall(SNESSetFromOptions(snes));
210: PetscCall(DMSNESCheckDiscretization(snes, dm, 0., u, PETSC_DETERMINE, error));
211: PetscCheck(PetscAbsReal(error[0]) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Exact solution does not fit into FEM space: %g should be zero", (double)error[0]);
212: if (user.divFree) {
213: PetscCall(DMSNESCheckResidual(snes, dm, u, PETSC_DETERMINE, &residual));
214: PetscCheck(PetscAbsReal(residual) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Exact solution is not divergence-free: %g should be zero", (double)residual);
215: } else {
216: PetscDS ds;
218: PetscCall(DMGetDS(dm, &ds));
219: PetscCall(PetscDSSetObjective(ds, 1, divergence));
220: PetscCall(DMPlexComputeIntegralFEM(dm, u, source, &user));
221: }
222: PetscCall(SNESDestroy(&snes));
224: PetscBdPointFunc funcs[] = {zero_bd, flux};
225: DMLabel label;
226: PetscInt id = 1;
228: PetscCall(DMGetLabel(dm, "marker", &label));
229: PetscCall(DMPlexComputeBdIntegral(dm, u, label, 1, &id, funcs, outflow, &user));
230: if (user.divFree) PetscCheck(PetscAbsScalar(outflow[1]) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Outflow %g should be zero for a divergence-free field", (double)PetscRealPart(outflow[1]));
231: else PetscCheck(PetscAbsScalar(source[1] + outflow[1]) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Outflow %g should oppose source %g", (double)PetscRealPart(outflow[1]), (double)PetscRealPart(source[1]));
233: PetscCall(DMRestoreGlobalVector(dm, &u));
234: PetscCall(DMDestroy(&dm));
235: PetscCall(PetscFinalize());
236: return 0;
237: }
239: /*TEST
241: testset:
242: suffix: p
243: requires: triangle ctetgen
244: args: -dm_plex_dim {{2 3}} -dm_plex_box_faces 2,2,2
246: test:
247: suffix: 1
248: args: -field_petscspace_degree 1 -pot_petscspace_degree 1
249: test:
250: suffix: 2
251: args: -field_petscspace_degree 2 -pot_petscspace_degree 2
252: test:
253: suffix: 3
254: args: -field_petscspace_degree 3 -pot_petscspace_degree 3
255: test:
256: suffix: 4
257: args: -field_petscspace_degree 4 -pot_petscspace_degree 4
259: testset:
260: suffix: q
261: args: -dm_plex_dim {{2 3}} -dm_plex_simplex 0 -dm_plex_box_faces 2,2
263: test:
264: suffix: 1
265: args: -field_petscspace_degree 1 -pot_petscspace_degree 1
266: test:
267: suffix: 2
268: args: -field_petscspace_degree 2 -pot_petscspace_degree 2
269: test:
270: suffix: 3
271: args: -field_petscspace_degree 3 -pot_petscspace_degree 3
272: test:
273: suffix: 4
274: args: -field_petscspace_degree 4 -pot_petscspace_degree 4
276: testset:
277: suffix: bdm
278: requires: triangle ctetgen
279: args: -dm_plex_dim 2 -dm_plex_box_faces 2,2
281: test:
282: suffix: 1
283: args: -pot_petscspace_degree 0 -pot_petscdualspace_lagrange_continuity 0 \
284: -field_petscspace_degree 1 -field_petscdualspace_type bdm \
285: -field_petscfe_default_quadrature_order 2
287: testset:
288: suffix: rt
289: requires: triangle ctetgen
290: args: -dm_plex_dim 2 -dm_plex_box_faces 2,2
292: test:
293: suffix: 1
294: args: -pot_petscspace_degree 0 -pot_petscdualspace_lagrange_continuity 0 \
295: -field_petscspace_type ptrimmed \
296: -field_petscspace_components 2 \
297: -field_petscspace_ptrimmed_form_degree -1 \
298: -field_petscdualspace_order 1 \
299: -field_petscdualspace_form_degree -1 \
300: -field_petscdualspace_lagrange_trimmed true \
301: -field_petscfe_default_quadrature_order 2
303: testset:
304: suffix: rtq
305: requires: triangle ctetgen
306: args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
308: test:
309: suffix: 1
310: args: -pot_petscspace_degree 0 -pot_petscdualspace_lagrange_continuity 0 \
311: -field_petscspace_degree 1 \
312: -field_petscspace_type sum \
313: -field_petscspace_variables 2 \
314: -field_petscspace_components 2 \
315: -field_petscspace_sum_spaces 2 \
316: -field_petscspace_sum_concatenate true \
317: -field_sumcomp_0_petscspace_variables 2 \
318: -field_sumcomp_0_petscspace_type tensor \
319: -field_sumcomp_0_petscspace_tensor_spaces 2 \
320: -field_sumcomp_0_petscspace_tensor_uniform false \
321: -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
322: -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
323: -field_sumcomp_1_petscspace_variables 2 \
324: -field_sumcomp_1_petscspace_type tensor \
325: -field_sumcomp_1_petscspace_tensor_spaces 2 \
326: -field_sumcomp_1_petscspace_tensor_uniform false \
327: -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
328: -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
329: -field_petscdualspace_order 1 \
330: -field_petscdualspace_form_degree -1 \
331: -field_petscdualspace_lagrange_trimmed true \
332: -field_petscfe_default_quadrature_order 2
334: TEST*/