1 static const char help[] = "Tests for determining whether a new finite element works"; 2 3 /* 4 Use -interpolation_view and -l2_projection_view to look at the interpolants. 5 */ 6 7 #include <petscdmplex.h> 8 #include <petscfe.h> 9 #include <petscds.h> 10 #include <petscsnes.h> 11 12 static void constant(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[]) { 13 const PetscInt Nc = uOff[1] - uOff[0]; 14 for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5.; 15 } 16 17 static void linear(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[]) { 18 const PetscInt Nc = uOff[1] - uOff[0]; 19 for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5. * x[c]; 20 } 21 22 static void quadratic(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[]) { 23 const PetscInt Nc = uOff[1] - uOff[0]; 24 for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5. * x[c] * x[c]; 25 } 26 27 static void trig(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[]) { 28 const PetscInt Nc = uOff[1] - uOff[0]; 29 for (PetscInt c = 0; c < Nc; ++c) f0[c] += PetscCosReal(2. * PETSC_PI * x[c]); 30 } 31 32 /* 33 The prime basis for the Wheeler-Yotov-Xue prism. 34 */ 35 static void prime(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[]) { 36 PetscReal x = X[0], y = X[1], z = X[2], b = 1 + x + y + z; 37 f0[0] += b + 2.0 * x * z + 2.0 * y * z + x * y + x * x; 38 f0[1] += b + 2.0 * x * z + 2.0 * y * z + x * y + y * y; 39 f0[2] += b - 3.0 * x * z - 3.0 * y * z - 2.0 * z * z; 40 } 41 42 static const char *names[] = {"constant", "linear", "quadratic", "trig", "prime"}; 43 static PetscPointFunc functions[] = {constant, linear, quadratic, trig, prime}; 44 45 typedef struct { 46 PetscPointFunc exactSol; 47 PetscReal shear, flatten; 48 } AppCtx; 49 50 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 51 char name[PETSC_MAX_PATH_LEN] = "constant"; 52 PetscInt Nfunc = PETSC_STATIC_ARRAY_LENGTH(names), i; 53 54 PetscFunctionBeginUser; 55 options->exactSol = NULL; 56 options->shear = 0.; 57 options->flatten = 1.; 58 59 PetscOptionsBegin(comm, "", "FE Test Options", "PETSCFE"); 60 PetscCall(PetscOptionsString("-func", "Function to project into space", "", name, name, PETSC_MAX_PATH_LEN, NULL)); 61 PetscCall(PetscOptionsReal("-shear", "Factor by which to shear along the x-direction", "", options->shear, &(options->shear), NULL)); 62 PetscCall(PetscOptionsReal("-flatten", "Factor by which to flatten", "", options->flatten, &(options->flatten), NULL)); 63 PetscOptionsEnd(); 64 65 for (i = 0; i < Nfunc; ++i) { 66 PetscBool flg; 67 68 PetscCall(PetscStrcmp(name, names[i], &flg)); 69 if (flg) { 70 options->exactSol = functions[i]; 71 break; 72 } 73 } 74 PetscCheck(options->exactSol, comm, PETSC_ERR_ARG_WRONG, "Invalid test function %s", name); 75 PetscFunctionReturn(0); 76 } 77 78 /* The exact solution is the negative of the f0 contribution */ 79 static PetscErrorCode exactSolution(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 80 AppCtx *user = (AppCtx *)ctx; 81 PetscInt uOff[2] = {0, Nc}; 82 83 user->exactSol(dim, 1, 0, uOff, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, time, x, 0, NULL, u); 84 for (PetscInt c = 0; c < Nc; ++c) u[c] *= -1.; 85 return 0; 86 } 87 88 static void f0(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[]) { 89 const PetscInt Nc = uOff[1] - uOff[0]; 90 for (PetscInt c = 0; c < Nc; ++c) f0[c] += u[c]; 91 } 92 93 static void g0(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[]) { 94 const PetscInt Nc = uOff[1] - uOff[0]; 95 for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0; 96 } 97 98 static PetscErrorCode SetupProblem(DM dm, AppCtx *user) { 99 PetscDS ds; 100 PetscWeakForm wf; 101 102 PetscFunctionBeginUser; 103 PetscCall(DMGetDS(dm, &ds)); 104 PetscCall(PetscDSGetWeakForm(ds, &wf)); 105 PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0, 0, NULL)); 106 PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, user->exactSol, 0, NULL)); 107 PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0, 0, NULL, 0, NULL, 0, NULL)); 108 PetscCall(PetscDSSetExactSolution(ds, 0, exactSolution, user)); 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode SetupDiscretization(DM dm, const char name[], AppCtx *user) { 113 DM cdm = dm; 114 PetscFE fe; 115 char prefix[PETSC_MAX_PATH_LEN]; 116 117 PetscFunctionBeginUser; 118 if (name) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%s_", name)); 119 PetscCall(DMCreateFEDefault(dm, 1, name ? prefix : NULL, -1, &fe)); 120 PetscCall(PetscObjectSetName((PetscObject)fe, name ? name : "Solution")); 121 /* Set discretization and boundary conditions for each mesh */ 122 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 123 PetscCall(DMCreateDS(dm)); 124 PetscCall(SetupProblem(dm, user)); 125 while (cdm) { 126 PetscCall(DMCopyDisc(dm, cdm)); 127 PetscCall(DMGetCoarseDM(cdm, &cdm)); 128 } 129 PetscCall(PetscFEDestroy(&fe)); 130 PetscFunctionReturn(0); 131 } 132 133 /* This test tells us whether the given function is contained in the approximation space */ 134 static PetscErrorCode CheckInterpolation(DM dm, AppCtx *user) { 135 PetscSimplePointFunc exactSol[1]; 136 void *exactCtx[1]; 137 PetscDS ds; 138 Vec u; 139 PetscReal error, tol = PETSC_SMALL; 140 MPI_Comm comm; 141 142 PetscFunctionBeginUser; 143 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 144 PetscCall(DMGetDS(dm, &ds)); 145 PetscCall(DMGetGlobalVector(dm, &u)); 146 PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0])); 147 PetscCall(DMProjectFunction(dm, 0.0, exactSol, exactCtx, INSERT_ALL_VALUES, u)); 148 PetscCall(VecViewFromOptions(u, NULL, "-interpolation_view")); 149 PetscCall(DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error)); 150 PetscCall(DMRestoreGlobalVector(dm, &u)); 151 if (error > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL at tolerance %g error %g\n", (double)tol, (double)error)); 152 else PetscCall(PetscPrintf(comm, "Interpolation tests pass at tolerance %g\n", (double)tol)); 153 PetscFunctionReturn(0); 154 } 155 156 /* This test tells us whether the element is unisolvent (the mass matrix has full rank), and what rate of convergence we achieve */ 157 static PetscErrorCode CheckL2Projection(DM dm, AppCtx *user) { 158 PetscSimplePointFunc exactSol[1]; 159 void *exactCtx[1]; 160 SNES snes; 161 PetscDS ds; 162 Vec u; 163 PetscReal error, tol = PETSC_SMALL; 164 MPI_Comm comm; 165 166 PetscFunctionBeginUser; 167 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 168 PetscCall(DMGetDS(dm, &ds)); 169 PetscCall(DMGetGlobalVector(dm, &u)); 170 PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0])); 171 PetscCall(SNESCreate(comm, &snes)); 172 PetscCall(SNESSetDM(snes, dm)); 173 PetscCall(VecSet(u, 0.0)); 174 PetscCall(PetscObjectSetName((PetscObject)u, "solution")); 175 PetscCall(DMPlexSetSNESLocalFEM(dm, user, user, user)); 176 PetscCall(SNESSetFromOptions(snes)); 177 PetscCall(DMSNESCheckFromOptions(snes, u)); 178 PetscCall(SNESSolve(snes, NULL, u)); 179 PetscCall(SNESDestroy(&snes)); 180 PetscCall(VecViewFromOptions(u, NULL, "-l2_projection_view")); 181 PetscCall(DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error)); 182 PetscCall(DMRestoreGlobalVector(dm, &u)); 183 if (error > tol) PetscCall(PetscPrintf(comm, "L2 projection tests FAIL at tolerance %g error %g\n", (double)tol, (double)error)); 184 else PetscCall(PetscPrintf(comm, "L2 projection tests pass at tolerance %g\n", (double)tol)); 185 PetscFunctionReturn(0); 186 } 187 188 /* Distorts the mesh by shearing in the x-direction and flattening, factors provided in the options. */ 189 static PetscErrorCode DistortMesh(DM dm, AppCtx *user) { 190 Vec coordinates; 191 PetscScalar *ca; 192 PetscInt dE, n, i; 193 194 PetscFunctionBeginUser; 195 PetscCall(DMGetCoordinateDim(dm, &dE)); 196 PetscCall(DMGetCoordinates(dm, &coordinates)); 197 PetscCall(VecGetLocalSize(coordinates, &n)); 198 PetscCall(VecGetArray(coordinates, &ca)); 199 for (i = 0; i < (n / dE); ++i) { 200 ca[i * dE + 0] += user->shear * ca[i * dE + 0]; 201 ca[i * dE + 1] *= user->flatten; 202 } 203 PetscCall(VecRestoreArray(coordinates, &ca)); 204 PetscFunctionReturn(0); 205 } 206 207 int main(int argc, char **argv) { 208 DM dm; 209 AppCtx user; 210 PetscMPIInt size; 211 212 PetscFunctionBeginUser; 213 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 214 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 215 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only."); 216 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 217 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 218 PetscCall(DMSetType(dm, DMPLEX)); 219 PetscCall(DMSetFromOptions(dm)); 220 PetscCall(DistortMesh(dm, &user)); 221 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 222 PetscCall(SetupDiscretization(dm, NULL, &user)); 223 224 PetscCall(CheckInterpolation(dm, &user)); 225 PetscCall(CheckL2Projection(dm, &user)); 226 227 PetscCall(DMDestroy(&dm)); 228 PetscCall(PetscFinalize()); 229 return 0; 230 } 231 232 /*TEST 233 234 testset: 235 args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -petscspace_degree 1\ 236 -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu 237 238 test: 239 suffix: p1_0 240 args: -func {{constant linear}} 241 242 # Using -dm_refine 2 -convest_num_refine 4 gives convergence rate 2.0 243 test: 244 suffix: p1_1 245 args: -func {{quadratic trig}} \ 246 -snes_convergence_estimate -convest_num_refine 2 247 248 testset: 249 requires: !complex double 250 args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \ 251 -petscspace_type sum \ 252 -petscspace_variables 3 \ 253 -petscspace_components 3 \ 254 -petscspace_sum_spaces 2 \ 255 -petscspace_sum_concatenate false \ 256 -sumcomp_0_petscspace_variables 3 \ 257 -sumcomp_0_petscspace_components 3 \ 258 -sumcomp_0_petscspace_degree 1 \ 259 -sumcomp_1_petscspace_variables 3 \ 260 -sumcomp_1_petscspace_components 3 \ 261 -sumcomp_1_petscspace_type wxy \ 262 -petscdualspace_form_degree 0 \ 263 -petscdualspace_order 1 \ 264 -petscdualspace_components 3 \ 265 -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu 266 267 test: 268 suffix: wxy_0 269 args: -func constant 270 271 test: 272 suffix: wxy_1 273 args: -func linear 274 275 test: 276 suffix: wxy_2 277 args: -func prime 278 279 test: 280 suffix: wxy_3 281 args: -func linear -shear 1 -flatten 1e-5 282 283 TEST*/ 284