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\

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 {
134:     PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, NULL, NULL, g3_uu));
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*/
```