1 static char help[] = "Poisson Problem with finite elements.\n\ 2 This example supports automatic convergence estimation for multilevel solvers\n\ 3 and solver adaptivity.\n\n\n"; 4 5 #include <petscdmplex.h> 6 #include <petscsnes.h> 7 #include <petscds.h> 8 #include <petscconvest.h> 9 10 /* Next steps: 11 12 - Show lowest eigenmodes using SLEPc code from my ex6 13 14 - Run CR example from Brannick's slides that looks like semicoarsening 15 - Show lowest modes 16 - Show CR convergence rate 17 - Show CR solution to show non-convergence 18 - Refine coarse grid around non-converged dofs 19 - Maybe use Barry's "more then Z% above the average" monitor to label bad dofs 20 - Mark coarse cells that contain bad dofs 21 - Run SBR on coarse grid 22 23 - Run Helmholtz example from Gander's writeup 24 25 - Run Low Mach example? 26 27 - Run subduction example? 28 */ 29 30 typedef struct { 31 PetscBool cr; /* Use compatible relaxation */ 32 } AppCtx; 33 34 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 35 { 36 PetscInt d; 37 u[0] = 0.0; 38 for (d = 0; d < dim; ++d) u[0] += PetscSinReal(2.0*PETSC_PI*x[d]); 39 return 0; 40 } 41 42 static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 43 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 44 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 45 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 46 { 47 PetscInt d; 48 for (d = 0; d < dim; ++d) f0[0] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]); 49 } 50 51 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 52 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 53 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 54 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 55 { 56 PetscInt d; 57 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 58 } 59 60 static void g3_uu(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, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 64 { 65 PetscInt d; 66 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 67 } 68 69 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 70 { 71 PetscErrorCode ierr; 72 73 PetscFunctionBeginUser; 74 options->cr = PETSC_FALSE; 75 76 ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr); 77 ierr = PetscOptionsBool("-cr", "Use compatible relaxarion", "ex20.c", options->cr, &options->cr, NULL);CHKERRQ(ierr); 78 ierr = PetscOptionsEnd(); 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 83 { 84 PetscErrorCode ierr; 85 86 PetscFunctionBeginUser; 87 ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 88 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 89 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 90 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 91 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 92 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 97 { 98 PetscDS prob; 99 const PetscInt id = 1; 100 PetscErrorCode ierr; 101 102 PetscFunctionBeginUser; 103 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 104 ierr = PetscDSSetResidual(prob, 0, f0_trig_u, f1_u);CHKERRQ(ierr); 105 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 106 ierr = PetscDSSetExactSolution(prob, 0, trig_u, user);CHKERRQ(ierr); 107 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) trig_u, NULL, 1, &id, user);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 112 { 113 DM cdm = dm; 114 PetscFE fe; 115 DMPolytopeType ct; 116 PetscBool simplex; 117 PetscInt dim, cStart; 118 char prefix[PETSC_MAX_PATH_LEN]; 119 PetscErrorCode ierr; 120 121 PetscFunctionBeginUser; 122 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 123 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 124 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr); 125 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 126 127 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 128 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 129 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 130 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 131 ierr = DMCreateDS(dm);CHKERRQ(ierr); 132 ierr = (*setup)(dm, user);CHKERRQ(ierr); 133 while (cdm) { 134 ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr); 135 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 136 } 137 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 /* 142 How to do CR in PETSc: 143 144 Loop over PCMG levels, coarse to fine: 145 Run smoother for 5 iterates 146 At each iterate, solve Inj u_f = u_c with LSQR to 1e-15 147 Suppose that e_k = c^k e_0, which means log e_k = log e_0 + k log c 148 Fit log of error to look at log c, the slope 149 Check R^2 for linearity (1 - square residual / variance) 150 Solve exactly 151 Prolong to next level 152 */ 153 154 int main(int argc, char **argv) 155 { 156 DM dm; /* Problem specification */ 157 SNES snes; /* Nonlinear solver */ 158 Vec u; /* Solutions */ 159 AppCtx user; /* User-defined work context */ 160 PetscErrorCode ierr; 161 162 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 163 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 164 /* Primal system */ 165 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 166 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 167 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 168 ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);CHKERRQ(ierr); 169 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 170 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 171 ierr = PetscObjectSetName((PetscObject) u, "potential");CHKERRQ(ierr); 172 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 173 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 174 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 175 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 176 ierr = VecViewFromOptions(u, NULL, "-potential_view");CHKERRQ(ierr); 177 /* Cleanup */ 178 ierr = VecDestroy(&u);CHKERRQ(ierr); 179 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 180 ierr = DMDestroy(&dm);CHKERRQ(ierr); 181 ierr = PetscFinalize(); 182 return ierr; 183 } 184 185 /*TEST 186 187 test: 188 suffix: 2d_p1_gmg_vcycle_rate 189 requires: triangle 190 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 191 -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg \ 192 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 193 -mg_levels_esteig_ksp_type cg \ 194 -mg_levels_esteig_ksp_max_it 10 \ 195 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 196 -mg_levels_pc_type jacobi 197 198 test: 199 suffix: 2d_p1_gmg_vcycle_cr 200 requires: triangle 201 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 202 -ksp_rtol 5e-10 -pc_type mg -pc_mg_adapt_cr \ 203 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned \ 204 -mg_levels_esteig_ksp_type cg \ 205 -mg_levels_esteig_ksp_max_it 10 \ 206 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 207 -mg_levels_cr_ksp_max_it 5 -mg_levels_cr_ksp_converged_rate -mg_levels_cr_ksp_converged_rate_type error 208 209 test: 210 suffix: 2d_p1_gmg_fcycle_rate 211 requires: triangle 212 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 213 -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg -pc_mg_type full \ 214 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 215 -mg_levels_esteig_ksp_type cg \ 216 -mg_levels_esteig_ksp_max_it 10 \ 217 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 218 -mg_levels_pc_type jacobi 219 test: 220 suffix: 2d_p1_gmg_vcycle_adapt_rate 221 requires: triangle bamg 222 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 223 -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg \ 224 -pc_mg_galerkin -pc_mg_adapt_interp -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \ 225 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 226 -mg_levels_esteig_ksp_type cg \ 227 -mg_levels_esteig_ksp_max_it 10 \ 228 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 229 -mg_levels_pc_type jacobi 230 test: 231 suffix: 2d_p1_scalable_rate 232 requires: triangle 233 args: -potential_petscspace_degree 1 -dm_refine 3 \ 234 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -ksp_converged_rate \ 235 -pc_type gamg \ 236 -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \ 237 -pc_gamg_coarse_eq_limit 1000 \ 238 -pc_gamg_square_graph 1 \ 239 -pc_gamg_threshold 0.05 \ 240 -pc_gamg_threshold_scale .0 \ 241 -mg_levels_ksp_type chebyshev -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 242 -mg_levels_ksp_max_it 5 \ 243 -mg_levels_esteig_ksp_type cg \ 244 -mg_levels_esteig_ksp_max_it 10 \ 245 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 246 -mg_levels_pc_type jacobi \ 247 -matptap_via scalable 248 249 TEST*/ 250