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, 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: }
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: 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
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
229: TEST*/