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 0;
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 0;
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;
94: PetscObjectGetComm((PetscObject)dm, &comm);
95: DMGetDimension(dm, &dim);
96: options->trig = PETSC_FALSE;
98: PetscOptionsBegin(comm, "", "Helmholtz Problem Options", "DMPLEX");
99: PetscOptionsBool("-exact_trig", "Use trigonometric exact solution (better for more complex finite elements)", "ex26.c", options->trig, &options->trig, NULL);
100: PetscOptionsEnd();
102: return 0;
103: }
105: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
106: {
108: DMCreate(comm, dm);
109: DMSetType(*dm, DMPLEX);
110: DMSetFromOptions(*dm);
112: PetscObjectSetName((PetscObject)*dm, "Mesh");
113: DMSetApplicationContext(*dm, user);
114: DMViewFromOptions(*dm, NULL, "-dm_view");
116: return 0;
117: }
119: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
120: {
121: PetscDS ds;
122: DMLabel label;
123: const PetscInt id = 1;
126: DMGetDS(dm, &ds);
127: DMGetLabel(dm, "marker", &label);
128: if (user->trig) {
129: PetscDSSetResidual(ds, 0, f0_trig_u, f1_u);
130: PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu);
131: PetscDSSetExactSolution(ds, 0, trig_u, user);
132: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL);
133: PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n");
134: } else {
135: PetscDSSetResidual(ds, 0, f0_quad_u, f1_u);
136: PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu);
137: PetscDSSetExactSolution(ds, 0, quad_u, user);
138: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL);
139: }
140: return 0;
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];
153: DMGetDimension(dm, &dim);
155: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
156: DMPlexGetCellType(dm, cStart, &ct);
157: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
158: /* Create finite element */
159: PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
160: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe);
161: PetscObjectSetName((PetscObject)fe, name);
162: /* Set discretization and boundary conditions for each mesh */
163: DMSetField(dm, 0, NULL, (PetscObject)fe);
164: DMCreateDS(dm);
165: (*setup)(dm, user);
166: while (cdm) {
167: DMCopyDisc(dm, cdm);
168: DMGetCoarseDM(cdm, &cdm);
169: }
170: PetscFEDestroy(&fe);
171: return 0;
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 */
183: PetscInitialize(&argc, &argv, NULL, help);
184: /* Primal system */
185: SNESCreate(PETSC_COMM_WORLD, &snes);
186: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
187: ProcessOptions(dm, &user);
188: SNESSetDM(snes, dm);
189: SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);
190: DMCreateGlobalVector(dm, &u);
191: VecSet(u, 0.0);
192: PetscObjectSetName((PetscObject)u, "potential");
193: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
194: SNESSetFromOptions(snes);
195: DMSNESCheckFromOptions(snes, u);
197: /* Looking for field error */
198: PetscInt Nfields;
199: DMGetDS(dm, &ds);
200: PetscDSGetNumFields(ds, &Nfields);
201: SNESSolve(snes, NULL, u);
202: SNESGetSolution(snes, &u);
203: VecViewFromOptions(u, NULL, "-potential_view");
205: /* Cleanup */
206: VecDestroy(&u);
207: SNESDestroy(&snes);
208: DMDestroy(&dm);
209: 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*/