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