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 than 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 PETSC_SUCCESS; 40 } 41 42 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[]) 43 { 44 PetscInt d; 45 for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]); 46 } 47 48 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[]) 49 { 50 PetscInt d; 51 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 52 } 53 54 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[]) 55 { 56 PetscInt d; 57 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 58 } 59 60 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 61 { 62 PetscFunctionBeginUser; 63 options->cr = PETSC_FALSE; 64 PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX"); 65 PetscCall(PetscOptionsBool("-cr", "Use compatible relaxarion", "ex20.c", options->cr, &options->cr, NULL)); 66 PetscOptionsEnd(); 67 PetscFunctionReturn(PETSC_SUCCESS); 68 } 69 70 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 71 { 72 PetscFunctionBeginUser; 73 PetscCall(DMCreate(comm, dm)); 74 PetscCall(DMSetType(*dm, DMPLEX)); 75 PetscCall(DMSetFromOptions(*dm)); 76 PetscCall(DMSetApplicationContext(*dm, user)); 77 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 78 PetscFunctionReturn(PETSC_SUCCESS); 79 } 80 81 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 82 { 83 PetscDS ds; 84 DMLabel label; 85 const PetscInt id = 1; 86 87 PetscFunctionBeginUser; 88 PetscCall(DMGetDS(dm, &ds)); 89 PetscCall(DMGetLabel(dm, "marker", &label)); 90 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u)); 91 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 92 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user)); 93 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL)); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 98 { 99 DM cdm = dm; 100 PetscFE fe; 101 DMPolytopeType ct; 102 PetscBool simplex; 103 PetscInt dim, cStart; 104 char prefix[PETSC_MAX_PATH_LEN]; 105 106 PetscFunctionBeginUser; 107 PetscCall(DMGetDimension(dm, &dim)); 108 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 109 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 110 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE; 111 112 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 113 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe)); 114 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 115 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 116 PetscCall(DMCreateDS(dm)); 117 PetscCall((*setup)(dm, user)); 118 while (cdm) { 119 PetscCall(DMCopyDisc(dm, cdm)); 120 PetscCall(DMGetCoarseDM(cdm, &cdm)); 121 } 122 PetscCall(PetscFEDestroy(&fe)); 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 /* 127 How to do CR in PETSc: 128 129 Loop over PCMG levels, coarse to fine: 130 Run smoother for 5 iterates 131 At each iterate, solve Inj u_f = u_c with LSQR to 1e-15 132 Suppose that e_k = c^k e_0, which means log e_k = log e_0 + k log c 133 Fit log of error to look at log c, the slope 134 Check R^2 for linearity (1 - square residual / variance) 135 Solve exactly 136 Prolong to next level 137 */ 138 139 int main(int argc, char **argv) 140 { 141 DM dm; /* Problem specification */ 142 SNES snes; /* Nonlinear solver */ 143 Vec u; /* Solutions */ 144 AppCtx user; /* User-defined work context */ 145 146 PetscFunctionBeginUser; 147 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 148 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 149 /* Primal system */ 150 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 151 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 152 PetscCall(SNESSetDM(snes, dm)); 153 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 154 PetscCall(DMCreateGlobalVector(dm, &u)); 155 PetscCall(VecSet(u, 0.0)); 156 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 157 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user)); 158 PetscCall(SNESSetFromOptions(snes)); 159 PetscCall(SNESSolve(snes, NULL, u)); 160 PetscCall(SNESGetSolution(snes, &u)); 161 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 162 /* Cleanup */ 163 PetscCall(VecDestroy(&u)); 164 PetscCall(SNESDestroy(&snes)); 165 PetscCall(DMDestroy(&dm)); 166 PetscCall(PetscFinalize()); 167 return 0; 168 } 169 170 /*TEST 171 172 test: 173 suffix: 2d_p1_gmg_vcycle_rate 174 requires: triangle 175 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 176 -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg \ 177 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 178 -mg_levels_esteig_ksp_type cg \ 179 -mg_levels_esteig_ksp_max_it 10 \ 180 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 181 -mg_levels_pc_type jacobi 182 183 test: 184 suffix: 2d_p1_gmg_vcycle_cr 185 requires: triangle 186 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 187 -ksp_rtol 5e-10 -pc_type mg -pc_mg_adapt_cr \ 188 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned \ 189 -mg_levels_esteig_ksp_type cg \ 190 -mg_levels_esteig_ksp_max_it 10 \ 191 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 192 -mg_levels_cr_ksp_max_it 5 -mg_levels_cr_ksp_converged_rate -mg_levels_cr_ksp_converged_rate_type error 193 194 test: 195 suffix: 2d_p1_gmg_fcycle_rate 196 requires: triangle 197 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 198 -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg -pc_mg_type full \ 199 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 200 -mg_levels_esteig_ksp_type cg \ 201 -mg_levels_esteig_ksp_max_it 10 \ 202 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 203 -mg_levels_pc_type jacobi 204 test: 205 suffix: 2d_p1_gmg_vcycle_adapt_rate 206 requires: triangle 207 args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 208 -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg \ 209 -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \ 210 -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 211 -mg_levels_esteig_ksp_type cg \ 212 -mg_levels_esteig_ksp_max_it 10 \ 213 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \ 214 -mg_levels_pc_type jacobi 215 test: 216 suffix: 2d_p1_scalable_rate 217 requires: triangle 218 args: -potential_petscspace_degree 1 -dm_refine 3 \ 219 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -ksp_converged_rate \ 220 -pc_type gamg -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_esteig_ksp_type cg \ 221 -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \ 222 -pc_gamg_coarse_eq_limit 1000 \ 223 -pc_gamg_threshold 0.05 \ 224 -pc_gamg_threshold_scale .0 \ 225 -mg_levels_ksp_type chebyshev -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \ 226 -mg_levels_ksp_max_it 5 \ 227 -matptap_via scalable 228 229 TEST*/ 230