Actual source code: ex3.c

  1: static const char help[] = "Tests for determining whether a new finite element works";

  3: /*
  4:   Use -interpolation_view and -l2_projection_view to look at the interpolants.
  5: */

  7: #include <petscdmplex.h>
  8: #include <petscfe.h>
  9: #include <petscds.h>
 10: #include <petscsnes.h>

 12: static void constant(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[])
 13: {
 14:   const PetscInt Nc = uOff[1] - uOff[0];
 15:   for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5.;
 16: }

 18: static void 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[])
 19: {
 20:   const PetscInt Nc = uOff[1] - uOff[0];
 21:   for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5. * x[c];
 22: }

 24: static void quadratic(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[])
 25: {
 26:   const PetscInt Nc = uOff[1] - uOff[0];
 27:   for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5. * x[c] * x[c];
 28: }

 30: static void trig(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[])
 31: {
 32:   const PetscInt Nc = uOff[1] - uOff[0];
 33:   for (PetscInt c = 0; c < Nc; ++c) f0[c] += PetscCosReal(2. * PETSC_PI * x[c]);
 34: }

 36: /*
 37:  The prime basis for the Wheeler-Yotov-Xue prism.
 38:  */
 39: static void prime(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[])
 40: {
 41:   PetscReal x = X[0], y = X[1], z = X[2], b = 1 + x + y + z;
 42:   f0[0] += b + 2.0 * x * z + 2.0 * y * z + x * y + x * x;
 43:   f0[1] += b + 2.0 * x * z + 2.0 * y * z + x * y + y * y;
 44:   f0[2] += b - 3.0 * x * z - 3.0 * y * z - 2.0 * z * z;
 45: }

 47: static const char    *names[]     = {"constant", "linear", "quadratic", "trig", "prime"};
 48: static PetscPointFunc functions[] = {constant, linear, quadratic, trig, prime};

 50: typedef struct {
 51:   PetscPointFunc exactSol;
 52:   PetscReal      shear, flatten;
 53: } AppCtx;

 55: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 56: {
 57:   char     name[PETSC_MAX_PATH_LEN] = "constant";
 58:   PetscInt Nfunc                    = PETSC_STATIC_ARRAY_LENGTH(names), i;

 61:   options->exactSol = NULL;
 62:   options->shear    = 0.;
 63:   options->flatten  = 1.;

 65:   PetscOptionsBegin(comm, "", "FE Test Options", "PETSCFE");
 66:   PetscOptionsString("-func", "Function to project into space", "", name, name, PETSC_MAX_PATH_LEN, NULL);
 67:   PetscOptionsReal("-shear", "Factor by which to shear along the x-direction", "", options->shear, &(options->shear), NULL);
 68:   PetscOptionsReal("-flatten", "Factor by which to flatten", "", options->flatten, &(options->flatten), NULL);
 69:   PetscOptionsEnd();

 71:   for (i = 0; i < Nfunc; ++i) {
 72:     PetscBool flg;

 74:     PetscStrcmp(name, names[i], &flg);
 75:     if (flg) {
 76:       options->exactSol = functions[i];
 77:       break;
 78:     }
 79:   }
 81:   return 0;
 82: }

 84: /* The exact solution is the negative of the f0 contribution */
 85: static PetscErrorCode exactSolution(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 86: {
 87:   AppCtx  *user    = (AppCtx *)ctx;
 88:   PetscInt uOff[2] = {0, Nc};

 90:   user->exactSol(dim, 1, 0, uOff, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, time, x, 0, NULL, u);
 91:   for (PetscInt c = 0; c < Nc; ++c) u[c] *= -1.;
 92:   return 0;
 93: }

 95: static void f0(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[])
 96: {
 97:   const PetscInt Nc = uOff[1] - uOff[0];
 98:   for (PetscInt c = 0; c < Nc; ++c) f0[c] += u[c];
 99: }

101: static void g0(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[])
102: {
103:   const PetscInt Nc = uOff[1] - uOff[0];
104:   for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0;
105: }

107: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
108: {
109:   PetscDS       ds;
110:   PetscWeakForm wf;

113:   DMGetDS(dm, &ds);
114:   PetscDSGetWeakForm(ds, &wf);
115:   PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0, 0, NULL);
116:   PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, user->exactSol, 0, NULL);
117:   PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0, 0, NULL, 0, NULL, 0, NULL);
118:   PetscDSSetExactSolution(ds, 0, exactSolution, user);
119:   return 0;
120: }

122: static PetscErrorCode SetupDiscretization(DM dm, const char name[], AppCtx *user)
123: {
124:   DM      cdm = dm;
125:   PetscFE fe;
126:   char    prefix[PETSC_MAX_PATH_LEN];

129:   if (name) PetscSNPrintf(prefix, sizeof(prefix), "%s_", name);
130:   DMCreateFEDefault(dm, 1, name ? prefix : NULL, -1, &fe);
131:   PetscObjectSetName((PetscObject)fe, name ? name : "Solution");
132:   /* Set discretization and boundary conditions for each mesh */
133:   DMSetField(dm, 0, NULL, (PetscObject)fe);
134:   DMCreateDS(dm);
135:   SetupProblem(dm, user);
136:   while (cdm) {
137:     DMCopyDisc(dm, cdm);
138:     DMGetCoarseDM(cdm, &cdm);
139:   }
140:   PetscFEDestroy(&fe);
141:   return 0;
142: }

