Actual source code: ex26.c
1: static char help[] = "'Good Cop' Helmholtz Problem in 2d and 3d with finite elements.\n\
2: We solve the 'Good Cop' Helmholtz problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4: This example supports automatic convergence estimation\n\
5: and coarse space adaptivity.\n\n\n";
7: /*
8: The model problem:
9: Solve "Good Cop" Helmholtz equation on the unit square: (0,1) x (0,1)
10: - \Delta u + u = f,
11: where \Delta = Laplace operator
12: Dirichlet b.c.'s on all sides
14: */
16: #include <petscdmplex.h>
17: #include <petscsnes.h>
18: #include <petscds.h>
19: #include <petscconvest.h>
21: typedef struct {
22: PetscBool trig; /* Use trig function as exact solution */
23: } AppCtx;
25: /* For Primal Problem */
26: static void g0_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 g0[])
27: {
28: PetscInt d;
29: for (d = 0; d < dim; ++d) g0[0] = 1.0;
30: }
32: 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[])
33: {
34: PetscInt d;
35: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
36: }
38: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
39: {
40: PetscInt d;
41: *u = 0.0;
42: for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
43: return PETSC_SUCCESS;
44: }
46: static PetscErrorCode quad_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
47: {
48: PetscInt d;
49: *u = 1.0;
50: for (d = 0; d < dim; ++d) *u += (d + 1) * PetscSqr(x[d]);
51: return PETSC_SUCCESS;
52: }
54: 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[])
55: {
56: PetscInt d;
57: f0[0] += u[0];
58: for (d = 0; d < dim; ++d) f0[0] -= 4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]) + PetscSinReal(2.0 * PETSC_PI * x[d]);
59: }
61: static void f0_quad_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[])
62: {
63: PetscInt d;
64: switch (dim) {
65: case 1:
66: f0[0] = 1.0;
67: break;
68: case 2:
69: f0[0] = 5.0;
70: break;
71: case 3:
72: f0[0] = 11.0;
73: break;
74: default:
75: f0[0] = 5.0;
76: break;
77: }
78: f0[0] += u[0];
79: for (d = 0; d < dim; ++d) f0[0] -= (d + 1) * PetscSqr(x[d]);
80: }
82: 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[])
83: {
84: PetscInt d;
85: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
86: }
88: static PetscErrorCode ProcessOptions(DM dm, AppCtx *options)
89: {
90: MPI_Comm comm;
91: PetscInt dim;
93: PetscFunctionBeginUser;
94: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
95: PetscCall(DMGetDimension(dm, &dim));
96: options->trig = PETSC_FALSE;
98: PetscOptionsBegin(comm, "", "Helmholtz Problem Options", "DMPLEX");
99: PetscCall(PetscOptionsBool("-exact_trig", "Use trigonometric exact solution (better for more complex finite elements)", "ex26.c", options->trig, &options->trig, NULL));
100: PetscOptionsEnd();
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
106: {
107: PetscFunctionBeginUser;
108: PetscCall(DMCreate(comm, dm));
109: PetscCall(DMSetType(*dm, DMPLEX));
110: PetscCall(DMSetFromOptions(*dm));
112: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
113: PetscCall(DMSetApplicationContext(*dm, user));
114: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
120: {
121: PetscDS ds;
122: DMLabel label;
123: const PetscInt id = 1;
125: PetscFunctionBeginUser;
126: PetscCall(DMGetDS(dm, &ds));
127: PetscCall(DMGetLabel(dm, "marker", &label));
128: if (user->trig) {
129: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
130: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
131: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
132: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL));
133: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n"));
134: } else {
135: PetscCall(PetscDSSetResidual(ds, 0, f0_quad_u, f1_u));
136: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
137: PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user));
138: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL));
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
144: {
145: DM cdm = dm;
146: PetscFE fe;
147: DMPolytopeType ct;
148: PetscBool simplex;
149: PetscInt dim, cStart;
150: char prefix[PETSC_MAX_PATH_LEN];
152: PetscFunctionBeginUser;
153: PetscCall(DMGetDimension(dm, &dim));
155: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
156: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
157: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
158: /* Create finite element */
159: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
160: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
161: PetscCall(PetscObjectSetName((PetscObject)fe, name));
162: /* Set discretization and boundary conditions for each mesh */
163: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
164: PetscCall(DMCreateDS(dm));
165: PetscCall((*setup)(dm, user));
166: while (cdm) {
167: PetscCall(DMCopyDisc(dm, cdm));
168: PetscCall(DMGetCoarseDM(cdm, &cdm));
169: }
170: PetscCall(PetscFEDestroy(&fe));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: int main(int argc, char **argv)
175: {
176: DM dm; /* Problem specification */
177: PetscDS ds;
178: SNES snes; /* Nonlinear solver */
179: Vec u; /* Solutions */
180: AppCtx user; /* User-defined work context */
182: PetscFunctionBeginUser;
183: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
184: /* Primal system */
185: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
186: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
187: PetscCall(ProcessOptions(dm, &user));
188: PetscCall(SNESSetDM(snes, dm));
189: PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
190: PetscCall(DMCreateGlobalVector(dm, &u));
191: PetscCall(VecSet(u, 0.0));
192: PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
193: PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user));
194: PetscCall(SNESSetFromOptions(snes));
195: PetscCall(DMSNESCheckFromOptions(snes, u));
197: /* Looking for field error */
198: PetscInt Nfields;
199: PetscCall(DMGetDS(dm, &ds));
200: PetscCall(PetscDSGetNumFields(ds, &Nfields));
201: PetscCall(SNESSolve(snes, NULL, u));
202: PetscCall(SNESGetSolution(snes, &u));
203: PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
205: /* Cleanup */
206: PetscCall(VecDestroy(&u));
207: PetscCall(SNESDestroy(&snes));
208: PetscCall(DMDestroy(&dm));
209: PetscCall(PetscFinalize());
210: return 0;
211: }
213: /*TEST
214: test:
215: # L_2 convergence rate: 1.9
216: suffix: 2d_p1_conv
217: requires: triangle
218: args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
219: test:
220: # L_2 convergence rate: 1.9
221: suffix: 2d_q1_conv
222: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
223: test:
224: # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
225: suffix: 3d_p1_conv
226: requires: ctetgen
227: args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu
228: test:
229: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
230: suffix: 3d_q1_conv
231: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu
232: test:
233: # L_2 convergence rate: 1.9
234: suffix: 2d_p1_trig_conv
235: requires: triangle
236: args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
237: test:
238: # L_2 convergence rate: 1.9
239: suffix: 2d_q1_trig_conv
240: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
241: test:
242: # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
243: suffix: 3d_p1_trig_conv
244: requires: ctetgen
245: args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig
246: test:
247: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
248: suffix: 3d_q1_trig_conv
249: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig
250: test:
251: suffix: 2d_p1_gmg_vcycle
252: requires: triangle
253: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
254: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
255: -mg_levels_ksp_max_it 1 \
256: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
257: test:
258: suffix: 2d_p1_gmg_fcycle
259: requires: triangle
260: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
261: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg -pc_mg_type full \
262: -mg_levels_ksp_max_it 2 \
263: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
264: test:
265: suffix: 2d_p1_gmg_vcycle_trig
266: requires: triangle
267: args: -exact_trig -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
268: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
269: -mg_levels_ksp_max_it 1 \
270: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
271: test:
272: suffix: 2d_q3_cgns
273: requires: cgns
274: args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_box_faces 3,3 -snes_view_solution cgns:sol.cgns -potential_petscspace_degree 3 -dm_coord_petscspace_degree 3
275: TEST*/