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