1 static char help[] = "Tests for Gauss' Law\n\n";
2
3 /* We want to check the weak version of Gauss' Law, namely that
4
5 \int_\Omega v div q - \int_\Gamma v (q \cdot n) = 0
6
7 */
8
9 #include <petscdmplex.h>
10 #include <petscsnes.h>
11 #include <petscds.h>
12
13 typedef struct {
14 PetscInt degree; // The degree of the discretization
15 PetscBool divFree; // True if the solution is divergence-free
16 } AppCtx;
17
zero(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)18 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
19 {
20 u[0] = 0.0;
21 return PETSC_SUCCESS;
22 }
23
24 // div <x^d y^{d-1}, -x^{d-1} y^d> = 0
solenoidal_2d(PetscInt n,const PetscReal x[],PetscScalar u[])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 <x^d y^{d-1} z^{d-1}, -2 x^{d-1} y^d z^{d-1}, x^{d-1} y^{d-1} z^d> = 0
solenoidal_3d(PetscInt n,const PetscReal x[],PetscScalar u[])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 }
37
solenoidal_totaldeg_2d(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)38 static PetscErrorCode solenoidal_totaldeg_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
39 {
40 const PetscInt deg = *(PetscInt *)ctx;
41 const PetscInt n = deg / 2 + deg % 2;
42
43 solenoidal_2d(n, x, u);
44 return PETSC_SUCCESS;
45 }
46
solenoidal_totaldeg_3d(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)47 static PetscErrorCode solenoidal_totaldeg_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
48 {
49 const PetscInt deg = *(PetscInt *)ctx;
50 const PetscInt n = deg / 3 + (deg % 3 ? 1 : 0);
51
52 solenoidal_3d(n, x, u);
53 return PETSC_SUCCESS;
54 }
55
56 // This is in P_n^{-}
source_totaldeg(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)57 static PetscErrorCode source_totaldeg(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
58 {
59 const PetscInt n = *(PetscInt *)ctx;
60
61 for (PetscInt d = 0; d < dim; ++d) u[d] = PetscPowRealInt(x[d], n + 1);
62 return PETSC_SUCCESS;
63 }
64
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[])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];
69
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 }
74
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[])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 }
79
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[])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 }
84
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[])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 }
89
CreateMesh(MPI_Comm comm,DM * dm)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 }
99
ProcessOptions(MPI_Comm comm,AppCtx * options)100 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
101 {
102 options->degree = -1;
103
104 PetscFunctionBeginUser;
105 PetscOptionsBegin(comm, "", "Gauss' Law Test Options", "DMPLEX");
106 PetscOptionsEnd();
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
CreateDiscretization(DM dm,AppCtx * user)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;
120
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;
137
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));
159
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));
171
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));
178
179 {
180 PetscScalar constants[1];
181
182 constants[0] = user->degree;
183 PetscCall(PetscDSSetConstants(ds, 1, constants));
184 }
185 PetscFunctionReturn(PETSC_SUCCESS);
186 }
187
main(int argc,char ** argv)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;
196
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));
205
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;
217
218 PetscCall(DMGetDS(dm, &ds));
219 PetscCall(PetscDSSetObjective(ds, 1, divergence));
220 PetscCall(DMPlexComputeIntegralFEM(dm, u, source, &user));
221 }
222 PetscCall(SNESDestroy(&snes));
223
224 PetscBdPointFn *funcs[] = {zero_bd, flux};
225 DMLabel label;
226 PetscInt id = 1;
227
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]));
232
233 PetscCall(DMRestoreGlobalVector(dm, &u));
234 PetscCall(DMDestroy(&dm));
235 PetscCall(PetscFinalize());
236 return 0;
237 }
238
239 /*TEST
240
241 testset:
242 suffix: p
243 requires: triangle ctetgen
244 args: -dm_plex_dim {{2 3}} -dm_plex_box_faces 2,2,2
245 output_file: output/empty.out
246
247 test:
248 suffix: 1
249 args: -field_petscspace_degree 1 -pot_petscspace_degree 1
250 test:
251 suffix: 2
252 args: -field_petscspace_degree 2 -pot_petscspace_degree 2
253 test:
254 suffix: 3
255 args: -field_petscspace_degree 3 -pot_petscspace_degree 3
256 test:
257 suffix: 4
258 args: -field_petscspace_degree 4 -pot_petscspace_degree 4
259
260 testset:
261 suffix: q
262 args: -dm_plex_dim {{2 3}} -dm_plex_simplex 0 -dm_plex_box_faces 2,2
263 output_file: output/empty.out
264
265 test:
266 suffix: 1
267 args: -field_petscspace_degree 1 -pot_petscspace_degree 1
268 test:
269 suffix: 2
270 args: -field_petscspace_degree 2 -pot_petscspace_degree 2
271 test:
272 suffix: 3
273 args: -field_petscspace_degree 3 -pot_petscspace_degree 3
274 test:
275 suffix: 4
276 args: -field_petscspace_degree 4 -pot_petscspace_degree 4
277
278 testset:
279 suffix: bdm
280 requires: triangle ctetgen
281 args: -dm_plex_dim 2 -dm_plex_box_faces 2,2
282 output_file: output/empty.out
283
284 test:
285 suffix: 1
286 args: -pot_petscspace_degree 0 -pot_petscdualspace_lagrange_continuity 0 \
287 -field_petscspace_degree 1 -field_petscdualspace_type bdm \
288 -field_petscfe_default_quadrature_order 2
289
290 testset:
291 suffix: rt
292 requires: triangle ctetgen
293 args: -dm_plex_dim 2 -dm_plex_box_faces 2,2
294 output_file: output/empty.out
295
296 test:
297 suffix: 1
298 args: -pot_petscspace_degree 0 -pot_petscdualspace_lagrange_continuity 0 \
299 -field_petscspace_type ptrimmed \
300 -field_petscspace_components 2 \
301 -field_petscspace_ptrimmed_form_degree -1 \
302 -field_petscdualspace_order 1 \
303 -field_petscdualspace_form_degree -1 \
304 -field_petscdualspace_lagrange_trimmed true \
305 -field_petscfe_default_quadrature_order 2
306
307 testset:
308 suffix: rtq
309 requires: triangle ctetgen
310 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
311 output_file: output/empty.out
312
313 test:
314 suffix: 1
315 args: -pot_petscspace_degree 0 -pot_petscdualspace_lagrange_continuity 0 \
316 -field_petscspace_degree 1 \
317 -field_petscspace_type sum \
318 -field_petscspace_variables 2 \
319 -field_petscspace_components 2 \
320 -field_petscspace_sum_spaces 2 \
321 -field_petscspace_sum_concatenate true \
322 -field_sumcomp_0_petscspace_variables 2 \
323 -field_sumcomp_0_petscspace_type tensor \
324 -field_sumcomp_0_petscspace_tensor_spaces 2 \
325 -field_sumcomp_0_petscspace_tensor_uniform false \
326 -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
327 -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
328 -field_sumcomp_1_petscspace_variables 2 \
329 -field_sumcomp_1_petscspace_type tensor \
330 -field_sumcomp_1_petscspace_tensor_spaces 2 \
331 -field_sumcomp_1_petscspace_tensor_uniform false \
332 -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
333 -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
334 -field_petscdualspace_order 1 \
335 -field_petscdualspace_form_degree -1 \
336 -field_petscdualspace_lagrange_trimmed true \
337 -field_petscfe_default_quadrature_order 2
338
339 TEST*/
340