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, PetscCtx 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, PetscCtx 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: f0[0] += u[0];
57: for (PetscInt 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]);
58: }
60: 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[])
61: {
62: switch (dim) {
63: case 1:
64: f0[0] = 1.0;
65: break;
66: case 2:
67: f0[0] = 5.0;
68: break;
69: case 3:
70: f0[0] = 11.0;
71: break;
72: default:
73: f0[0] = 5.0;
74: break;
75: }
76: f0[0] += u[0];
77: for (PetscInt d = 0; d < dim; ++d) f0[0] -= (d + 1) * PetscSqr(x[d]);
78: }
80: 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[])
81: {
82: for (PetscInt d = 0; d < dim; ++d) f1[d] = u_x[d];
83: }
85: static PetscErrorCode ProcessOptions(DM dm, AppCtx *options)
86: {
87: MPI_Comm comm;
88: PetscInt dim;
90: PetscFunctionBeginUser;
91: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
92: PetscCall(DMGetDimension(dm, &dim));
93: options->trig = PETSC_FALSE;
95: PetscOptionsBegin(comm, "", "Helmholtz Problem Options", "DMPLEX");
96: PetscCall(PetscOptionsBool("-exact_trig", "Use trigonometric exact solution (better for more complex finite elements)", "ex26.c", options->trig, &options->trig, NULL));
97: PetscOptionsEnd();
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
102: {
103: PetscFunctionBeginUser;
104: PetscCall(DMCreate(comm, dm));
105: PetscCall(DMSetType(*dm, DMPLEX));
106: PetscCall(DMSetFromOptions(*dm));
108: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
109: PetscCall(DMSetApplicationContext(*dm, user));
110: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
115: {
116: PetscDS ds;
117: DMLabel label;
118: const PetscInt id = 1;
120: PetscFunctionBeginUser;
121: PetscCall(DMGetDS(dm, &ds));
122: PetscCall(DMGetLabel(dm, "marker", &label));
123: if (user->trig) {
124: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
125: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
126: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
127: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
128: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Trig Exact Solution\n"));
129: } else {
130: PetscCall(PetscDSSetResidual(ds, 0, f0_quad_u, f1_u));
131: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
132: PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user));
133: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quad_u, NULL, user, NULL));
134: }
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
139: {
140: DM cdm = dm;
141: PetscFE fe;
142: DMPolytopeType ct;
143: PetscBool simplex;
144: PetscInt dim, cStart;
145: char prefix[PETSC_MAX_PATH_LEN];
147: PetscFunctionBeginUser;
148: PetscCall(DMGetDimension(dm, &dim));
150: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
151: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
152: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
153: /* Create finite element */
154: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
155: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
156: PetscCall(PetscObjectSetName((PetscObject)fe, name));
157: /* Set discretization and boundary conditions for each mesh */
158: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
159: PetscCall(DMCreateDS(dm));
160: PetscCall((*setup)(dm, user));
161: while (cdm) {
162: PetscCall(DMCopyDisc(dm, cdm));
163: PetscCall(DMGetCoarseDM(cdm, &cdm));
164: }
165: PetscCall(PetscFEDestroy(&fe));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: int main(int argc, char **argv)
170: {
171: DM dm; /* Problem specification */
172: PetscDS ds;
173: SNES snes; /* Nonlinear solver */
174: Vec u; /* Solutions */
175: AppCtx user; /* User-defined work context */
177: PetscFunctionBeginUser;
178: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
179: /* Primal system */
180: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
181: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
182: PetscCall(ProcessOptions(dm, &user));
183: PetscCall(SNESSetDM(snes, dm));
184: PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
185: PetscCall(DMCreateGlobalVector(dm, &u));
186: PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
187: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
188: PetscCall(SNESSetFromOptions(snes));
189: PetscCall(DMSNESCheckFromOptions(snes, u));
191: /* Looking for field error */
192: PetscInt Nfields;
193: PetscCall(DMGetDS(dm, &ds));
194: PetscCall(PetscDSGetNumFields(ds, &Nfields));
195: PetscCall(SNESSolve(snes, NULL, u));
196: PetscCall(SNESGetSolution(snes, &u));
197: PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
199: /* Cleanup */
200: PetscCall(VecDestroy(&u));
201: PetscCall(SNESDestroy(&snes));
202: PetscCall(DMDestroy(&dm));
203: PetscCall(PetscFinalize());
204: return 0;
205: }
207: /*TEST
208: test:
209: # L_2 convergence rate: 1.9
210: suffix: 2d_p1_conv
211: requires: triangle
212: args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
213: test:
214: # L_2 convergence rate: 1.9
215: suffix: 2d_q1_conv
216: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu
217: test:
218: # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
219: suffix: 3d_p1_conv
220: requires: ctetgen
221: args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu
222: test:
223: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
224: suffix: 3d_q1_conv
225: 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
226: test:
227: # L_2 convergence rate: 1.9
228: suffix: 2d_p1_trig_conv
229: requires: triangle
230: args: -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
231: test:
232: # L_2 convergence rate: 1.9
233: suffix: 2d_q1_trig_conv
234: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -dm_refine 2 -convest_num_refine 3 -pc_type lu -exact_trig
235: test:
236: # Using -convest_num_refine 3 we get L_2 convergence rate: -1.5
237: suffix: 3d_p1_trig_conv
238: requires: ctetgen
239: args: -dm_plex_dim 3 -dm_refine 2 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 -pc_type lu -exact_trig
240: test:
241: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: -1.2
242: suffix: 3d_q1_trig_conv
243: 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
244: test:
245: suffix: 2d_p1_gmg_vcycle
246: requires: triangle
247: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
248: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
249: -mg_levels_ksp_max_it 1 \
250: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
251: test:
252: suffix: 2d_p1_gmg_fcycle
253: requires: triangle
254: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
255: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg -pc_mg_type full \
256: -mg_levels_ksp_max_it 2 \
257: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
258: test:
259: suffix: 2d_p1_gmg_vcycle_trig
260: requires: triangle
261: args: -exact_trig -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
262: -ksp_type cg -ksp_rtol 1e-10 -pc_type mg \
263: -mg_levels_ksp_max_it 1 \
264: -mg_levels_pc_type jacobi -snes_monitor -ksp_monitor
265: test:
266: suffix: 2d_q3_cgns
267: requires: cgns
268: 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
269: output_file: output/empty.out
270: TEST*/