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*/