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