144: /* This test tells us whether the given function is contained in the approximation space */
145: static PetscErrorCode CheckInterpolation(DM dm, AppCtx *user)
146: {
147:   PetscSimplePointFunc exactSol[1];
148:   void                *exactCtx[1];
149:   PetscDS              ds;
150:   Vec                  u;
151:   PetscReal            error, tol = PETSC_SMALL;
152:   MPI_Comm             comm;

155:   PetscObjectGetComm((PetscObject)dm, &comm);
156:   DMGetDS(dm, &ds);
157:   DMGetGlobalVector(dm, &u);
158:   PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]);
159:   DMProjectFunction(dm, 0.0, exactSol, exactCtx, INSERT_ALL_VALUES, u);
160:   VecViewFromOptions(u, NULL, "-interpolation_view");
161:   DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);
162:   DMRestoreGlobalVector(dm, &u);
163:   if (error > tol) PetscPrintf(comm, "Interpolation tests FAIL at tolerance %g error %g\n", (double)tol, (double)error);
164:   else PetscPrintf(comm, "Interpolation tests pass at tolerance %g\n", (double)tol);
165:   return 0;
166: }

168: /* This test tells us whether the element is unisolvent (the mass matrix has full rank), and what rate of convergence we achieve */
169: static PetscErrorCode CheckL2Projection(DM dm, AppCtx *user)
170: {
171:   PetscSimplePointFunc exactSol[1];
172:   void                *exactCtx[1];
173:   SNES                 snes;
174:   PetscDS              ds;
175:   Vec                  u;
176:   PetscReal            error, tol = PETSC_SMALL;
177:   MPI_Comm             comm;

180:   PetscObjectGetComm((PetscObject)dm, &comm);
181:   DMGetDS(dm, &ds);
182:   DMGetGlobalVector(dm, &u);
183:   PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]);
184:   SNESCreate(comm, &snes);
185:   SNESSetDM(snes, dm);
186:   VecSet(u, 0.0);
187:   PetscObjectSetName((PetscObject)u, "solution");
188:   DMPlexSetSNESLocalFEM(dm, user, user, user);
189:   SNESSetFromOptions(snes);
190:   DMSNESCheckFromOptions(snes, u);
191:   SNESSolve(snes, NULL, u);
192:   SNESDestroy(&snes);
193:   VecViewFromOptions(u, NULL, "-l2_projection_view");
194:   DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);
195:   DMRestoreGlobalVector(dm, &u);
196:   if (error > tol) PetscPrintf(comm, "L2 projection tests FAIL at tolerance %g error %g\n", (double)tol, (double)error);
197:   else PetscPrintf(comm, "L2 projection tests pass at tolerance %g\n", (double)tol);
198:   return 0;
199: }

201: /* Distorts the mesh by shearing in the x-direction and flattening, factors provided in the options. */
202: static PetscErrorCode DistortMesh(DM dm, AppCtx *user)
203: {
204:   Vec          coordinates;
205:   PetscScalar *ca;
206:   PetscInt     dE, n, i;

209:   DMGetCoordinateDim(dm, &dE);
210:   DMGetCoordinates(dm, &coordinates);
211:   VecGetLocalSize(coordinates, &n);
212:   VecGetArray(coordinates, &ca);
213:   for (i = 0; i < (n / dE); ++i) {
214:     ca[i * dE + 0] += user->shear * ca[i * dE + 0];
215:     ca[i * dE + 1] *= user->flatten;
216:   }
217:   VecRestoreArray(coordinates, &ca);
218:   return 0;
219: }

221: int main(int argc, char **argv)
222: {
223:   DM          dm;
224:   AppCtx      user;
225:   PetscMPIInt size;

228:   PetscInitialize(&argc, &argv, NULL, help);
229:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
231:   ProcessOptions(PETSC_COMM_WORLD, &user);
232:   DMCreate(PETSC_COMM_WORLD, &dm);
233:   DMSetType(dm, DMPLEX);
234:   DMSetFromOptions(dm);
235:   DistortMesh(dm, &user);
236:   DMViewFromOptions(dm, NULL, "-dm_view");
237:   SetupDiscretization(dm, NULL, &user);

239:   CheckInterpolation(dm, &user);
240:   CheckL2Projection(dm, &user);

242:   DMDestroy(&dm);
243:   PetscFinalize();
244:   return 0;
245: }

247: /*TEST

249:   testset:
250:     args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -petscspace_degree 1\
251:           -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu

253:     test:
254:       suffix: p1_0
255:       args: -func {{constant linear}}

257:     # Using -dm_refine 2 -convest_num_refine 4 gives convergence rate 2.0
258:     test:
259:       suffix: p1_1
260:       args: -func {{quadratic trig}} \
261:             -snes_convergence_estimate -convest_num_refine 2

263:   testset:
264:     requires: !complex double
265:     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
266:             -petscspace_type sum \
267:             -petscspace_variables 3 \
268:             -petscspace_components 3 \
269:             -petscspace_sum_spaces 2 \
270:             -petscspace_sum_concatenate false \
271:               -sumcomp_0_petscspace_variables 3 \
272:               -sumcomp_0_petscspace_components 3 \
273:               -sumcomp_0_petscspace_degree 1 \
274:               -sumcomp_1_petscspace_variables 3 \
275:               -sumcomp_1_petscspace_components 3 \
276:               -sumcomp_1_petscspace_type wxy \
277:             -petscdualspace_form_degree 0 \
278:             -petscdualspace_order 1 \
279:             -petscdualspace_components 3 \
280:           -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu

282:     test:
283:       suffix: wxy_0
284:       args: -func constant

286:     test:
287:       suffix: wxy_1
288:       args: -func linear

290:     test:
291:       suffix: wxy_2
292:       args: -func prime

294:     test:
295:       suffix: wxy_3
296:       args: -func linear -shear 1 -flatten 1e-5

298: TEST*/