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