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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
372: }

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;

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: }

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: }

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;

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: }

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};

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: }

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;

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: }

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;

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(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, 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];

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