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();
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
105: {
106: PetscFunctionBeginUser;
107: PetscCall(DMCreate(comm, dm));
108: PetscCall(DMSetType(*dm, DMPLEX));
109: PetscCall(DMSetFromOptions(*dm));
111: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
112: PetscCall(DMSetApplicationContext(*dm, user));
113: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
118: {
119: PetscDS ds;
120: DMLabel label;
121: const PetscInt id = 1;
123: PetscFunctionBeginUser;
124: PetscCall(DMGetDS(dm, &ds));
125: PetscCall(DMGetLabel(dm, "marker", &label));
126: if (user->trig) {
127: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
128: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
129: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
130: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL));
131: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n"));
132: } else {
133: PetscCall(PetscDSSetResidual(ds, 0, f0_quad_u, f1_u));
134: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
135: PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user));
136: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL));
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
142: {
143: DM cdm = dm;
144: PetscFE fe;
145: DMPolytopeType ct;
146: PetscBool simplex;
147: PetscInt dim, cStart;
148: char prefix[PETSC_MAX_PATH_LEN];
150: PetscFunctionBeginUser;
151: PetscCall(DMGetDimension(dm, &dim));
153: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
154: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
155: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
156: /* Create finite element */
157: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
158: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
159: PetscCall(PetscObjectSetName((PetscObject)fe, name));
160: /* Set discretization and boundary conditions for each mesh */
161: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
162: PetscCall(DMCreateDS(dm));
163: PetscCall((*setup)(dm, user));
164: while (cdm) {
165: PetscCall(DMCopyDisc(dm, cdm));
166: PetscCall(DMGetCoarseDM(cdm, &cdm));
167: }
168: PetscCall(PetscFEDestroy(&fe));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: int main(int argc, char **argv)
173: {
174: DM dm; /* Problem specification */
175: PetscDS ds;
176: SNES snes; /* Nonlinear solver */
177: Vec u; /* Solutions */
178: AppCtx user; /* User-defined work context */
180: PetscFunctionBeginUser;
181: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
182: /* Primal system */
183: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
184: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
185: PetscCall(ProcessOptions(dm, &user));
186: PetscCall(SNESSetDM(snes, dm));
187: PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
188: PetscCall(DMCreateGlobalVector(dm, &u));
189: PetscCall(VecSet(u, 0.0));
190: PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
191: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
192: PetscCall(SNESSetFromOptions(snes));
193: PetscCall(DMSNESCheckFromOptions(snes, u));
195: /* Looking for field error */
196: PetscInt Nfields;
197: PetscCall(DMGetDS(dm, &ds));
198: PetscCall(PetscDSGetNumFields(ds, &Nfields));
199: PetscCall(SNESSolve(snes, NULL, u));
200: PetscCall(SNESGetSolution(snes, &u));
201: PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
203: /* Cleanup */
204: PetscCall(VecDestroy(&u));
205: PetscCall(SNESDestroy(&snes));
206: PetscCall(DMDestroy(&dm));
207: PetscCall(PetscFinalize());
208: return 0;
209: }
211: /*TEST
212: test:
213: # L_2 convergence rate: 1.9
214: suffix: 2d_p1_conv
215: requires: triangle
216: args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
217: test:
218: # L_2 convergence rate: 1.9
219: suffix: 2d_q1_conv
220: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
221: test:
222: # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
223: suffix: 3d_p1_conv
224: requires: ctetgen
225: args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu
226: test:
227: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
228: suffix: 3d_q1_conv
229: 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
230: test:
231: # L_2 convergence rate: 1.9
232: suffix: 2d_p1_trig_conv
233: requires: triangle
234: args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
235: test:
236: # L_2 convergence rate: 1.9
237: suffix: 2d_q1_trig_conv
238: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
239: test:
240: # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
241: suffix: 3d_p1_trig_conv
242: requires: ctetgen
243: args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig
244: test:
245: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
246: suffix: 3d_q1_trig_conv
247: 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
248: test:
249: suffix: 2d_p1_gmg_vcycle
250: requires: triangle
251: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
252: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
253: -mg_levels_ksp_max_it 1 \
254: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
255: test:
256: suffix: 2d_p1_gmg_fcycle
257: requires: triangle
258: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
259: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg -pc_mg_type full \
260: -mg_levels_ksp_max_it 2 \
261: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
262: test:
263: suffix: 2d_p1_gmg_vcycle_trig
264: requires: triangle
265: args: -exact_trig -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
266: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
267: -mg_levels_ksp_max_it 1 \
268: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
269: test:
270: suffix: 2d_q3_cgns
271: requires: cgns
272: 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
273: TEST*/