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 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 } 23 24 // div <x^d y^{d-1}, -x^{d-1} y^d> = 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 <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 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 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; 42 43 solenoidal_2d(n, x, u); 44 return PETSC_SUCCESS; 45 } 46 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); 51 52 solenoidal_3d(n, x, u); 53 return PETSC_SUCCESS; 54 } 55 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; 60 61 for (PetscInt d = 0; d < dim; ++d) u[d] = PetscPowRealInt(x[d], n + 1); 62 return PETSC_SUCCESS; 63 } 64 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 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 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 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 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 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 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 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 PetscBdPointFunc 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 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 258 259 testset: 260 suffix: q 261 args: -dm_plex_dim {{2 3}} -dm_plex_simplex 0 -dm_plex_box_faces 2,2 262 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 275 276 testset: 277 suffix: bdm 278 requires: triangle ctetgen 279 args: -dm_plex_dim 2 -dm_plex_box_faces 2,2 280 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 286 287 testset: 288 suffix: rt 289 requires: triangle ctetgen 290 args: -dm_plex_dim 2 -dm_plex_box_faces 2,2 291 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 302 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 307 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 333 334 TEST*/ 335