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, void *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 0;
160: }
161: PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *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 0;
172: }
174: /* u = x */
175: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *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 0;
186: }
187: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx)
188: {
189: AppCtx *user = (AppCtx *)ctx;
190: PetscInt d = user->dir;
192: if (Nc > 1) {
193: PetscInt e;
194: for (d = 0; d < Nc; ++d) {
195: u[d] = 0.0;
196: for (e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
197: }
198: } else {
199: u[0] = n[d];
200: }
201: return 0;
202: }
204: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
205: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx)
206: {
207: AppCtx *user = (AppCtx *)ctx;
208: PetscInt d = user->dir;
210: if (Nc > 1) {
211: if (Nc > 2) {
212: u[0] = coords[0] * coords[1];
213: u[1] = coords[1] * coords[2];
214: u[2] = coords[2] * coords[0];
215: } else {
216: u[0] = coords[0] * coords[0];
217: u[1] = coords[0] * coords[1];
218: }
219: } else {
220: u[0] = coords[d] * coords[d];
221: }
222: return 0;
223: }
224: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx)
225: {
226: AppCtx *user = (AppCtx *)ctx;
227: PetscInt d = user->dir;
229: if (Nc > 1) {
230: if (Nc > 2) {
231: u[0] = coords[1] * n[0] + coords[0] * n[1];
232: u[1] = coords[2] * n[1] + coords[1] * n[2];
233: u[2] = coords[2] * n[0] + coords[0] * n[2];
234: } else {
235: u[0] = 2.0 * coords[0] * n[0];
236: u[1] = coords[1] * n[0] + coords[0] * n[1];
237: }
238: } else {
239: u[0] = 2.0 * coords[d] * n[d];
240: }
241: return 0;
242: }
244: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
245: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx)
246: {
247: AppCtx *user = (AppCtx *)ctx;
248: PetscInt d = user->dir;
250: if (Nc > 1) {
251: if (Nc > 2) {
252: u[0] = coords[0] * coords[0] * coords[1];
253: u[1] = coords[1] * coords[1] * coords[2];
254: u[2] = coords[2] * coords[2] * coords[0];
255: } else {
256: u[0] = coords[0] * coords[0] * coords[0];
257: u[1] = coords[0] * coords[0] * coords[1];
258: }
259: } else {
260: u[0] = coords[d] * coords[d] * coords[d];
261: }
262: return 0;
263: }
264: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx)
265: {
266: AppCtx *user = (AppCtx *)ctx;
267: PetscInt d = user->dir;
269: if (Nc > 1) {
270: if (Nc > 2) {
271: u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
272: u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2];
273: u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0];
274: } else {
275: u[0] = 3.0 * coords[0] * coords[0] * n[0];
276: u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
277: }
278: } else {
279: u[0] = 3.0 * coords[d] * coords[d] * n[d];
280: }
281: return 0;
282: }
284: /* u = x^4 or u = (x^4, x^2y^2) or u = (x^2y^2, y^2z^2, z^2x^2) */
285: PetscErrorCode quartic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx)
286: {
287: AppCtx *user = (AppCtx *)ctx;
288: PetscInt d = user->dir;
290: if (Nc > 1) {
291: if (Nc > 2) {
292: u[0] = coords[0] * coords[0] * coords[1] * coords[1];
293: u[1] = coords[1] * coords[1] * coords[2] * coords[2];
294: u[2] = coords[2] * coords[2] * coords[0] * coords[0];
295: } else {
296: u[0] = coords[0] * coords[0] * coords[0] * coords[0];
297: u[1] = coords[0] * coords[0] * coords[1] * coords[1];
298: }
299: } else {
300: u[0] = coords[d] * coords[d] * coords[d] * coords[d];
301: }
302: return 0;
303: }
304: PetscErrorCode quarticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx)
305: {
306: AppCtx *user = (AppCtx *)ctx;
307: PetscInt d = user->dir;
309: if (Nc > 1) {
310: if (Nc > 2) {
311: u[0] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1];
312: u[1] = 2.0 * coords[1] * coords[2] * coords[2] * n[1] + 2.0 * coords[1] * coords[1] * coords[2] * n[2];
313: u[2] = 2.0 * coords[2] * coords[0] * coords[0] * n[2] + 2.0 * coords[2] * coords[2] * coords[0] * n[0];
314: } else {
315: u[0] = 4.0 * coords[0] * coords[0] * coords[0] * n[0];
316: u[1] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1];
317: }
318: } else {
319: u[0] = 4.0 * coords[d] * coords[d] * coords[d] * n[d];
320: }
321: return 0;
322: }
324: PetscErrorCode mytanh(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx)
325: {
326: AppCtx *user = (AppCtx *)ctx;
327: PetscInt d = user->dir;
329: if (Nc > 1) {
330: for (d = 0; d < Nc; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
331: } else {
332: u[0] = PetscTanhReal(coords[d] - 0.5);
333: }
334: return 0;
335: }
336: PetscErrorCode mytanhDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx)
337: {
338: AppCtx *user = (AppCtx *)ctx;
339: PetscInt d = user->dir;
341: if (Nc > 1) {
342: for (d = 0; d < Nc; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
343: } else {
344: u[0] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
345: }
346: return 0;
347: }
349: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx)
350: {
351: AppCtx *user = (AppCtx *)ctx;
352: PetscInt m = user->m, d = user->dir;
354: if (Nc > 1) {
355: for (d = 0; d < Nc; ++d) u[d] = PetscSinReal(PETSC_PI * m * coords[d]);
356: } else {
357: u[0] = PetscSinReal(PETSC_PI * m * coords[d]);
358: }
359: return 0;
360: }
361: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx)
362: {
363: AppCtx *user = (AppCtx *)ctx;
364: PetscInt m = user->m, d = user->dir;
366: if (Nc > 1) {
367: for (d = 0; d < Nc; ++d) u[d] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d];
368: } else {
369: u[0] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d];
370: }
371: return 0;
372: }
374: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
375: {
377: options->qorder = 0;
378: options->Nc = PETSC_DEFAULT;
379: options->porder = 0;
380: options->m = 1;
381: options->dir = 0;
382: options->K = 0;
383: options->usePoly = PETSC_TRUE;
385: PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
386: PetscOptionsInt("-qorder", "The quadrature order", "ex8.c", options->qorder, &options->qorder, NULL);
387: PetscOptionsInt("-num_comp", "The number of field components", "ex8.c", options->Nc, &options->Nc, NULL);
388: PetscOptionsInt("-porder", "The order of polynomials to test", "ex8.c", options->porder, &options->porder, NULL);
389: PetscOptionsInt("-K", "The number of coarse modes used in optimization", "ex8.c", options->K, &options->K, NULL);
390: PetscOptionsBool("-use_poly", "Use polynomials (or harmonics) to adapt interpolator", "ex8.c", options->usePoly, &options->usePoly, NULL);
391: PetscOptionsEnd();
392: return 0;
393: }
395: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
396: {
398: DMCreate(comm, dm);
399: DMSetType(*dm, DMPLEX);
400: DMSetFromOptions(*dm);
401: DMViewFromOptions(*dm, NULL, "-dm_view");
402: return 0;
403: }
405: /* Setup functions to approximate */
406: 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)
407: {
408: PetscInt dim;
411: user->dir = dir;
412: if (usePoly) {
413: switch (order) {
414: case 0:
415: exactFuncs[0] = constant;
416: exactFuncDers[0] = constantDer;
417: break;
418: case 1:
419: exactFuncs[0] = linear;
420: exactFuncDers[0] = linearDer;
421: break;
422: case 2:
423: exactFuncs[0] = quadratic;
424: exactFuncDers[0] = quadraticDer;
425: break;
426: case 3:
427: exactFuncs[0] = cubic;
428: exactFuncDers[0] = cubicDer;
429: break;
430: case 4:
431: exactFuncs[0] = quartic;
432: exactFuncDers[0] = quarticDer;
433: break;
434: default:
435: DMGetDimension(dm, &dim);
436: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order);
437: }
438: } else {
439: user->m = order;
440: exactFuncs[0] = trig;
441: exactFuncDers[0] = trigDer;
442: }
443: return 0;
444: }
446: 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)
447: {
448: Vec u;
449: PetscReal n[3] = {1.0, 1.0, 1.0};
452: DMGetGlobalVector(dm, &u);
453: /* Project function into FE function space */
454: DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
455: VecViewFromOptions(u, NULL, "-projection_view");
456: /* Compare approximation to exact in L_2 */
457: DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error);
458: DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer);
459: DMRestoreGlobalVector(dm, &u);
460: return 0;
461: }
463: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
464: {
465: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
466: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
467: void *exactCtxs[3];
468: MPI_Comm comm;
469: PetscReal error, errorDer, tol = PETSC_SMALL;
472: exactCtxs[0] = user;
473: exactCtxs[1] = user;
474: exactCtxs[2] = user;
475: user->constants[0] = 1.0;
476: user->constants[1] = 2.0;
477: user->constants[2] = 3.0;
478: PetscObjectGetComm((PetscObject)dm, &comm);
479: SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user);
480: ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
481: /* Report result */
482: if (error > tol) PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error);
483: else PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol);
484: if (errorDer > tol) PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);
485: else PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol);
486: return 0;
487: }
489: /* Compare approximation to exact in L_2 */
490: static PetscErrorCode CheckTransferError(DM fdm, PetscBool usePoly, PetscInt order, PetscInt dir, const char *testname, Vec fu, AppCtx *user)
491: {
492: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
493: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
494: PetscReal n[3] = {1.0, 1.0, 1.0};
495: void *exactCtxs[3];
496: MPI_Comm comm;
497: PetscReal error, errorDer, tol = PETSC_SMALL;
500: exactCtxs[0] = user;
501: exactCtxs[1] = user;
502: exactCtxs[2] = user;
503: user->constants[0] = 1.0;
504: user->constants[1] = 2.0;
505: user->constants[2] = 3.0;
506: PetscObjectGetComm((PetscObject)fdm, &comm);
507: SetupFunctions(fdm, usePoly, order, dir, exactFuncs, exactFuncDers, user);
508: DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error);
509: DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer);
510: /* Report result */
511: if (error > tol) PetscPrintf(comm, "%s tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", testname, order, (double)tol, (double)error);
512: else PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " at tolerance %g\n", testname, order, (double)tol);
513: if (errorDer > tol) PetscPrintf(comm, "%s tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", testname, order, (double)tol, (double)errorDer);
514: else PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", testname, order, (double)tol);
515: return 0;
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, void *ctx);
521: PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *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];
531: exactCtxs[0] = exactCtxs[1] = exactCtxs[2] = user;
532: PetscObjectGetComm((PetscObject)dm, &comm);
533: DMRefine(dm, comm, &rdm);
534: DMViewFromOptions(rdm, NULL, "-ref_dm_view");
535: DMSetCoarseDM(rdm, dm);
536: 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: DMGetGlobalVector(idm, &iu);
555: DMGetGlobalVector(fdm, &fu);
556: DMSetApplicationContext(dm, user);
557: DMSetApplicationContext(rdm, user);
558: /* Project function into initial FE function space */
559: SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user);
560: 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: DMCreateInterpolation(dm, rdm, &Interp, &scaling);
565: MatInterpolate(Interp, iu, fu);
566: break;
567: case RESTRICTION:
568: DMCreateInterpolation(dm, rdm, &Interp, &scaling);
569: MatRestrict(Interp, iu, fu);
570: VecPointwiseMult(fu, scaling, fu);
571: break;
572: case INJECTION:
573: DMCreateInjection(dm, rdm, &Interp);
574: MatRestrict(Interp, iu, fu);
575: break;
576: }
577: 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: PetscPrintf(comm, " Adapting interpolator using %s\n", user->usePoly ? "polynomials" : "harmonics");
585: DMGetDimension(dm, &dim);
586: /* Project coarse modes into initial and final FE function space */
587: DMGetGlobalVector(idm, &iV);
588: DMGetGlobalVector(fdm, &fV);
589: VecGetLocalSize(iV, &im);
590: VecGetLocalSize(fV, &fm);
591: MatCreateDense(PetscObjectComm((PetscObject)dm), im, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &iVM);
592: MatCreateDense(PetscObjectComm((PetscObject)dm), fm, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &fVM);
593: DMRestoreGlobalVector(idm, &iV);
594: DMRestoreGlobalVector(fdm, &fV);
595: for (k = 0; k < user->K; ++k) {
596: for (d = 0; d < dim; ++d) {
597: MatDenseGetColumnVecWrite(iVM, k * dim + d, &iV);
598: MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV);
599: SetupFunctions(idm, user->usePoly, user->usePoly ? k : k + 1, d, exactFuncs, exactFuncDers, user);
600: DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iV);
601: DMProjectFunction(fdm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, fV);
602: MatDenseRestoreColumnVecWrite(iVM, k * dim + d, &iV);
603: MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV);
604: }
605: }
607: /* Adapt interpolator */
608: DMCreateMatrix(rdm, &A);
609: MatShift(A, 1.0);
610: KSPCreate(comm, &smoother);
611: KSPSetFromOptions(smoother);
612: KSPSetOperators(smoother, A, A);
613: DMAdaptInterpolator(dm, rdm, Interp, smoother, fVM, iVM, &InterpAdapt, user);
614: /* Interpolate function into final FE function space */
615: PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s poly", testname);
616: MatInterpolate(InterpAdapt, iu, fu);
617: 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: PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s trig (%" PetscInt_FMT ", %" PetscInt_FMT ")", testname, k, d);
621: MatDenseGetColumnVecRead(iVM, k * dim + d, &iV);
622: MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV);
623: MatInterpolate(InterpAdapt, iV, fV);
624: CheckTransferError(fdm, PETSC_FALSE, k + 1, d, checkname, fV, user);
625: MatDenseRestoreColumnVecRead(iVM, k * dim + d, &iV);
626: MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV);
627: }
628: }
629: /* Cleanup */
630: KSPDestroy(&smoother);
631: MatDestroy(&A);
632: MatDestroy(&InterpAdapt);
633: MatDestroy(&iVM);
634: MatDestroy(&fVM);
635: }
636: DMRestoreGlobalVector(idm, &iu);
637: DMRestoreGlobalVector(fdm, &fu);
638: MatDestroy(&Interp);
639: VecDestroy(&scaling);
640: DMDestroy(&rdm);
641: return 0;
642: }
644: int main(int argc, char **argv)
645: {
646: DM dm;
647: PetscFE fe;
648: AppCtx user;
649: PetscInt dim;
650: PetscBool simplex;
653: PetscInitialize(&argc, &argv, NULL, help);
654: ProcessOptions(PETSC_COMM_WORLD, &user);
655: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
657: DMGetDimension(dm, &dim);
658: DMPlexIsSimplex(dm, &simplex);
659: PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.Nc < 0 ? dim : user.Nc, simplex, NULL, user.qorder, &fe);
660: DMSetField(dm, 0, NULL, (PetscObject)fe);
661: PetscFEDestroy(&fe);
662: DMCreateDS(dm);
664: CheckFunctions(dm, user.porder, &user);
665: CheckTransfer(dm, INTERPOLATION, user.porder, &user);
666: CheckTransfer(dm, INJECTION, user.porder, &user);
667: DMDestroy(&dm);
668: 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*/