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 0; 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 0; 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 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 106 { 107 PetscFunctionBeginUser; 108 PetscCall(DMCreate(comm, dm)); 109 PetscCall(DMSetType(*dm, DMPLEX)); 110 PetscCall(DMSetFromOptions(*dm)); 111 112 PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh")); 113 PetscCall(DMSetApplicationContext(*dm, user)); 114 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 115 116 PetscFunctionReturn(0); 117 } 118 119 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 120 { 121 PetscDS ds; 122 DMLabel label; 123 const PetscInt id = 1; 124 125 PetscFunctionBeginUser; 126 PetscCall(DMGetDS(dm, &ds)); 127 PetscCall(DMGetLabel(dm, "marker", &label)); 128 if (user->trig) { 129 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u)); 130 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu)); 131 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user)); 132 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL)); 133 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n")); 134 } else { 135 PetscCall(PetscDSSetResidual(ds, 0, f0_quad_u, f1_u)); 136 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu)); 137 PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user)); 138 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL)); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 144 { 145 DM cdm = dm; 146 PetscFE fe; 147 DMPolytopeType ct; 148 PetscBool simplex; 149 PetscInt dim, cStart; 150 char prefix[PETSC_MAX_PATH_LEN]; 151 152 PetscFunctionBeginUser; 153 PetscCall(DMGetDimension(dm, &dim)); 154 155 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 156 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 157 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE; 158 /* Create finite element */ 159 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 160 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe)); 161 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 162 /* Set discretization and boundary conditions for each mesh */ 163 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 164 PetscCall(DMCreateDS(dm)); 165 PetscCall((*setup)(dm, user)); 166 while (cdm) { 167 PetscCall(DMCopyDisc(dm, cdm)); 168 PetscCall(DMGetCoarseDM(cdm, &cdm)); 169 } 170 PetscCall(PetscFEDestroy(&fe)); 171 PetscFunctionReturn(0); 172 } 173 174 int main(int argc, char **argv) 175 { 176 DM dm; /* Problem specification */ 177 PetscDS ds; 178 SNES snes; /* Nonlinear solver */ 179 Vec u; /* Solutions */ 180 AppCtx user; /* User-defined work context */ 181 182 PetscFunctionBeginUser; 183 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 184 /* Primal system */ 185 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 186 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 187 PetscCall(ProcessOptions(dm, &user)); 188 PetscCall(SNESSetDM(snes, dm)); 189 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 190 PetscCall(DMCreateGlobalVector(dm, &u)); 191 PetscCall(VecSet(u, 0.0)); 192 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 193 PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 194 PetscCall(SNESSetFromOptions(snes)); 195 PetscCall(DMSNESCheckFromOptions(snes, u)); 196 197 /* Looking for field error */ 198 PetscInt Nfields; 199 PetscCall(DMGetDS(dm, &ds)); 200 PetscCall(PetscDSGetNumFields(ds, &Nfields)); 201 PetscCall(SNESSolve(snes, NULL, u)); 202 PetscCall(SNESGetSolution(snes, &u)); 203 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 204 205 /* Cleanup */ 206 PetscCall(VecDestroy(&u)); 207 PetscCall(SNESDestroy(&snes)); 208 PetscCall(DMDestroy(&dm)); 209 PetscCall(PetscFinalize()); 210 return 0; 211 } 212 213 /*TEST 214 test: 215 # L_2 convergence rate: 1.9 216 suffix: 2d_p1_conv 217 requires: triangle 218 args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu 219 test: 220 # L_2 convergence rate: 1.9 221 suffix: 2d_q1_conv 222 args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu 223 test: 224 # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5 225 suffix: 3d_p1_conv 226 requires: ctetgen 227 args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu 228 test: 229 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2 230 suffix: 3d_q1_conv 231 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 232 test: 233 # L_2 convergence rate: 1.9 234 suffix: 2d_p1_trig_conv 235 requires: triangle 236 args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig 237 test: 238 # L_2 convergence rate: 1.9 239 suffix: 2d_q1_trig_conv 240 args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig 241 test: 242 # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5 243 suffix: 3d_p1_trig_conv 244 requires: ctetgen 245 args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig 246 test: 247 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2 248 suffix: 3d_q1_trig_conv 249 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 250 test: 251 suffix: 2d_p1_gmg_vcycle 252 requires: triangle 253 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 254 -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \ 255 -mg_levels_ksp_max_it 1 \ 256 -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor 257 test: 258 suffix: 2d_p1_gmg_fcycle 259 requires: triangle 260 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 261 -ksp_type cg -ksp_rtol 1e-10 -pc_type mg -pc_mg_type full \ 262 -mg_levels_ksp_max_it 2 \ 263 -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor 264 test: 265 suffix: 2d_p1_gmg_vcycle_trig 266 requires: triangle 267 args: -exact_trig -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 268 -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \ 269 -mg_levels_ksp_max_it 1 \ 270 -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor 271 test: 272 suffix: 2d_q3_cgns 273 requires: cgns 274 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 275 TEST*/ 276