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