1 static char help[] = "Tests implementation of PetscSpace_Sum by solving the Poisson equations using a PetscSpace_Poly and a PetscSpace_Sum and checking that \
2 solutions agree up to machine precision.\n\n";
3
4 #include <petscdmplex.h>
5 #include <petscds.h>
6 #include <petscfe.h>
7 #include <petscsnes.h>
8 /* We are solving the system of equations:
9 * \vec{u} = -\grad{p}
10 * \div{u} = f
11 */
12
13 /* Exact solutions for linear velocity
14 \vec{u} = \vec{x};
15 p = -0.5*(\vec{x} \cdot \vec{x});
16 */
linear_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)17 static PetscErrorCode linear_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
18 {
19 PetscInt c;
20
21 for (c = 0; c < Nc; ++c) u[c] = x[c];
22 return PETSC_SUCCESS;
23 }
24
linear_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)25 static PetscErrorCode linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
26 {
27 PetscInt d;
28
29 u[0] = 0.;
30 for (d = 0; d < dim; ++d) u[0] += -0.5 * x[d] * x[d];
31 return PETSC_SUCCESS;
32 }
33
linear_divu(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)34 static PetscErrorCode linear_divu(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
35 {
36 u[0] = dim;
37 return PETSC_SUCCESS;
38 }
39
40 /* fx_v are the residual functions for the equation \vec{u} = \grad{p}. f0_v is the term <v,u>.*/
f0_v(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[])41 static void f0_v(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[])
42 {
43 PetscInt i;
44
45 for (i = 0; i < dim; ++i) f0[i] = u[uOff[0] + i];
46 }
47
48 /* f1_v is the term <v,-\grad{p}> but we integrate by parts to get <\grad{v}, -p*I> */
f1_v(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[])49 static void f1_v(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[])
50 {
51 PetscInt c;
52
53 for (c = 0; c < dim; ++c) {
54 PetscInt d;
55
56 for (d = 0; d < dim; ++d) f1[c * dim + d] = (c == d) ? -u[uOff[1]] : 0;
57 }
58 }
59
60 /* Residual function for enforcing \div{u} = f. */
f0_q_linear(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 static void f0_q_linear(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 PetscScalar rhs, divu = 0;
64 PetscInt i;
65
66 (void)linear_divu(dim, t, x, dim, &rhs, NULL);
67 for (i = 0; i < dim; ++i) divu += u_x[uOff_x[0] + i * dim + i];
68 f0[0] = divu - rhs;
69 }
70
71 /* Boundary residual. Dirichlet boundary for u means u_bdy=p*n */
f0_bd_u_linear(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[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])72 static void f0_bd_u_linear(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
73 {
74 PetscScalar pressure;
75 PetscInt d;
76
77 (void)linear_p(dim, t, x, dim, &pressure, NULL);
78 for (d = 0; d < dim; ++d) f0[d] = pressure * n[d];
79 }
80
81 /* gx_yz are the jacobian functions obtained by taking the derivative of the y residual w.r.t z*/
g0_vu(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[])82 static void g0_vu(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[])
83 {
84 PetscInt c;
85
86 for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
87 }
88
g1_qu(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 g1[])89 static void g1_qu(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 g1[])
90 {
91 PetscInt c;
92
93 for (c = 0; c < dim; ++c) g1[c * dim + c] = 1.0;
94 }
95
g2_vp(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 g2[])96 static void g2_vp(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 g2[])
97 {
98 PetscInt c;
99
100 for (c = 0; c < dim; ++c) g2[c * dim + c] = -1.0;
101 }
102
103 typedef struct {
104 PetscInt dummy;
105 } AppCtx;
106
CreateMesh(MPI_Comm comm,AppCtx * ctx,DM * mesh)107 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *mesh)
108 {
109 PetscFunctionBegin;
110 PetscCall(DMCreate(comm, mesh));
111 PetscCall(DMSetType(*mesh, DMPLEX));
112 PetscCall(DMSetFromOptions(*mesh));
113 PetscCall(DMSetApplicationContext(*mesh, ctx));
114 PetscCall(DMViewFromOptions(*mesh, NULL, "-dm_view"));
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /* Setup the system of equations that we wish to solve */
SetupProblem(DM dm,AppCtx * ctx)119 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
120 {
121 PetscDS ds;
122 DMLabel label;
123 PetscWeakForm wf;
124 const PetscInt id = 1;
125 PetscInt bd;
126
127 PetscFunctionBegin;
128 PetscCall(DMGetDS(dm, &ds));
129 /* All of these are independent of the user's choice of solution */
130 PetscCall(PetscDSSetResidual(ds, 0, f0_v, f1_v));
131 PetscCall(PetscDSSetResidual(ds, 1, f0_q_linear, NULL));
132 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_vu, NULL, NULL, NULL));
133 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_vp, NULL));
134 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_qu, NULL, NULL));
135
136 PetscCall(DMGetLabel(dm, "marker", &label));
137 PetscCall(PetscDSAddBoundary(ds, DM_BC_NATURAL, "Boundary Integral", label, 1, &id, 0, 0, NULL, (PetscFortranCallbackFn *)NULL, NULL, ctx, &bd));
138 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
139 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_linear, 0, NULL));
140
141 PetscCall(PetscDSSetExactSolution(ds, 0, linear_u, NULL));
142 PetscCall(PetscDSSetExactSolution(ds, 1, linear_p, NULL));
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
146 /* Create the finite element spaces we will use for this system */
SetupDiscretization(DM mesh,DM mesh_sum,PetscErrorCode (* setup)(DM,AppCtx *),AppCtx * ctx)147 static PetscErrorCode SetupDiscretization(DM mesh, DM mesh_sum, PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *ctx)
148 {
149 DM cdm = mesh, cdm_sum = mesh_sum;
150 PetscDS ds;
151 PetscFE u, divu, u_sum, divu_sum;
152 PetscInt dim;
153 PetscBool simplex;
154
155 PetscFunctionBegin;
156 PetscCall(DMGetDimension(mesh, &dim));
157 PetscCall(DMPlexIsSimplex(mesh, &simplex));
158
159 {
160 PetscBool force;
161 // Turn off automatic quadrature selection as a test
162 PetscCall(DMGetDS(mesh_sum, &ds));
163 PetscCall(PetscDSGetForceQuad(ds, &force));
164 if (force) PetscCall(PetscDSSetForceQuad(ds, PETSC_FALSE));
165 }
166
167 /* Create FE objects and give them names so that options can be set from
168 * command line. Each field gets 2 instances (i.e. velocity and velocity_sum)created twice so that we can compare between approaches. */
169 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, dim, simplex, "velocity_", -1, &u));
170 PetscCall(PetscObjectSetName((PetscObject)u, "velocity"));
171 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum), dim, dim, simplex, "velocity_sum_", -1, &u_sum));
172 PetscCall(PetscObjectSetName((PetscObject)u_sum, "velocity_sum"));
173 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "divu_", -1, &divu));
174 PetscCall(PetscObjectSetName((PetscObject)divu, "divu"));
175 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum), dim, 1, simplex, "divu_sum_", -1, &divu_sum));
176 PetscCall(PetscObjectSetName((PetscObject)divu_sum, "divu_sum"));
177
178 PetscCall(PetscFECopyQuadrature(u, divu));
179 PetscCall(PetscFECopyQuadrature(u_sum, divu_sum));
180
181 /* Associate the FE objects with the mesh and setup the system */
182 PetscCall(DMSetField(mesh, 0, NULL, (PetscObject)u));
183 PetscCall(DMSetField(mesh, 1, NULL, (PetscObject)divu));
184 PetscCall(DMCreateDS(mesh));
185 PetscCall((*setup)(mesh, ctx));
186
187 PetscCall(DMSetField(mesh_sum, 0, NULL, (PetscObject)u_sum));
188 PetscCall(DMSetField(mesh_sum, 1, NULL, (PetscObject)divu_sum));
189 PetscCall(DMCreateDS(mesh_sum));
190 PetscCall((*setup)(mesh_sum, ctx));
191
192 while (cdm) {
193 PetscCall(DMCopyDisc(mesh, cdm));
194 PetscCall(DMGetCoarseDM(cdm, &cdm));
195 }
196
197 while (cdm_sum) {
198 PetscCall(DMCopyDisc(mesh_sum, cdm_sum));
199 PetscCall(DMGetCoarseDM(cdm_sum, &cdm_sum));
200 }
201
202 /* The Mesh now owns the fields, so we can destroy the FEs created in this
203 * function */
204 PetscCall(PetscFEDestroy(&u));
205 PetscCall(PetscFEDestroy(&divu));
206 PetscCall(PetscFEDestroy(&u_sum));
207 PetscCall(PetscFEDestroy(&divu_sum));
208 PetscCall(DMDestroy(&cdm));
209 PetscCall(DMDestroy(&cdm_sum));
210 PetscFunctionReturn(PETSC_SUCCESS);
211 }
212
main(int argc,char ** argv)213 int main(int argc, char **argv)
214 {
215 AppCtx ctx;
216 DM dm, dm_sum;
217 SNES snes, snes_sum;
218 Vec u, u_sum;
219 PetscReal errNorm;
220 const PetscReal errTol = PETSC_SMALL;
221
222 PetscFunctionBeginUser;
223 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
224
225 /* Set up a snes for the standard approach, one space with 2 components */
226 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
227 PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
228 PetscCall(SNESSetDM(snes, dm));
229
230 /* Set up a snes for the sum space approach, where each subspace of the sum space represents one component */
231 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_sum));
232 PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm_sum));
233 PetscCall(SNESSetDM(snes_sum, dm_sum));
234 PetscCall(SetupDiscretization(dm, dm_sum, SetupProblem, &ctx));
235
236 /* Set up and solve the system using standard approach. */
237 PetscCall(DMCreateGlobalVector(dm, &u));
238 PetscCall(VecSet(u, 0.0));
239 PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
240 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &ctx));
241 PetscCall(SNESSetFromOptions(snes));
242 PetscCall(DMSNESCheckFromOptions(snes, u));
243 PetscCall(SNESSolve(snes, NULL, u));
244 PetscCall(SNESGetSolution(snes, &u));
245 PetscCall(VecViewFromOptions(u, NULL, "-solution_view"));
246
247 /* Set up and solve the sum space system */
248 PetscCall(DMCreateGlobalVector(dm_sum, &u_sum));
249 PetscCall(VecSet(u_sum, 0.0));
250 PetscCall(PetscObjectSetName((PetscObject)u_sum, "solution_sum"));
251 PetscCall(DMPlexSetSNESLocalFEM(dm_sum, PETSC_FALSE, &ctx));
252 PetscCall(SNESSetFromOptions(snes_sum));
253 PetscCall(DMSNESCheckFromOptions(snes_sum, u_sum));
254 PetscCall(SNESSolve(snes_sum, NULL, u_sum));
255 PetscCall(SNESGetSolution(snes_sum, &u_sum));
256 PetscCall(VecViewFromOptions(u_sum, NULL, "-solution_sum_view"));
257
258 /* Check if standard solution and sum space solution match to machine precision */
259 PetscCall(VecAXPY(u_sum, -1, u));
260 PetscCall(VecNorm(u_sum, NORM_2, &errNorm));
261 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Sum space provides the same solution as a regular space: %s", (errNorm < errTol) ? "true" : "false"));
262
263 /* Cleanup */
264 PetscCall(VecDestroy(&u_sum));
265 PetscCall(VecDestroy(&u));
266 PetscCall(SNESDestroy(&snes_sum));
267 PetscCall(SNESDestroy(&snes));
268 PetscCall(DMDestroy(&dm_sum));
269 PetscCall(DMDestroy(&dm));
270 PetscCall(PetscFinalize());
271 return 0;
272 }
273
274 /*TEST
275 test:
276 suffix: 2d_lagrange
277 requires: triangle
278 args: -velocity_petscspace_degree 1 \
279 -velocity_petscspace_type poly \
280 -velocity_petscspace_components 2\
281 -velocity_petscdualspace_type lagrange \
282 -divu_petscspace_degree 0 \
283 -divu_petscspace_type poly \
284 -divu_petscdualspace_lagrange_continuity false \
285 -velocity_sum_petscfe_default_quadrature_order 1 \
286 -velocity_sum_petscspace_degree 1 \
287 -velocity_sum_petscspace_type sum \
288 -velocity_sum_petscspace_variables 2 \
289 -velocity_sum_petscspace_components 2 \
290 -velocity_sum_petscspace_sum_spaces 2 \
291 -velocity_sum_petscspace_sum_concatenate true \
292 -velocity_sum_petscdualspace_type lagrange \
293 -velocity_sum_sumcomp_0_petscspace_type poly \
294 -velocity_sum_sumcomp_0_petscspace_degree 1 \
295 -velocity_sum_sumcomp_0_petscspace_variables 2 \
296 -velocity_sum_sumcomp_0_petscspace_components 1 \
297 -velocity_sum_sumcomp_1_petscspace_type poly \
298 -velocity_sum_sumcomp_1_petscspace_degree 1 \
299 -velocity_sum_sumcomp_1_petscspace_variables 2 \
300 -velocity_sum_sumcomp_1_petscspace_components 1 \
301 -divu_sum_petscspace_degree 0 \
302 -divu_sum_petscspace_type sum \
303 -divu_sum_petscspace_variables 2 \
304 -divu_sum_petscspace_components 1 \
305 -divu_sum_petscspace_sum_spaces 1 \
306 -divu_sum_petscspace_sum_concatenate true \
307 -divu_sum_petscdualspace_lagrange_continuity false \
308 -divu_sum_sumcomp_0_petscspace_type poly \
309 -divu_sum_sumcomp_0_petscspace_degree 0 \
310 -divu_sum_sumcomp_0_petscspace_variables 2 \
311 -divu_sum_sumcomp_0_petscspace_components 1 \
312 -dm_refine 0 \
313 -snes_error_if_not_converged \
314 -ksp_rtol 1e-10 \
315 -ksp_error_if_not_converged \
316 -pc_type fieldsplit\
317 -pc_fieldsplit_type schur\
318 -divu_sum_petscdualspace_lagrange_continuity false \
319 -pc_fieldsplit_schur_precondition full
320 TEST*/
321