Actual source code: ex8.c
1: static char help[] = "Test adaptive interpolation of functions of a given polynomial order\n\n";
3: #include <petscdmplex.h>
4: #include <petscsnes.h>
6: /*
7: What properties does the adapted interpolator have?
9: 1) If we adapt to quadratics, we can get lower interpolation error for quadratics (than local interpolation) when using a linear basis
11: $ ./ex8 -dm_refine 2 -petscspace_degree 1 -qorder 1 -dim 2 -porder 2 -K 2 -num_comp 1 -use_poly 1
12: Function tests FAIL for order 2 at tolerance 1e-10 error 0.00273757
13: Function tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.0721688
14: Interpolation tests FAIL for order 2 at tolerance 1e-10 error 0.00284555
15: Interpolation tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.0721688
16: Adapting interpolator using polynomials
17: The number of input vectors 4 < 7 the maximum number of column entries
18: Interpolation poly tests FAIL for order 2 at tolerance 1e-10 error 0.00659864
19: Interpolation poly tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.0836582
20: Interpolation trig (0, 0) tests FAIL for order 1 at tolerance 1e-10 error 0.476194
21: Interpolation trig (0, 0) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22144
22: Interpolation trig (0, 1) tests FAIL for order 1 at tolerance 1e-10 error 1.39768
23: Interpolation trig (0, 1) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22144
24: Interpolation trig (1, 0) tests FAIL for order 2 at tolerance 1e-10 error 1.07315
25: Interpolation trig (1, 0) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55403
26: Interpolation trig (1, 1) tests FAIL for order 2 at tolerance 1e-10 error 1.07315
27: Interpolation trig (1, 1) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55403
29: $ ./ex8 -dm_refine 2 -petscspace_degree 1 -qorder 1 -dim 2 -porder 2 -K 3 -num_comp 1 -use_poly 1
30: Function tests FAIL for order 2 at tolerance 1e-10 error 0.00273757
31: Function tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.0721688
32: Interpolation tests FAIL for order 2 at tolerance 1e-10 error 0.00284555
33: Interpolation tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.0721688
34: Adapting interpolator using polynomials
35: The number of input vectors 6 < 7 the maximum number of column entries
36: Interpolation poly tests FAIL for order 2 at tolerance 1e-10 error 0.00194055
37: Interpolation poly tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.0525591
38: Interpolation trig (0, 0) tests FAIL for order 1 at tolerance 1e-10 error 0.476255
39: Interpolation trig (0, 0) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22132
40: Interpolation trig (0, 1) tests FAIL for order 1 at tolerance 1e-10 error 1.39785
41: Interpolation trig (0, 1) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22119
42: Interpolation trig (1, 0) tests FAIL for order 2 at tolerance 1e-10 error 1.0727
43: Interpolation trig (1, 0) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55364
44: Interpolation trig (1, 1) tests FAIL for order 2 at tolerance 1e-10 error 1.0727
45: Interpolation trig (1, 1) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55364
46: Interpolation trig (2, 0) tests FAIL for order 3 at tolerance 1e-10 error 0.705258
47: Interpolation trig (2, 0) tests FAIL for order 3 derivatives at tolerance 1e-10 error 6.82037
48: Interpolation trig (2, 1) tests FAIL for order 3 at tolerance 1e-10 error 0.705258
49: Interpolation trig (2, 1) tests FAIL for order 3 derivatives at tolerance 1e-10 error 6.82037
51: 2) We can more accurately capture low harmonics
53: If we adapt polynomials, we can be exact
55: $ ./ex8 -dm_refine 2 -petscspace_degree 1 -qorder 1 -dim 2 -porder 1 -K 2 -num_comp 1 -use_poly 1
56: Function tests pass for order 1 at tolerance 1e-10
57: Function tests pass for order 1 derivatives at tolerance 1e-10
58: Interpolation tests pass for order 1 at tolerance 1e-10
59: Interpolation tests pass for order 1 derivatives at tolerance 1e-10
60: Adapting interpolator using polynomials
61: The number of input vectors 4 < 7 the maximum number of column entries
62: Interpolation poly tests pass for order 1 at tolerance 1e-10
63: Interpolation poly tests pass for order 1 derivatives at tolerance 1e-10
64: Interpolation trig (0, 0) tests FAIL for order 1 at tolerance 1e-10 error 0.476194
65: Interpolation trig (0, 0) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22144
66: Interpolation trig (0, 1) tests FAIL for order 1 at tolerance 1e-10 error 1.39768
67: Interpolation trig (0, 1) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22144
68: Interpolation trig (1, 0) tests FAIL for order 2 at tolerance 1e-10 error 1.07315
69: Interpolation trig (1, 0) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55403
70: Interpolation trig (1, 1) tests FAIL for order 2 at tolerance 1e-10 error 1.07315
71: Interpolation trig (1, 1) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55403
73: and least for small K,
75: $ ./ex8 -dm_refine 2 -petscspace_degree 1 -qorder 1 -dim 2 -porder 1 -K 4 -num_comp 1 -use_poly 1
76: Function tests pass for order 1 at tolerance 1e-10
77: Function tests pass for order 1 derivatives at tolerance 1e-10
78: Interpolation tests pass for order 1 at tolerance 1e-10
79: Interpolation tests pass for order 1 derivatives at tolerance 1e-10
80: Adapting interpolator using polynomials
81: Interpolation poly tests FAIL for order 1 at tolerance 1e-10 error 0.0015351
82: Interpolation poly tests FAIL for order 1 derivatives at tolerance 1e-10 error 0.0427369
83: Interpolation trig (0, 0) tests FAIL for order 1 at tolerance 1e-10 error 0.476359
84: Interpolation trig (0, 0) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22115
85: Interpolation trig (0, 1) tests FAIL for order 1 at tolerance 1e-10 error 1.3981
86: Interpolation trig (0, 1) tests FAIL for order 1 derivatives at tolerance 1e-10 error 2.22087
87: Interpolation trig (1, 0) tests FAIL for order 2 at tolerance 1e-10 error 1.07228
88: Interpolation trig (1, 0) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55238
89: Interpolation trig (1, 1) tests FAIL for order 2 at tolerance 1e-10 error 1.07228
90: Interpolation trig (1, 1) tests FAIL for order 2 derivatives at tolerance 1e-10 error 4.55238
91: Interpolation trig (2, 0) tests FAIL for order 3 at tolerance 1e-10 error 0.704947
92: Interpolation trig (2, 0) tests FAIL for order 3 derivatives at tolerance 1e-10 error 6.82254
93: Interpolation trig (2, 1) tests FAIL for order 3 at tolerance 1e-10 error 0.704948
94: Interpolation trig (2, 1) tests FAIL for order 3 derivatives at tolerance 1e-10 error 6.82254
95: Interpolation trig (3, 0) tests FAIL for order 4 at tolerance 1e-10 error 0.893279
96: Interpolation trig (3, 0) tests FAIL for order 4 derivatives at tolerance 1e-10 error 8.93718
97: Interpolation trig (3, 1) tests FAIL for order 4 at tolerance 1e-10 error 0.89328
98: Interpolation trig (3, 1) tests FAIL for order 4 derivatives at tolerance 1e-10 error 8.93717
100: but adapting to harmonics gives alright polynomials errors and much better harmonics errors.
102: $ ./ex8 -dm_refine 2 -petscspace_degree 1 -qorder 1 -dim 2 -porder 1 -K 4 -num_comp 1 -use_poly 0
103: Function tests pass for order 1 at tolerance 1e-10
104: Function tests pass for order 1 derivatives at tolerance 1e-10
105: Interpolation tests pass for order 1 at tolerance 1e-10
106: Interpolation tests pass for order 1 derivatives at tolerance 1e-10
107: Adapting interpolator using harmonics
108: Interpolation poly tests FAIL for order 1 at tolerance 1e-10 error 0.0720606
109: Interpolation poly tests FAIL for order 1 derivatives at tolerance 1e-10 error 1.97779
110: Interpolation trig (0, 0) tests FAIL for order 1 at tolerance 1e-10 error 0.0398055
111: Interpolation trig (0, 0) tests FAIL for order 1 derivatives at tolerance 1e-10 error 0.995963
112: Interpolation trig (0, 1) tests FAIL for order 1 at tolerance 1e-10 error 0.0398051
113: Interpolation trig (0, 1) tests FAIL for order 1 derivatives at tolerance 1e-10 error 0.995964
114: Interpolation trig (1, 0) tests FAIL for order 2 at tolerance 1e-10 error 0.0238441
115: Interpolation trig (1, 0) tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.888611
116: Interpolation trig (1, 1) tests FAIL for order 2 at tolerance 1e-10 error 0.0238346
117: Interpolation trig (1, 1) tests FAIL for order 2 derivatives at tolerance 1e-10 error 0.888612
118: Interpolation trig (2, 0) tests FAIL for order 3 at tolerance 1e-10 error 0.0537968
119: Interpolation trig (2, 0) tests FAIL for order 3 derivatives at tolerance 1e-10 error 1.57665
120: Interpolation trig (2, 1) tests FAIL for order 3 at tolerance 1e-10 error 0.0537779
121: Interpolation trig (2, 1) tests FAIL for order 3 derivatives at tolerance 1e-10 error 1.57666
122: Interpolation trig (3, 0) tests FAIL for order 4 at tolerance 1e-10 error 0.0775838
123: Interpolation trig (3, 0) tests FAIL for order 4 derivatives at tolerance 1e-10 error 2.36926
124: Interpolation trig (3, 1) tests FAIL for order 4 at tolerance 1e-10 error 0.0775464
125: Interpolation trig (3, 1) tests FAIL for order 4 derivatives at tolerance 1e-10 error 2.36929
126: */
128: typedef struct {
129: /* Element definition */
130: PetscInt qorder; /* Order of the quadrature */
131: PetscInt Nc; /* Number of field components */
132: /* Testing space */
133: PetscInt porder; /* Order of polynomials to test */
134: PetscReal constants[3]; /* Constant values for each dimension */
135: PetscInt m; /* The frequency of sinusoids to use */
136: PetscInt dir; /* The direction of sinusoids to use */
137: /* Adaptation */
138: PetscInt K; /* Number of coarse modes used for optimization */
139: PetscBool usePoly; /* Use polynomials, or harmonics, to adapt interpolator */
140: } AppCtx;
142: typedef enum {
143: INTERPOLATION,
144: RESTRICTION,
145: INJECTION
146: } InterpType;
148: /* u = 1 */
149: PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
150: {
151: AppCtx *user = (AppCtx *)ctx;
152: PetscInt d = user->dir;
154: if (Nc > 1) {
155: for (d = 0; d < Nc; ++d) u[d] = user->constants[d];
156: } else {
157: u[0] = user->constants[d];
158: }
159: return PETSC_SUCCESS;
160: }
161: PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
162: {
163: AppCtx *user = (AppCtx *)ctx;
164: PetscInt d = user->dir;
166: if (Nc > 1) {
167: for (d = 0; d < Nc; ++d) u[d] = 0.0;
168: } else {
169: u[0] = user->constants[d];
170: }
171: return PETSC_SUCCESS;
172: }
174: /* u = x */
175: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
176: {
177: AppCtx *user = (AppCtx *)ctx;
178: PetscInt d = user->dir;
180: if (Nc > 1) {
181: for (d = 0; d < Nc; ++d) u[d] = coords[d];
182: } else {
183: u[0] = coords[d];
184: }
185: return PETSC_SUCCESS;
186: }
187: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
188: {
189: AppCtx *user = (AppCtx *)ctx;
190: PetscInt d = user->dir;
192: if (Nc > 1) {
193: for (d = 0; d < Nc; ++d) {
194: u[d] = 0.0;
195: for (PetscInt e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
196: }
197: } else {
198: u[0] = n[d];
199: }
200: return PETSC_SUCCESS;
201: }
203: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
204: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
205: {
206: AppCtx *user = (AppCtx *)ctx;
207: PetscInt d = user->dir;
209: if (Nc > 1) {
210: if (Nc > 2) {
211: u[0] = coords[0] * coords[1];
212: u[1] = coords[1] * coords[2];
213: u[2] = coords[2] * coords[0];
214: } else {
215: u[0] = coords[0] * coords[0];
216: u[1] = coords[0] * coords[1];
217: }
218: } else {
219: u[0] = coords[d] * coords[d];
220: }
221: return PETSC_SUCCESS;
222: }
223: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
224: {
225: AppCtx *user = (AppCtx *)ctx;
226: PetscInt d = user->dir;
228: if (Nc > 1) {
229: if (Nc > 2) {
230: u[0] = coords[1] * n[0] + coords[0] * n[1];
231: u[1] = coords[2] * n[1] + coords[1] * n[2];
232: u[2] = coords[2] * n[0] + coords[0] * n[2];
233: } else {
234: u[0] = 2.0 * coords[0] * n[0];
235: u[1] = coords[1] * n[0] + coords[0] * n[1];
236: }
237: } else {
238: u[0] = 2.0 * coords[d] * n[d];
239: }
240: return PETSC_SUCCESS;
241: }
243: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
244: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
245: {
246: AppCtx *user = (AppCtx *)ctx;
247: PetscInt d = user->dir;
249: if (Nc > 1) {
250: if (Nc > 2) {
251: u[0] = coords[0] * coords[0] * coords[1];
252: u[1] = coords[1] * coords[1] * coords[2];
253: u[2] = coords[2] * coords[2] * coords[0];
254: } else {
255: u[0] = coords[0] * coords[0] * coords[0];
256: u[1] = coords[0] * coords[0] * coords[1];
257: }
258: } else {
259: u[0] = coords[d] * coords[d] * coords[d];
260: }
261: return PETSC_SUCCESS;
262: }
263: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
264: {
265: AppCtx *user = (AppCtx *)ctx;
266: PetscInt d = user->dir;
268: if (Nc > 1) {
269: if (Nc > 2) {
270: u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
271: u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2];
272: u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0];
273: } else {
274: u[0] = 3.0 * coords[0] * coords[0] * n[0];
275: u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
276: }
277: } else {
278: u[0] = 3.0 * coords[d] * coords[d] * n[d];
279: }
280: return PETSC_SUCCESS;
281: }
283: /* u = x^4 or u = (x^4, x^2y^2) or u = (x^2y^2, y^2z^2, z^2x^2) */
284: PetscErrorCode quartic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
285: {
286: AppCtx *user = (AppCtx *)ctx;
287: PetscInt d = user->dir;
289: if (Nc > 1) {
290: if (Nc > 2) {
291: u[0] = coords[0] * coords[0] * coords[1] * coords[1];
292: u[1] = coords[1] * coords[1] * coords[2] * coords[2];
293: u[2] = coords[2] * coords[2] * coords[0] * coords[0];
294: } else {
295: u[0] = coords[0] * coords[0] * coords[0] * coords[0];
296: u[1] = coords[0] * coords[0] * coords[1] * coords[1];
297: }
298: } else {
299: u[0] = coords[d] * coords[d] * coords[d] * coords[d];
300: }
301: return PETSC_SUCCESS;
302: }
303: PetscErrorCode quarticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
304: {
305: AppCtx *user = (AppCtx *)ctx;
306: PetscInt d = user->dir;
308: if (Nc > 1) {
309: if (Nc > 2) {
310: u[0] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1];
311: u[1] = 2.0 * coords[1] * coords[2] * coords[2] * n[1] + 2.0 * coords[1] * coords[1] * coords[2] * n[2];
312: u[2] = 2.0 * coords[2] * coords[0] * coords[0] * n[2] + 2.0 * coords[2] * coords[2] * coords[0] * n[0];
313: } else {
314: u[0] = 4.0 * coords[0] * coords[0] * coords[0] * n[0];
315: u[1] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1];
316: }
317: } else {
318: u[0] = 4.0 * coords[d] * coords[d] * coords[d] * n[d];
319: }
320: return PETSC_SUCCESS;
321: }
323: PetscErrorCode mytanh(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
324: {
325: AppCtx *user = (AppCtx *)ctx;
326: PetscInt d = user->dir;
328: if (Nc > 1) {
329: for (d = 0; d < Nc; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
330: } else {
331: u[0] = PetscTanhReal(coords[d] - 0.5);
332: }
333: return PETSC_SUCCESS;
334: }
335: PetscErrorCode mytanhDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
336: {
337: AppCtx *user = (AppCtx *)ctx;
338: PetscInt d = user->dir;
340: if (Nc > 1) {
341: for (d = 0; d < Nc; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
342: } else {
343: u[0] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
344: }
345: return PETSC_SUCCESS;
346: }
348: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
349: {
350: AppCtx *user = (AppCtx *)ctx;
351: PetscInt m = user->m, d = user->dir;
353: if (Nc > 1) {
354: for (d = 0; d < Nc; ++d) u[d] = PetscSinReal(PETSC_PI * m * coords[d]);
355: } else {
356: u[0] = PetscSinReal(PETSC_PI * m * coords[d]);
357: }
358: return PETSC_SUCCESS;
359: }
360: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
361: {
362: AppCtx *user = (AppCtx *)ctx;
363: PetscInt m = user->m, d = user->dir;
365: if (Nc > 1) {
366: for (d = 0; d < Nc; ++d) u[d] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d];
367: } else {
368: u[0] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d];
369: }
370: return PETSC_SUCCESS;
371: }
373: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
374: {
375: PetscFunctionBeginUser;
376: options->qorder = 0;
377: options->Nc = PETSC_DEFAULT;
378: options->porder = 0;
379: options->m = 1;
380: options->dir = 0;
381: options->K = 0;
382: options->usePoly = PETSC_TRUE;
384: PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
385: PetscCall(PetscOptionsInt("-qorder", "The quadrature order", "ex8.c", options->qorder, &options->qorder, NULL));
386: PetscCall(PetscOptionsInt("-num_comp", "The number of field components", "ex8.c", options->Nc, &options->Nc, NULL));
387: PetscCall(PetscOptionsInt("-porder", "The order of polynomials to test", "ex8.c", options->porder, &options->porder, NULL));
388: PetscCall(PetscOptionsInt("-K", "The number of coarse modes used in optimization", "ex8.c", options->K, &options->K, NULL));
389: PetscCall(PetscOptionsBool("-use_poly", "Use polynomials (or harmonics) to adapt interpolator", "ex8.c", options->usePoly, &options->usePoly, NULL));
390: PetscOptionsEnd();
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
395: {
396: PetscFunctionBeginUser;
397: PetscCall(DMCreate(comm, dm));
398: PetscCall(DMSetType(*dm, DMPLEX));
399: PetscCall(DMSetFromOptions(*dm));
400: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /* Setup functions to approximate */
405: static PetscErrorCode SetupFunctions(DM dm, PetscBool usePoly, PetscInt order, PetscInt dir, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), AppCtx *user)
406: {
407: PetscInt dim;
409: PetscFunctionBeginUser;
410: user->dir = dir;
411: if (usePoly) {
412: switch (order) {
413: case 0:
414: exactFuncs[0] = constant;
415: exactFuncDers[0] = constantDer;
416: break;
417: case 1:
418: exactFuncs[0] = linear;
419: exactFuncDers[0] = linearDer;
420: break;
421: case 2:
422: exactFuncs[0] = quadratic;
423: exactFuncDers[0] = quadraticDer;
424: break;
425: case 3:
426: exactFuncs[0] = cubic;
427: exactFuncDers[0] = cubicDer;
428: break;
429: case 4:
430: exactFuncs[0] = quartic;
431: exactFuncDers[0] = quarticDer;
432: break;
433: default:
434: PetscCall(DMGetDimension(dm, &dim));
435: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order);
436: }
437: } else {
438: user->m = order;
439: exactFuncs[0] = trig;
440: exactFuncDers[0] = trigDer;
441: }
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
446: {
447: Vec u;
448: PetscReal n[3] = {1.0, 1.0, 1.0};
450: PetscFunctionBeginUser;
451: PetscCall(DMGetGlobalVector(dm, &u));
452: /* Project function into FE function space */
453: PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
454: PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
455: /* Compare approximation to exact in L_2 */
456: PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
457: PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
458: PetscCall(DMRestoreGlobalVector(dm, &u));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
463: {
464: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
465: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
466: void *exactCtxs[3];
467: MPI_Comm comm;
468: PetscReal error, errorDer, tol = PETSC_SMALL;
470: PetscFunctionBeginUser;
471: exactCtxs[0] = user;
472: exactCtxs[1] = user;
473: exactCtxs[2] = user;
474: user->constants[0] = 1.0;
475: user->constants[1] = 2.0;
476: user->constants[2] = 3.0;
477: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
478: PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user));
479: PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
480: /* Report result */
481: if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
482: else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
483: if (errorDer > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer));
484: else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: /* Compare approximation to exact in L_2 */
489: static PetscErrorCode CheckTransferError(DM fdm, PetscBool usePoly, PetscInt order, PetscInt dir, const char *testname, Vec fu, AppCtx *user)
490: {
491: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
492: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
493: PetscReal n[3] = {1.0, 1.0, 1.0};
494: void *exactCtxs[3];
495: MPI_Comm comm;
496: PetscReal error, errorDer, tol = PETSC_SMALL;
498: PetscFunctionBeginUser;
499: exactCtxs[0] = user;
500: exactCtxs[1] = user;
501: exactCtxs[2] = user;
502: user->constants[0] = 1.0;
503: user->constants[1] = 2.0;
504: user->constants[2] = 3.0;
505: PetscCall(PetscObjectGetComm((PetscObject)fdm, &comm));
506: PetscCall(SetupFunctions(fdm, usePoly, order, dir, exactFuncs, exactFuncDers, user));
507: PetscCall(DMGetCoordinatesLocalSetUp(fdm));
508: PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error));
509: PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer));
510: /* Report result */
511: if (error > tol) PetscCall(PetscPrintf(comm, "%s tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", testname, order, (double)tol, (double)error));
512: else PetscCall(PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " at tolerance %g\n", testname, order, (double)tol));
513: if (errorDer > tol) PetscCall(PetscPrintf(comm, "%s tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", testname, order, (double)tol, (double)errorDer));
514: else PetscCall(PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", testname, order, (double)tol));
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static PetscErrorCode CheckTransfer(DM dm, InterpType inType, PetscInt order, AppCtx *user)
519: {
520: PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
521: PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, PetscCtx ctx);
522: void *exactCtxs[3];
523: DM rdm = NULL, idm = NULL, fdm = NULL;
524: Mat Interp, InterpAdapt = NULL;
525: Vec iu, fu, scaling = NULL;
526: MPI_Comm comm;
527: const char *testname = "Unknown";
528: char checkname[PETSC_MAX_PATH_LEN];
530: PetscFunctionBeginUser;
531: exactCtxs[0] = exactCtxs[1] = exactCtxs[2] = user;
532: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
533: PetscCall(DMRefine(dm, comm, &rdm));
534: PetscCall(DMViewFromOptions(rdm, NULL, "-ref_dm_view"));
535: PetscCall(DMSetCoarseDM(rdm, dm));
536: PetscCall(DMCopyDisc(dm, rdm));
537: switch (inType) {
538: case INTERPOLATION:
539: testname = "Interpolation";
540: idm = dm;
541: fdm = rdm;
542: break;
543: case RESTRICTION:
544: testname = "Restriction";
545: idm = rdm;
546: fdm = dm;
547: break;
548: case INJECTION:
549: testname = "Injection";
550: idm = rdm;
551: fdm = dm;
552: break;
553: }
554: PetscCall(DMGetGlobalVector(idm, &iu));
555: PetscCall(DMGetGlobalVector(fdm, &fu));
556: PetscCall(DMSetApplicationContext(dm, user));
557: PetscCall(DMSetApplicationContext(rdm, user));
558: /* Project function into initial FE function space */
559: PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user));
560: PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu));
561: /* Interpolate function into final FE function space */
562: switch (inType) {
563: case INTERPOLATION:
564: PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
565: PetscCall(MatInterpolate(Interp, iu, fu));
566: break;
567: case RESTRICTION:
568: PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
569: PetscCall(MatRestrict(Interp, iu, fu));
570: PetscCall(VecPointwiseMult(fu, scaling, fu));
571: break;
572: case INJECTION:
573: PetscCall(DMCreateInjection(dm, rdm, &Interp));
574: PetscCall(MatRestrict(Interp, iu, fu));
575: break;
576: }
577: PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, testname, fu, user));
578: if (user->K && (inType == INTERPOLATION)) {
579: KSP smoother;
580: Mat A, iVM, fVM;
581: Vec iV, fV;
582: PetscInt k, dim, d, im, fm;
584: PetscCall(PetscPrintf(comm, " Adapting interpolator using %s\n", user->usePoly ? "polynomials" : "harmonics"));
585: PetscCall(DMGetDimension(dm, &dim));
586: /* Project coarse modes into initial and final FE function space */
587: PetscCall(DMGetGlobalVector(idm, &iV));
588: PetscCall(DMGetGlobalVector(fdm, &fV));
589: PetscCall(VecGetLocalSize(iV, &im));
590: PetscCall(VecGetLocalSize(fV, &fm));
591: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)dm), im, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &iVM));
592: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)dm), fm, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &fVM));
593: PetscCall(DMRestoreGlobalVector(idm, &iV));
594: PetscCall(DMRestoreGlobalVector(fdm, &fV));
595: for (k = 0; k < user->K; ++k) {
596: for (d = 0; d < dim; ++d) {
597: PetscCall(MatDenseGetColumnVecWrite(iVM, k * dim + d, &iV));
598: PetscCall(MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV));
599: PetscCall(SetupFunctions(idm, user->usePoly, user->usePoly ? k : k + 1, d, exactFuncs, exactFuncDers, user));
600: PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iV));
601: PetscCall(DMProjectFunction(fdm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, fV));
602: PetscCall(MatDenseRestoreColumnVecWrite(iVM, k * dim + d, &iV));
603: PetscCall(MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV));
604: }
605: }
607: /* Adapt interpolator */
608: PetscCall(DMCreateMatrix(rdm, &A));
609: PetscCall(MatShift(A, 1.0));
610: PetscCall(KSPCreate(comm, &smoother));
611: PetscCall(KSPSetFromOptions(smoother));
612: PetscCall(KSPSetOperators(smoother, A, A));
613: PetscCall(DMAdaptInterpolator(dm, rdm, Interp, smoother, fVM, iVM, &InterpAdapt, user));
614: /* Interpolate function into final FE function space */
615: PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s poly", testname));
616: PetscCall(MatInterpolate(InterpAdapt, iu, fu));
617: PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, checkname, fu, user));
618: for (k = 0; k < user->K; ++k) {
619: for (d = 0; d < dim; ++d) {
620: PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s trig (%" PetscInt_FMT ", %" PetscInt_FMT ")", testname, k, d));
621: PetscCall(MatDenseGetColumnVecRead(iVM, k * dim + d, &iV));
622: PetscCall(MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV));
623: PetscCall(MatInterpolate(InterpAdapt, iV, fV));
624: PetscCall(CheckTransferError(fdm, PETSC_FALSE, k + 1, d, checkname, fV, user));
625: PetscCall(MatDenseRestoreColumnVecRead(iVM, k * dim + d, &iV));
626: PetscCall(MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV));
627: }
628: }
629: /* Cleanup */
630: PetscCall(KSPDestroy(&smoother));
631: PetscCall(MatDestroy(&A));
632: PetscCall(MatDestroy(&InterpAdapt));
633: PetscCall(MatDestroy(&iVM));
634: PetscCall(MatDestroy(&fVM));
635: }
636: PetscCall(DMRestoreGlobalVector(idm, &iu));
637: PetscCall(DMRestoreGlobalVector(fdm, &fu));
638: PetscCall(MatDestroy(&Interp));
639: PetscCall(VecDestroy(&scaling));
640: PetscCall(DMDestroy(&rdm));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: int main(int argc, char **argv)
645: {
646: DM dm;
647: PetscFE fe;
648: AppCtx user;
649: PetscInt dim;
650: PetscBool simplex;
652: PetscFunctionBeginUser;
653: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
654: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
655: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
657: PetscCall(DMGetDimension(dm, &dim));
658: PetscCall(DMPlexIsSimplex(dm, &simplex));
659: PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.Nc < 0 ? dim : user.Nc, simplex, NULL, user.qorder, &fe));
660: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
661: PetscCall(PetscFEDestroy(&fe));
662: PetscCall(DMCreateDS(dm));
664: PetscCall(CheckFunctions(dm, user.porder, &user));
665: PetscCall(CheckTransfer(dm, INTERPOLATION, user.porder, &user));
666: PetscCall(CheckTransfer(dm, INJECTION, user.porder, &user));
667: PetscCall(DMDestroy(&dm));
668: PetscCall(PetscFinalize());
669: return 0;
670: }
672: /*TEST
674: # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34)
675: # 2D/3D P_1 on a simplex
676: test:
677: suffix: p1
678: requires: triangle ctetgen
679: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -petscspace_degree 1 -num_comp 1 -qorder 1 -porder {{1}separate output}
680: test:
681: suffix: p1_pragmatic
682: requires: triangle ctetgen pragmatic
683: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder {{1 2}separate output}
684: test:
685: suffix: p1_adapt
686: requires: triangle ctetgen
687: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -dm_refine 3 -petscspace_degree 1 -qorder 1 -porder {{1 2}separate output}
689: # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34)
690: # 2D/3D P_2 on a simplex
691: test:
692: suffix: p2
693: requires: triangle ctetgen
694: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -petscspace_degree 2 -qorder 2 -porder {{1 2 3}separate output}
695: test:
696: suffix: p2_pragmatic
697: requires: triangle ctetgen pragmatic
698: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder {{1 2 3}separate output}
700: # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34)
701: # TODO This is broken. Check ex3 which worked
702: # 2D/3D P_3 on a simplex
703: test:
704: TODO: gll Lagrange nodes break this
705: suffix: p3
706: requires: triangle ctetgen !single
707: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -petscspace_degree 3 -qorder 3 -porder {{1 2 3 4}separate output}
708: test:
709: TODO: gll Lagrange nodes break this
710: suffix: p3_pragmatic
711: requires: triangle ctetgen pragmatic !single
712: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder {{1 2 3 4}separate output}
714: # 2D/3D Q_1 on a tensor cell
715: test:
716: suffix: q1
717: args: -dm_plex_dim {{2 3}separate output} -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder {{1 2}separate output}
719: # 2D/3D Q_2 on a tensor cell
720: test:
721: suffix: q2
722: requires: !single
723: args: -dm_plex_dim {{2 3}separate output} -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder {{1 2 3}separate output}
725: # 2D/3D Q_3 on a tensor cell
726: test:
727: TODO: gll Lagrange nodes break this
728: suffix: q3
729: requires: !single
730: args: -dm_plex_dim {{2 3}separate output} -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder {{1 2 3 4}separate output}
732: # 2D/3D P_1disc on a triangle/quadrilateral
733: # TODO Missing injection functional for simplices
734: test:
735: suffix: p1d
736: requires: triangle ctetgen
737: args: -dm_plex_dim {{2}separate output} -dm_plex_box_faces 2,2,2 -dm_plex_simplex {{0}separate output} -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder {{1 2}separate output}
739: TEST*/