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: PetscSimplePointFunc 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: 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;
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, user, user, 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_SUP, "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*/