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 static const char *names[] = {"constant", "linear", "quadratic", "trig"}; 49 static PetscPointFunc functions[] = { constant, linear, quadratic, trig}; 50 51 typedef struct { 52 PetscPointFunc exactSol; 53 } AppCtx; 54 55 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 56 { 57 char name[PETSC_MAX_PATH_LEN] = "constant"; 58 PetscInt Nfunc = sizeof(names)/sizeof(char *), i; 59 PetscErrorCode ierr; 60 61 PetscFunctionBeginUser; 62 options->exactSol = NULL; 63 64 ierr = PetscOptionsBegin(comm, "", "FE Test Options", "PETSCFE");CHKERRQ(ierr); 65 ierr = PetscOptionsString("-func", "Function to project into space", "", name, name, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 66 ierr = PetscOptionsEnd();CHKERRQ(ierr); 67 68 for (i = 0; i < Nfunc; ++i) { 69 PetscBool flg; 70 71 ierr = PetscStrcmp(name, names[i], &flg);CHKERRQ(ierr); 72 if (flg) {options->exactSol = functions[i]; break;} 73 } 74 if (!options->exactSol) SETERRQ1(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 { 81 AppCtx *user = (AppCtx *) ctx; 82 PetscInt uOff[2] = {0., Nc}; 83 84 user->exactSol(dim, 1, 0, uOff, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, time, x, 0, NULL, u); 85 for (PetscInt c = 0; c < Nc; ++c) u[c] *= -1.; 86 return 0; 87 } 88 89 static void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux, 90 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 91 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 92 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 93 { 94 const PetscInt Nc = uOff[1] - uOff[0]; 95 for (PetscInt c = 0; c < Nc; ++c) f0[c] += u[c]; 96 } 97 98 static void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux, 99 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 100 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 101 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 PetscErrorCode ierr; 112 113 PetscFunctionBeginUser; 114 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 115 ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr); 116 ierr = PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0, 0, NULL);CHKERRQ(ierr); 117 ierr = PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, user->exactSol, 0, NULL);CHKERRQ(ierr); 118 ierr = PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0, 0, NULL, 0, NULL, 0, NULL);CHKERRQ(ierr); 119 ierr = PetscDSSetExactSolution(ds, 0, exactSolution, user);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 static PetscErrorCode SetupDiscretization(DM dm, const char name[], AppCtx *user) 124 { 125 DM cdm = dm; 126 PetscFE fe; 127 char prefix[PETSC_MAX_PATH_LEN]; 128 PetscErrorCode ierr; 129 130 PetscFunctionBeginUser; 131 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 132 ierr = DMCreateFEDefault(dm, 1, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 133 ierr = PetscObjectSetName((PetscObject) fe, name ? name : "Solution");CHKERRQ(ierr); 134 /* Set discretization and boundary conditions for each mesh */ 135 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 136 ierr = DMCreateDS(dm);CHKERRQ(ierr); 137 ierr = SetupProblem(dm, user);CHKERRQ(ierr); 138 while (cdm) { 139 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 140 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 141 } 142 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 146 /* This test tells us whether the given function is contained in the approximation space */ 147 static PetscErrorCode CheckInterpolation(DM dm, AppCtx *user) 148 { 149 PetscSimplePointFunc exactSol[1]; 150 void *exactCtx[1]; 151 PetscDS ds; 152 Vec u; 153 PetscReal error, tol = PETSC_SMALL; 154 MPI_Comm comm; 155 PetscErrorCode ierr; 156 157 PetscFunctionBeginUser; 158 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 159 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 160 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 161 ierr = PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]);CHKERRQ(ierr); 162 ierr = DMProjectFunction(dm, 0.0, exactSol, exactCtx, INSERT_ALL_VALUES, u);CHKERRQ(ierr); 163 ierr = VecViewFromOptions(u, NULL, "-interpolation_view");CHKERRQ(ierr); 164 ierr = DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);CHKERRQ(ierr); 165 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 166 if (error > tol) {ierr = PetscPrintf(comm, "Interpolation tests FAIL at tolerance %g error %g\n", (double) tol, (double) error);CHKERRQ(ierr);} 167 else {ierr = PetscPrintf(comm, "Interpolation tests pass at tolerance %g\n", (double) tol);CHKERRQ(ierr);} 168 PetscFunctionReturn(0); 169 } 170 171 /* This test tells us whether the element is unisolvent (the mass matrix has full rank), and what rate of convergence we achieve */ 172 static PetscErrorCode CheckL2Projection(DM dm, AppCtx *user) 173 { 174 PetscSimplePointFunc exactSol[1]; 175 void *exactCtx[1]; 176 SNES snes; 177 PetscDS ds; 178 Vec u; 179 PetscReal error, tol = PETSC_SMALL; 180 MPI_Comm comm; 181 PetscErrorCode ierr; 182 183 PetscFunctionBeginUser; 184 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 185 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 186 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 187 ierr = PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]);CHKERRQ(ierr); 188 ierr = SNESCreate(comm, &snes);CHKERRQ(ierr); 189 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 190 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 191 ierr = PetscObjectSetName((PetscObject) u, "solution");CHKERRQ(ierr); 192 ierr = DMPlexSetSNESLocalFEM(dm, user, user, user);CHKERRQ(ierr); 193 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 194 ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr); 195 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 196 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 197 ierr = VecViewFromOptions(u, NULL, "-l2_projection_view");CHKERRQ(ierr); 198 ierr = DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);CHKERRQ(ierr); 199 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 200 if (error > tol) {ierr = PetscPrintf(comm, "L2 projection tests FAIL at tolerance %g error %g\n", (double) tol, (double) error);CHKERRQ(ierr);} 201 else {ierr = PetscPrintf(comm, "L2 projection tests pass at tolerance %g\n", (double) tol);CHKERRQ(ierr);} 202 PetscFunctionReturn(0); 203 } 204 205 int main(int argc, char **argv) 206 { 207 DM dm; 208 AppCtx user; 209 PetscMPIInt size; 210 PetscErrorCode ierr; 211 212 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 213 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 214 if (size > 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only."); 215 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 216 ierr = DMCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr); 217 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 218 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 219 ierr = SetupDiscretization(dm, NULL, &user);CHKERRQ(ierr); 220 221 ierr = CheckInterpolation(dm, &user);CHKERRQ(ierr); 222 ierr = CheckL2Projection(dm, &user);CHKERRQ(ierr); 223 224 ierr = DMDestroy(&dm);CHKERRQ(ierr); 225 ierr = PetscFinalize(); 226 return ierr; 227 } 228 229 /*TEST 230 231 testset: 232 args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -petscspace_degree 1\ 233 -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu 234 235 test: 236 suffix: p1_0 237 args: -func {{constant linear}} 238 239 # Using -dm_refine 2 -convest_num_refine 4 gives convergence rate 2.0 240 test: 241 suffix: p1_1 242 args: -func {{quadratic trig}} \ 243 -snes_convergence_estimate -convest_num_refine 2 244 245 testset: 246 args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \ 247 -petscspace_type sum \ 248 -petscspace_variables 3 \ 249 -petscspace_components 3 \ 250 -petscspace_sum_spaces 2 \ 251 -petscspace_sum_concatenate false \ 252 -sumcomp_0_petscspace_variables 3 \ 253 -sumcomp_0_petscspace_components 3 \ 254 -sumcomp_0_petscspace_degree 1 \ 255 -sumcomp_1_petscspace_variables 3 \ 256 -sumcomp_1_petscspace_components 3 \ 257 -sumcomp_1_petscspace_type wxy \ 258 -petscdualspace_form_degree 0 \ 259 -petscdualspace_order 1 \ 260 -petscdualspace_components 3 \ 261 -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu 262 263 test: 264 suffix: wxy_0 265 args: -func constant 266 267 test: 268 suffix: wxy_1 269 args: -func linear 270 271 TEST*/ 272