1 static char help[] = "Test adaptive interpolation of functions of a given polynomial order\n\n";
2
3 #include <petscdmplex.h>
4 #include <petscsnes.h>
5
6 /*
7 What properties does the adapted interpolator have?
8
9 1) If we adapt to quadratics, we can get lower interpolation error for quadratics (than local interpolation) when using a linear basis
10
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
28
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
50
51 2) We can more accurately capture low harmonics
52
53 If we adapt polynomials, we can be exact
54
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
72
73 and least for small K,
74
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
99
100 but adapting to harmonics gives alright polynomials errors and much better harmonics errors.
101
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 */
127
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;
141
142 typedef enum {
143 INTERPOLATION,
144 RESTRICTION,
145 INJECTION
146 } InterpType;
147
148 /* u = 1 */
constant(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)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;
153
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 }
constantDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)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;
165
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 }
173
174 /* u = x */
linear(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)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;
179
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 }
linearDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)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;
191
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 PETSC_SUCCESS;
202 }
203
204 /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
quadratic(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)205 PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
206 {
207 AppCtx *user = (AppCtx *)ctx;
208 PetscInt d = user->dir;
209
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 PETSC_SUCCESS;
223 }
quadraticDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)224 PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
225 {
226 AppCtx *user = (AppCtx *)ctx;
227 PetscInt d = user->dir;
228
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 PETSC_SUCCESS;
242 }
243
244 /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
cubic(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)245 PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
246 {
247 AppCtx *user = (AppCtx *)ctx;
248 PetscInt d = user->dir;
249
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 PETSC_SUCCESS;
263 }
cubicDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)264 PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
265 {
266 AppCtx *user = (AppCtx *)ctx;
267 PetscInt d = user->dir;
268
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 PETSC_SUCCESS;
282 }
283
284 /* u = x^4 or u = (x^4, x^2y^2) or u = (x^2y^2, y^2z^2, z^2x^2) */
quartic(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)285 PetscErrorCode quartic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
286 {
287 AppCtx *user = (AppCtx *)ctx;
288 PetscInt d = user->dir;
289
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 PETSC_SUCCESS;
303 }
quarticDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)304 PetscErrorCode quarticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
305 {
306 AppCtx *user = (AppCtx *)ctx;
307 PetscInt d = user->dir;
308
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 PETSC_SUCCESS;
322 }
323
mytanh(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)324 PetscErrorCode mytanh(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
325 {
326 AppCtx *user = (AppCtx *)ctx;
327 PetscInt d = user->dir;
328
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 PETSC_SUCCESS;
335 }
mytanhDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)336 PetscErrorCode mytanhDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
337 {
338 AppCtx *user = (AppCtx *)ctx;
339 PetscInt d = user->dir;
340
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 PETSC_SUCCESS;
347 }
348
trig(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)349 PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
350 {
351 AppCtx *user = (AppCtx *)ctx;
352 PetscInt m = user->m, d = user->dir;
353
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 PETSC_SUCCESS;
360 }
trigDer(PetscInt dim,PetscReal time,const PetscReal coords[],const PetscReal n[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)361 PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
362 {
363 AppCtx *user = (AppCtx *)ctx;
364 PetscInt m = user->m, d = user->dir;
365
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 PETSC_SUCCESS;
372 }
373
ProcessOptions(MPI_Comm comm,AppCtx * options)374 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
375 {
376 PetscFunctionBeginUser;
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;
384
385 PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
386 PetscCall(PetscOptionsInt("-qorder", "The quadrature order", "ex8.c", options->qorder, &options->qorder, NULL));
387 PetscCall(PetscOptionsInt("-num_comp", "The number of field components", "ex8.c", options->Nc, &options->Nc, NULL));
388 PetscCall(PetscOptionsInt("-porder", "The order of polynomials to test", "ex8.c", options->porder, &options->porder, NULL));
389 PetscCall(PetscOptionsInt("-K", "The number of coarse modes used in optimization", "ex8.c", options->K, &options->K, NULL));
390 PetscCall(PetscOptionsBool("-use_poly", "Use polynomials (or harmonics) to adapt interpolator", "ex8.c", options->usePoly, &options->usePoly, NULL));
391 PetscOptionsEnd();
392 PetscFunctionReturn(PETSC_SUCCESS);
393 }
394
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)395 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
396 {
397 PetscFunctionBeginUser;
398 PetscCall(DMCreate(comm, dm));
399 PetscCall(DMSetType(*dm, DMPLEX));
400 PetscCall(DMSetFromOptions(*dm));
401 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
402 PetscFunctionReturn(PETSC_SUCCESS);
403 }
404
405 /* Setup functions to approximate */
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 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;
409
410 PetscFunctionBeginUser;
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 PetscCall(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 PetscFunctionReturn(PETSC_SUCCESS);
444 }
445
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 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};
450
451 PetscFunctionBeginUser;
452 PetscCall(DMGetGlobalVector(dm, &u));
453 /* Project function into FE function space */
454 PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
455 PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
456 /* Compare approximation to exact in L_2 */
457 PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
458 PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
459 PetscCall(DMRestoreGlobalVector(dm, &u));
460 PetscFunctionReturn(PETSC_SUCCESS);
461 }
462
CheckFunctions(DM dm,PetscInt order,AppCtx * user)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, PetscCtx ctx);
466 PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
467 void *exactCtxs[3];
468 MPI_Comm comm;
469 PetscReal error, errorDer, tol = PETSC_SMALL;
470
471 PetscFunctionBeginUser;
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 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
479 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user));
480 PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
481 /* Report result */
482 if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
483 else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
484 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));
485 else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
486 PetscFunctionReturn(PETSC_SUCCESS);
487 }
488
489 /* Compare approximation to exact in L_2 */
CheckTransferError(DM fdm,PetscBool usePoly,PetscInt order,PetscInt dir,const char * testname,Vec fu,AppCtx * user)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, PetscCtx ctx);
493 PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, PetscCtx 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;
498
499 PetscFunctionBeginUser;
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 PetscCall(PetscObjectGetComm((PetscObject)fdm, &comm));
507 PetscCall(SetupFunctions(fdm, usePoly, order, dir, exactFuncs, exactFuncDers, user));
508 PetscCall(DMGetCoordinatesLocalSetUp(fdm));
509 PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error));
510 PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer));
511 /* Report result */
512 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));
513 else PetscCall(PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " at tolerance %g\n", testname, order, (double)tol));
514 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));
515 else PetscCall(PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", testname, order, (double)tol));
516 PetscFunctionReturn(PETSC_SUCCESS);
517 }
518
CheckTransfer(DM dm,InterpType inType,PetscInt order,AppCtx * user)519 static PetscErrorCode CheckTransfer(DM dm, InterpType inType, PetscInt order, AppCtx *user)
520 {
521 PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
522 PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, PetscCtx ctx);
523 void *exactCtxs[3];
524 DM rdm = NULL, idm = NULL, fdm = NULL;
525 Mat Interp, InterpAdapt = NULL;
526 Vec iu, fu, scaling = NULL;
527 MPI_Comm comm;
528 const char *testname = "Unknown";
529 char checkname[PETSC_MAX_PATH_LEN];
530
531 PetscFunctionBeginUser;
532 exactCtxs[0] = exactCtxs[1] = exactCtxs[2] = user;
533 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
534 PetscCall(DMRefine(dm, comm, &rdm));
535 PetscCall(DMViewFromOptions(rdm, NULL, "-ref_dm_view"));
536 PetscCall(DMSetCoarseDM(rdm, dm));
537 PetscCall(DMCopyDisc(dm, rdm));
538 switch (inType) {
539 case INTERPOLATION:
540 testname = "Interpolation";
541 idm = dm;
542 fdm = rdm;
543 break;
544 case RESTRICTION:
545 testname = "Restriction";
546 idm = rdm;
547 fdm = dm;
548 break;
549 case INJECTION:
550 testname = "Injection";
551 idm = rdm;
552 fdm = dm;
553 break;
554 }
555 PetscCall(DMGetGlobalVector(idm, &iu));
556 PetscCall(DMGetGlobalVector(fdm, &fu));
557 PetscCall(DMSetApplicationContext(dm, user));
558 PetscCall(DMSetApplicationContext(rdm, user));
559 /* Project function into initial FE function space */
560 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user));
561 PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu));
562 /* Interpolate function into final FE function space */
563 switch (inType) {
564 case INTERPOLATION:
565 PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
566 PetscCall(MatInterpolate(Interp, iu, fu));
567 break;
568 case RESTRICTION:
569 PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
570 PetscCall(MatRestrict(Interp, iu, fu));
571 PetscCall(VecPointwiseMult(fu, scaling, fu));
572 break;
573 case INJECTION:
574 PetscCall(DMCreateInjection(dm, rdm, &Interp));
575 PetscCall(MatRestrict(Interp, iu, fu));
576 break;
577 }
578 PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, testname, fu, user));
579 if (user->K && (inType == INTERPOLATION)) {
580 KSP smoother;
581 Mat A, iVM, fVM;
582 Vec iV, fV;
583 PetscInt k, dim, d, im, fm;
584
585 PetscCall(PetscPrintf(comm, " Adapting interpolator using %s\n", user->usePoly ? "polynomials" : "harmonics"));
586 PetscCall(DMGetDimension(dm, &dim));
587 /* Project coarse modes into initial and final FE function space */
588 PetscCall(DMGetGlobalVector(idm, &iV));
589 PetscCall(DMGetGlobalVector(fdm, &fV));
590 PetscCall(VecGetLocalSize(iV, &im));
591 PetscCall(VecGetLocalSize(fV, &fm));
592 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)dm), im, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &iVM));
593 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)dm), fm, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &fVM));
594 PetscCall(DMRestoreGlobalVector(idm, &iV));
595 PetscCall(DMRestoreGlobalVector(fdm, &fV));
596 for (k = 0; k < user->K; ++k) {
597 for (d = 0; d < dim; ++d) {
598 PetscCall(MatDenseGetColumnVecWrite(iVM, k * dim + d, &iV));
599 PetscCall(MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV));
600 PetscCall(SetupFunctions(idm, user->usePoly, user->usePoly ? k : k + 1, d, exactFuncs, exactFuncDers, user));
601 PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iV));
602 PetscCall(DMProjectFunction(fdm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, fV));
603 PetscCall(MatDenseRestoreColumnVecWrite(iVM, k * dim + d, &iV));
604 PetscCall(MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV));
605 }
606 }
607
608 /* Adapt interpolator */
609 PetscCall(DMCreateMatrix(rdm, &A));
610 PetscCall(MatShift(A, 1.0));
611 PetscCall(KSPCreate(comm, &smoother));
612 PetscCall(KSPSetFromOptions(smoother));
613 PetscCall(KSPSetOperators(smoother, A, A));
614 PetscCall(DMAdaptInterpolator(dm, rdm, Interp, smoother, fVM, iVM, &InterpAdapt, user));
615 /* Interpolate function into final FE function space */
616 PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s poly", testname));
617 PetscCall(MatInterpolate(InterpAdapt, iu, fu));
618 PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, checkname, fu, user));
619 for (k = 0; k < user->K; ++k) {
620 for (d = 0; d < dim; ++d) {
621 PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s trig (%" PetscInt_FMT ", %" PetscInt_FMT ")", testname, k, d));
622 PetscCall(MatDenseGetColumnVecRead(iVM, k * dim + d, &iV));
623 PetscCall(MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV));
624 PetscCall(MatInterpolate(InterpAdapt, iV, fV));
625 PetscCall(CheckTransferError(fdm, PETSC_FALSE, k + 1, d, checkname, fV, user));
626 PetscCall(MatDenseRestoreColumnVecRead(iVM, k * dim + d, &iV));
627 PetscCall(MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV));
628 }
629 }
630 /* Cleanup */
631 PetscCall(KSPDestroy(&smoother));
632 PetscCall(MatDestroy(&A));
633 PetscCall(MatDestroy(&InterpAdapt));
634 PetscCall(MatDestroy(&iVM));
635 PetscCall(MatDestroy(&fVM));
636 }
637 PetscCall(DMRestoreGlobalVector(idm, &iu));
638 PetscCall(DMRestoreGlobalVector(fdm, &fu));
639 PetscCall(MatDestroy(&Interp));
640 PetscCall(VecDestroy(&scaling));
641 PetscCall(DMDestroy(&rdm));
642 PetscFunctionReturn(PETSC_SUCCESS);
643 }
644
main(int argc,char ** argv)645 int main(int argc, char **argv)
646 {
647 DM dm;
648 PetscFE fe;
649 AppCtx user;
650 PetscInt dim;
651 PetscBool simplex;
652
653 PetscFunctionBeginUser;
654 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
655 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
656 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
657
658 PetscCall(DMGetDimension(dm, &dim));
659 PetscCall(DMPlexIsSimplex(dm, &simplex));
660 PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.Nc < 0 ? dim : user.Nc, simplex, NULL, user.qorder, &fe));
661 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
662 PetscCall(PetscFEDestroy(&fe));
663 PetscCall(DMCreateDS(dm));
664
665 PetscCall(CheckFunctions(dm, user.porder, &user));
666 PetscCall(CheckTransfer(dm, INTERPOLATION, user.porder, &user));
667 PetscCall(CheckTransfer(dm, INJECTION, user.porder, &user));
668 PetscCall(DMDestroy(&dm));
669 PetscCall(PetscFinalize());
670 return 0;
671 }
672
673 /*TEST
674
675 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34)
676 # 2D/3D P_1 on a simplex
677 test:
678 suffix: p1
679 requires: triangle ctetgen
680 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}
681 test:
682 suffix: p1_pragmatic
683 requires: triangle ctetgen pragmatic
684 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}
685 test:
686 suffix: p1_adapt
687 requires: triangle ctetgen
688 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
690 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34)
691 # 2D/3D P_2 on a simplex
692 test:
693 suffix: p2
694 requires: triangle ctetgen
695 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}
696 test:
697 suffix: p2_pragmatic
698 requires: triangle ctetgen pragmatic
699 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
701 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34)
702 # TODO This is broken. Check ex3 which worked
703 # 2D/3D P_3 on a simplex
704 test:
705 TODO: gll Lagrange nodes break this
706 suffix: p3
707 requires: triangle ctetgen !single
708 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}
709 test:
710 TODO: gll Lagrange nodes break this
711 suffix: p3_pragmatic
712 requires: triangle ctetgen pragmatic !single
713 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
715 # 2D/3D Q_1 on a tensor cell
716 test:
717 suffix: q1
718 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
720 # 2D/3D Q_2 on a tensor cell
721 test:
722 suffix: q2
723 requires: !single
724 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
726 # 2D/3D Q_3 on a tensor cell
727 test:
728 TODO: gll Lagrange nodes break this
729 suffix: q3
730 requires: !single
731 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
733 # 2D/3D P_1disc on a triangle/quadrilateral
734 # TODO Missing injection functional for simplices
735 test:
736 suffix: p1d
737 requires: triangle ctetgen
738 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
740 TEST*/
741