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