1 static char help[] = "'Good Cop' Helmholtz Problem in 2d and 3d with finite elements.\n\ 2 We solve the 'Good Cop' Helmholtz problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports automatic convergence estimation\n\ 5 and coarse space adaptivity.\n\n\n"; 6 7 /* 8 The model problem: 9 Solve "Good Cop" Helmholtz equation on the unit square: (0,1) x (0,1) 10 - \Delta u + u = f, 11 where \Delta = Laplace operator 12 Dirichlet b.c.'s on all sides 13 14 */ 15 16 #include <petscdmplex.h> 17 #include <petscsnes.h> 18 #include <petscds.h> 19 #include <petscconvest.h> 20 21 typedef struct { 22 PetscBool trig; /* Use trig function as exact solution */ 23 } AppCtx; 24 25 /* For Primal Problem */ 26 static void g0_uu(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[]) { 27 PetscInt d; 28 for (d = 0; d < dim; ++d) g0[0] = 1.0; 29 } 30 31 static void g3_uu(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 g3[]) { 32 PetscInt d; 33 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 34 } 35 36 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 37 PetscInt d; 38 *u = 0.0; 39 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]); 40 return 0; 41 } 42 43 static PetscErrorCode quad_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 44 PetscInt d; 45 *u = 1.0; 46 for (d = 0; d < dim; ++d) *u += (d + 1) * PetscSqr(x[d]); 47 return 0; 48 } 49 50 static void f0_trig_u(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[]) { 51 PetscInt d; 52 f0[0] += u[0]; 53 for (d = 0; d < dim; ++d) f0[0] -= 4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]) + PetscSinReal(2.0 * PETSC_PI * x[d]); 54 } 55 56 static void f0_quad_u(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[]) { 57 PetscInt d; 58 switch (dim) { 59 case 1: f0[0] = 1.0; break; 60 case 2: f0[0] = 5.0; break; 61 case 3: f0[0] = 11.0; break; 62 default: f0[0] = 5.0; break; 63 } 64 f0[0] += u[0]; 65 for (d = 0; d < dim; ++d) f0[0] -= (d + 1) * PetscSqr(x[d]); 66 } 67 68 static void f1_u(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 f1[]) { 69 PetscInt d; 70 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 71 } 72 73 static PetscErrorCode ProcessOptions(DM dm, AppCtx *options) { 74 MPI_Comm comm; 75 PetscInt dim; 76 77 PetscFunctionBeginUser; 78 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 79 PetscCall(DMGetDimension(dm, &dim)); 80 options->trig = PETSC_FALSE; 81 82 PetscOptionsBegin(comm, "", "Helmholtz Problem Options", "DMPLEX"); 83 PetscCall(PetscOptionsBool("-exact_trig", "Use trigonometric exact solution (better for more complex finite elements)", "ex26.c", options->trig, &options->trig, NULL)); 84 PetscOptionsEnd(); 85 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) { 90 PetscFunctionBeginUser; 91 PetscCall(DMCreate(comm, dm)); 92 PetscCall(DMSetType(*dm, DMPLEX)); 93 PetscCall(DMSetFromOptions(*dm)); 94 95 PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh")); 96 PetscCall(DMSetApplicationContext(*dm, user)); 97 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 98 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) { 103 PetscDS ds; 104 DMLabel label; 105 const PetscInt id = 1; 106 107 PetscFunctionBeginUser; 108 PetscCall(DMGetDS(dm, &ds)); 109 PetscCall(DMGetLabel(dm, "marker", &label)); 110 if (user->trig) { 111 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u)); 112 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu)); 113 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user)); 114 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL)); 115 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n")); 116 } else { 117 PetscCall(PetscDSSetResidual(ds, 0, f0_quad_u, f1_u)); 118 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu)); 119 PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user)); 120 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL)); 121 } 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) { 126 DM cdm = dm; 127 PetscFE fe; 128 DMPolytopeType ct; 129 PetscBool simplex; 130 PetscInt dim, cStart; 131 char prefix[PETSC_MAX_PATH_LEN]; 132 133 PetscFunctionBeginUser; 134 PetscCall(DMGetDimension(dm, &dim)); 135 136 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 137 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 138 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE; 139 /* Create finite element */ 140 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 141 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe)); 142 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 143 /* Set discretization and boundary conditions for each mesh */ 144 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 145 PetscCall(DMCreateDS(dm)); 146 PetscCall((*setup)(dm, user)); 147 while (cdm) { 148 PetscCall(DMCopyDisc(dm, cdm)); 149 PetscCall(DMGetCoarseDM(cdm, &cdm)); 150 } 151 PetscCall(PetscFEDestroy(&fe)); 152 PetscFunctionReturn(0); 153 } 154 155 int main(int argc, char **argv) { 156 DM dm; /* Problem specification */ 157 PetscDS ds; 158 SNES snes; /* Nonlinear solver */ 159 Vec u; /* Solutions */ 160 AppCtx user; /* User-defined work context */ 161 162 PetscFunctionBeginUser; 163 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 164 /* Primal system */ 165 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 166 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 167 PetscCall(ProcessOptions(dm, &user)); 168 PetscCall(SNESSetDM(snes, dm)); 169 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 170 PetscCall(DMCreateGlobalVector(dm, &u)); 171 PetscCall(VecSet(u, 0.0)); 172 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 173 PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 174 PetscCall(SNESSetFromOptions(snes)); 175 PetscCall(DMSNESCheckFromOptions(snes, u)); 176 177 /* Looking for field error */ 178 PetscInt Nfields; 179 PetscCall(DMGetDS(dm, &ds)); 180 PetscCall(PetscDSGetNumFields(ds, &Nfields)); 181 PetscCall(SNESSolve(snes, NULL, u)); 182 PetscCall(SNESGetSolution(snes, &u)); 183 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 184 185 /* Cleanup */ 186 PetscCall(VecDestroy(&u)); 187 PetscCall(SNESDestroy(&snes)); 188 PetscCall(DMDestroy(&dm)); 189 PetscCall(PetscFinalize()); 190 return 0; 191 } 192 193 /*TEST 194 test: 195 # L_2 convergence rate: 1.9 196 suffix: 2d_p1_conv 197 requires: triangle 198 args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu 199 test: 200 # L_2 convergence rate: 1.9 201 suffix: 2d_q1_conv 202 args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu 203 test: 204 # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5 205 suffix: 3d_p1_conv 206 requires: ctetgen 207 args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu 208 test: 209 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2 210 suffix: 3d_q1_conv 211 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu 212 test: 213 # L_2 convergence rate: 1.9 214 suffix: 2d_p1_trig_conv 215 requires: triangle 216 args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig 217 test: 218 # L_2 convergence rate: 1.9 219 suffix: 2d_q1_trig_conv 220 args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig 221 test: 222 # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5 223 suffix: 3d_p1_trig_conv 224 requires: ctetgen 225 args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig 226 test: 227 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2 228 suffix: 3d_q1_trig_conv 229 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig 230 test: 231 suffix: 2d_p1_gmg_vcycle 232 requires: triangle 233 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 234 -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \ 235 -mg_levels_ksp_max_it 1 \ 236 -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor 237 test: 238 suffix: 2d_p1_gmg_fcycle 239 requires: triangle 240 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 241 -ksp_type cg -ksp_rtol 1e-10 -pc_type mg -pc_mg_type full \ 242 -mg_levels_ksp_max_it 2 \ 243 -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor 244 test: 245 suffix: 2d_p1_gmg_vcycle_trig 246 requires: triangle 247 args: -exact_trig -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 248 -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \ 249 -mg_levels_ksp_max_it 1 \ 250 -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor 251 test: 252 suffix: 2d_q3_cgns 253 requires: cgns 254 args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_box_faces 3,3 -snes_view_solution cgns:sol.cgns -potential_petscspace_degree 3 -dm_coord_petscspace_degree 3 255 TEST*/ 256