Actual source code: ex20.c

  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";

  5: #include <petscdmplex.h>
  6: #include <petscsnes.h>
  7: #include <petscds.h>
  8: #include <petscconvest.h>

 10: /* Next steps:

 12: - Show lowest eigenmodes using SLEPc code from my ex6

 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

 23: - Run Helmholtz example from Gander's writeup

 25: - Run Low Mach example?

 27: - Run subduction example?
 28: */

 30: typedef struct {
 31:   PetscBool cr; /* Use compatible relaxation */
 32: } AppCtx;

 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: }

 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: }

 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: }

 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: }

 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: }

 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: }

 81: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
 82: {
 83:   PetscDS        ds;
 84:   DMLabel        label;
 85:   const PetscInt id = 1;

 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: }

 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];

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;

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: }

126: /*
127:   How to do CR in PETSc:

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: */

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 */

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, &user, &user, &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: }

170: /*TEST

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

183:   test:
184:     suffix: 2d_p1_gmg_vcycle_cr
185:     TODO: broken
186:     # cannot MatShift() a MATNORMAL until this MatType inherits from MATSHELL, cf. https://gitlab.com/petsc/petsc/-/issues/972
187:     requires: triangle
188:     args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
189:           -ksp_rtol 5e-10 -pc_type mg  -pc_mg_adapt_cr \
190:             -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned \
191:             -mg_levels_esteig_ksp_type cg \
192:             -mg_levels_esteig_ksp_max_it 10 \
193:             -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \
194:             -mg_levels_cr_ksp_max_it 5 -mg_levels_cr_ksp_converged_rate -mg_levels_cr_ksp_converged_rate_type error

196:   test:
197:     suffix: 2d_p1_gmg_fcycle_rate
198:     requires: triangle
199:     args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
200:           -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg -pc_mg_type full \
201:             -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \
202:             -mg_levels_esteig_ksp_type cg \
203:             -mg_levels_esteig_ksp_max_it 10 \
204:             -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \
205:             -mg_levels_pc_type jacobi
206:   test:
207:     suffix: 2d_p1_gmg_vcycle_adapt_rate
208:     requires: triangle
209:     args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
210:           -ksp_rtol 5e-10 -ksp_converged_rate -pc_type mg \
211:             -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \
212:             -mg_levels_ksp_max_it 5 -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \
213:             -mg_levels_esteig_ksp_type cg \
214:             -mg_levels_esteig_ksp_max_it 10 \
215:             -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 \
216:             -mg_levels_pc_type jacobi
217:   test:
218:     suffix: 2d_p1_scalable_rate
219:     requires: triangle
220:     args: -potential_petscspace_degree 1 -dm_refine 3 \
221:       -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -ksp_converged_rate \
222:       -pc_type gamg -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_esteig_ksp_type cg \
223:         -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \
224:         -pc_gamg_coarse_eq_limit 1000 \
225:         -pc_gamg_threshold 0.05 \
226:         -pc_gamg_threshold_scale .0 \
227:         -mg_levels_ksp_type chebyshev -mg_levels_ksp_norm_type preconditioned -mg_levels_ksp_converged_rate \
228:         -mg_levels_ksp_max_it 5                                                \
229:       -matptap_via scalable

231: TEST*/