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