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 *)&deg));
 71:   else PetscCallVoid(solenoidal_totaldeg_3d(dim, t, x, uOff[1] - uOff[0], p, (void *)&deg));
 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*/