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