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