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;

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

 65:   PetscOptionsBegin(comm, "", "FE Test Options", "PETSCFE");
 66:   PetscCall(PetscOptionsString("-func", "Function to project into space", "", name, name, PETSC_MAX_PATH_LEN, NULL));
 67:   PetscCall(PetscOptionsReal("-shear", "Factor by which to shear along the x-direction", "", options->shear, &options->shear, NULL));
 68:   PetscCall(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:     PetscCall(PetscStrcmp(name, names[i], &flg));
 75:     if (flg) {
 76:       options->exactSol = functions[i];
 77:       break;
 78:     }
 79:   }
 80:   PetscCheck(options->exactSol, comm, PETSC_ERR_ARG_WRONG, "Invalid test function %s", name);
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 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 PETSC_SUCCESS;
 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;

112:   PetscFunctionBeginUser;
113:   PetscCall(DMGetDS(dm, &ds));
114:   PetscCall(PetscDSGetWeakForm(ds, &wf));
115:   PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0, 0, NULL));
116:   PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, user->exactSol, 0, NULL));
117:   PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0, 0, NULL, 0, NULL, 0, NULL));
118:   PetscCall(PetscDSSetExactSolution(ds, 0, exactSolution, user));
119:   PetscFunctionReturn(PETSC_SUCCESS);
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];

128:   PetscFunctionBeginUser;
129:   if (name) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%s_", name));
130:   PetscCall(DMCreateFEDefault(dm, 1, name ? prefix : NULL, -1, &fe));
131:   PetscCall(PetscObjectSetName((PetscObject)fe, name ? name : "Solution"));
132:   /* Set discretization and boundary conditions for each mesh */
133:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
134:   PetscCall(DMCreateDS(dm));
135:   PetscCall(SetupProblem(dm, user));
136:   while (cdm) {
137:     PetscCall(DMCopyDisc(dm, cdm));
138:     PetscCall(DMGetCoarseDM(cdm, &cdm));
139:   }
140:   PetscCall(PetscFEDestroy(&fe));
141:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscSimplePointFn *exactSol[1];
148:   void               *exactCtx[1];
149:   PetscDS             ds;
150:   Vec                 u;
151:   PetscReal           error, tol = PETSC_SMALL;
152:   MPI_Comm            comm;

154:   PetscFunctionBeginUser;
155:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
156:   PetscCall(DMGetDS(dm, &ds));
157:   PetscCall(DMGetGlobalVector(dm, &u));
158:   PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]));
159:   PetscCall(DMProjectFunction(dm, 0.0, exactSol, exactCtx, INSERT_ALL_VALUES, u));
160:   PetscCall(VecViewFromOptions(u, NULL, "-interpolation_view"));
161:   PetscCall(DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error));
162:   PetscCall(DMRestoreGlobalVector(dm, &u));
163:   if (error > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL at tolerance %g error %g\n", (double)tol, (double)error));
164:   else PetscCall(PetscPrintf(comm, "Interpolation tests pass at tolerance %g\n", (double)tol));
165:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscSimplePointFn *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;

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

208:   PetscFunctionBeginUser;
209:   PetscCall(DMGetCoordinateDim(dm, &dE));
210:   PetscCall(DMGetCoordinates(dm, &coordinates));
211:   PetscCall(VecGetLocalSize(coordinates, &n));
212:   PetscCall(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:   PetscCall(VecRestoreArray(coordinates, &ca));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

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

227:   PetscFunctionBeginUser;
228:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
229:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
230:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only.");
231:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
232:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
233:   PetscCall(DMSetType(dm, DMPLEX));
234:   PetscCall(DMSetFromOptions(dm));
235:   PetscCall(DistortMesh(dm, &user));
236:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
237:   PetscCall(SetupDiscretization(dm, NULL, &user));

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

242:   PetscCall(DMDestroy(&dm));
243:   PetscCall(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 -dm_refine 1

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