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