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 */ 149 PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 150 AppCtx *user = (AppCtx *)ctx; 151 PetscInt d = user->dir; 152 153 if (Nc > 1) { 154 for (d = 0; d < Nc; ++d) u[d] = user->constants[d]; 155 } else { 156 u[0] = user->constants[d]; 157 } 158 return 0; 159 } 160 PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 161 AppCtx *user = (AppCtx *)ctx; 162 PetscInt d = user->dir; 163 164 if (Nc > 1) { 165 for (d = 0; d < Nc; ++d) u[d] = 0.0; 166 } else { 167 u[0] = user->constants[d]; 168 } 169 return 0; 170 } 171 172 /* u = x */ 173 PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 174 AppCtx *user = (AppCtx *)ctx; 175 PetscInt d = user->dir; 176 177 if (Nc > 1) { 178 for (d = 0; d < Nc; ++d) u[d] = coords[d]; 179 } else { 180 u[0] = coords[d]; 181 } 182 return 0; 183 } 184 PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 185 AppCtx *user = (AppCtx *)ctx; 186 PetscInt d = user->dir; 187 188 if (Nc > 1) { 189 PetscInt e; 190 for (d = 0; d < Nc; ++d) { 191 u[d] = 0.0; 192 for (e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e]; 193 } 194 } else { 195 u[0] = n[d]; 196 } 197 return 0; 198 } 199 200 /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */ 201 PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 202 AppCtx *user = (AppCtx *)ctx; 203 PetscInt d = user->dir; 204 205 if (Nc > 1) { 206 if (Nc > 2) { 207 u[0] = coords[0] * coords[1]; 208 u[1] = coords[1] * coords[2]; 209 u[2] = coords[2] * coords[0]; 210 } else { 211 u[0] = coords[0] * coords[0]; 212 u[1] = coords[0] * coords[1]; 213 } 214 } else { 215 u[0] = coords[d] * coords[d]; 216 } 217 return 0; 218 } 219 PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 220 AppCtx *user = (AppCtx *)ctx; 221 PetscInt d = user->dir; 222 223 if (Nc > 1) { 224 if (Nc > 2) { 225 u[0] = coords[1] * n[0] + coords[0] * n[1]; 226 u[1] = coords[2] * n[1] + coords[1] * n[2]; 227 u[2] = coords[2] * n[0] + coords[0] * n[2]; 228 } else { 229 u[0] = 2.0 * coords[0] * n[0]; 230 u[1] = coords[1] * n[0] + coords[0] * n[1]; 231 } 232 } else { 233 u[0] = 2.0 * coords[d] * n[d]; 234 } 235 return 0; 236 } 237 238 /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */ 239 PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 240 AppCtx *user = (AppCtx *)ctx; 241 PetscInt d = user->dir; 242 243 if (Nc > 1) { 244 if (Nc > 2) { 245 u[0] = coords[0] * coords[0] * coords[1]; 246 u[1] = coords[1] * coords[1] * coords[2]; 247 u[2] = coords[2] * coords[2] * coords[0]; 248 } else { 249 u[0] = coords[0] * coords[0] * coords[0]; 250 u[1] = coords[0] * coords[0] * coords[1]; 251 } 252 } else { 253 u[0] = coords[d] * coords[d] * coords[d]; 254 } 255 return 0; 256 } 257 PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 258 AppCtx *user = (AppCtx *)ctx; 259 PetscInt d = user->dir; 260 261 if (Nc > 1) { 262 if (Nc > 2) { 263 u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1]; 264 u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2]; 265 u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0]; 266 } else { 267 u[0] = 3.0 * coords[0] * coords[0] * n[0]; 268 u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1]; 269 } 270 } else { 271 u[0] = 3.0 * coords[d] * coords[d] * n[d]; 272 } 273 return 0; 274 } 275 276 /* u = x^4 or u = (x^4, x^2y^2) or u = (x^2y^2, y^2z^2, z^2x^2) */ 277 PetscErrorCode quartic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 278 AppCtx *user = (AppCtx *)ctx; 279 PetscInt d = user->dir; 280 281 if (Nc > 1) { 282 if (Nc > 2) { 283 u[0] = coords[0] * coords[0] * coords[1] * coords[1]; 284 u[1] = coords[1] * coords[1] * coords[2] * coords[2]; 285 u[2] = coords[2] * coords[2] * coords[0] * coords[0]; 286 } else { 287 u[0] = coords[0] * coords[0] * coords[0] * coords[0]; 288 u[1] = coords[0] * coords[0] * coords[1] * coords[1]; 289 } 290 } else { 291 u[0] = coords[d] * coords[d] * coords[d] * coords[d]; 292 } 293 return 0; 294 } 295 PetscErrorCode quarticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 296 AppCtx *user = (AppCtx *)ctx; 297 PetscInt d = user->dir; 298 299 if (Nc > 1) { 300 if (Nc > 2) { 301 u[0] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1]; 302 u[1] = 2.0 * coords[1] * coords[2] * coords[2] * n[1] + 2.0 * coords[1] * coords[1] * coords[2] * n[2]; 303 u[2] = 2.0 * coords[2] * coords[0] * coords[0] * n[2] + 2.0 * coords[2] * coords[2] * coords[0] * n[0]; 304 } else { 305 u[0] = 4.0 * coords[0] * coords[0] * coords[0] * n[0]; 306 u[1] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1]; 307 } 308 } else { 309 u[0] = 4.0 * coords[d] * coords[d] * coords[d] * n[d]; 310 } 311 return 0; 312 } 313 314 PetscErrorCode mytanh(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 315 AppCtx *user = (AppCtx *)ctx; 316 PetscInt d = user->dir; 317 318 if (Nc > 1) { 319 for (d = 0; d < Nc; ++d) u[d] = PetscTanhReal(coords[d] - 0.5); 320 } else { 321 u[0] = PetscTanhReal(coords[d] - 0.5); 322 } 323 return 0; 324 } 325 PetscErrorCode mytanhDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 326 AppCtx *user = (AppCtx *)ctx; 327 PetscInt d = user->dir; 328 329 if (Nc > 1) { 330 for (d = 0; d < Nc; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d]; 331 } else { 332 u[0] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d]; 333 } 334 return 0; 335 } 336 337 PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) { 338 AppCtx *user = (AppCtx *)ctx; 339 PetscInt m = user->m, d = user->dir; 340 341 if (Nc > 1) { 342 for (d = 0; d < Nc; ++d) u[d] = PetscSinReal(PETSC_PI * m * coords[d]); 343 } else { 344 u[0] = PetscSinReal(PETSC_PI * m * coords[d]); 345 } 346 return 0; 347 } 348 PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nc, PetscScalar *u, void *ctx) { 349 AppCtx *user = (AppCtx *)ctx; 350 PetscInt m = user->m, d = user->dir; 351 352 if (Nc > 1) { 353 for (d = 0; d < Nc; ++d) u[d] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d]; 354 } else { 355 u[0] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d]; 356 } 357 return 0; 358 } 359 360 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 361 PetscFunctionBeginUser; 362 options->qorder = 0; 363 options->Nc = PETSC_DEFAULT; 364 options->porder = 0; 365 options->m = 1; 366 options->dir = 0; 367 options->K = 0; 368 options->usePoly = PETSC_TRUE; 369 370 PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex"); 371 PetscCall(PetscOptionsInt("-qorder", "The quadrature order", "ex8.c", options->qorder, &options->qorder, NULL)); 372 PetscCall(PetscOptionsInt("-num_comp", "The number of field components", "ex8.c", options->Nc, &options->Nc, NULL)); 373 PetscCall(PetscOptionsInt("-porder", "The order of polynomials to test", "ex8.c", options->porder, &options->porder, NULL)); 374 PetscCall(PetscOptionsInt("-K", "The number of coarse modes used in optimization", "ex8.c", options->K, &options->K, NULL)); 375 PetscCall(PetscOptionsBool("-use_poly", "Use polynomials (or harmonics) to adapt interpolator", "ex8.c", options->usePoly, &options->usePoly, NULL)); 376 PetscOptionsEnd(); 377 PetscFunctionReturn(0); 378 } 379 380 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) { 381 PetscFunctionBeginUser; 382 PetscCall(DMCreate(comm, dm)); 383 PetscCall(DMSetType(*dm, DMPLEX)); 384 PetscCall(DMSetFromOptions(*dm)); 385 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 386 PetscFunctionReturn(0); 387 } 388 389 /* Setup functions to approximate */ 390 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) { 391 PetscInt dim; 392 393 PetscFunctionBeginUser; 394 user->dir = dir; 395 if (usePoly) { 396 switch (order) { 397 case 0: 398 exactFuncs[0] = constant; 399 exactFuncDers[0] = constantDer; 400 break; 401 case 1: 402 exactFuncs[0] = linear; 403 exactFuncDers[0] = linearDer; 404 break; 405 case 2: 406 exactFuncs[0] = quadratic; 407 exactFuncDers[0] = quadraticDer; 408 break; 409 case 3: 410 exactFuncs[0] = cubic; 411 exactFuncDers[0] = cubicDer; 412 break; 413 case 4: 414 exactFuncs[0] = quartic; 415 exactFuncDers[0] = quarticDer; 416 break; 417 default: PetscCall(DMGetDimension(dm, &dim)); SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order); 418 } 419 } else { 420 user->m = order; 421 exactFuncs[0] = trig; 422 exactFuncDers[0] = trigDer; 423 } 424 PetscFunctionReturn(0); 425 } 426 427 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) { 428 Vec u; 429 PetscReal n[3] = {1.0, 1.0, 1.0}; 430 431 PetscFunctionBeginUser; 432 PetscCall(DMGetGlobalVector(dm, &u)); 433 /* Project function into FE function space */ 434 PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u)); 435 PetscCall(VecViewFromOptions(u, NULL, "-projection_view")); 436 /* Compare approximation to exact in L_2 */ 437 PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error)); 438 PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer)); 439 PetscCall(DMRestoreGlobalVector(dm, &u)); 440 PetscFunctionReturn(0); 441 } 442 443 static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user) { 444 PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 445 PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx); 446 void *exactCtxs[3]; 447 MPI_Comm comm; 448 PetscReal error, errorDer, tol = PETSC_SMALL; 449 450 PetscFunctionBeginUser; 451 exactCtxs[0] = user; 452 exactCtxs[1] = user; 453 exactCtxs[2] = user; 454 user->constants[0] = 1.0; 455 user->constants[1] = 2.0; 456 user->constants[2] = 3.0; 457 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 458 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user)); 459 PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user)); 460 /* Report result */ 461 if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error)); 462 else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol)); 463 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)); 464 else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol)); 465 PetscFunctionReturn(0); 466 } 467 468 /* Compare approximation to exact in L_2 */ 469 static PetscErrorCode CheckTransferError(DM fdm, PetscBool usePoly, PetscInt order, PetscInt dir, const char *testname, Vec fu, AppCtx *user) { 470 PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 471 PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx); 472 PetscReal n[3] = {1.0, 1.0, 1.0}; 473 void *exactCtxs[3]; 474 MPI_Comm comm; 475 PetscReal error, errorDer, tol = PETSC_SMALL; 476 477 PetscFunctionBeginUser; 478 exactCtxs[0] = user; 479 exactCtxs[1] = user; 480 exactCtxs[2] = user; 481 user->constants[0] = 1.0; 482 user->constants[1] = 2.0; 483 user->constants[2] = 3.0; 484 PetscCall(PetscObjectGetComm((PetscObject)fdm, &comm)); 485 PetscCall(SetupFunctions(fdm, usePoly, order, dir, exactFuncs, exactFuncDers, user)); 486 PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error)); 487 PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer)); 488 /* Report result */ 489 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)); 490 else PetscCall(PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " at tolerance %g\n", testname, order, (double)tol)); 491 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)); 492 else PetscCall(PetscPrintf(comm, "%s tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", testname, order, (double)tol)); 493 PetscFunctionReturn(0); 494 } 495 496 static PetscErrorCode CheckTransfer(DM dm, InterpType inType, PetscInt order, AppCtx *user) { 497 PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 498 PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx); 499 void *exactCtxs[3]; 500 DM rdm = NULL, idm = NULL, fdm = NULL; 501 Mat Interp, InterpAdapt = NULL; 502 Vec iu, fu, scaling = NULL; 503 MPI_Comm comm; 504 const char *testname = "Unknown"; 505 char checkname[PETSC_MAX_PATH_LEN]; 506 507 PetscFunctionBeginUser; 508 exactCtxs[0] = exactCtxs[1] = exactCtxs[2] = user; 509 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 510 PetscCall(DMRefine(dm, comm, &rdm)); 511 PetscCall(DMViewFromOptions(rdm, NULL, "-ref_dm_view")); 512 PetscCall(DMSetCoarseDM(rdm, dm)); 513 PetscCall(DMCopyDisc(dm, rdm)); 514 switch (inType) { 515 case INTERPOLATION: 516 testname = "Interpolation"; 517 idm = dm; 518 fdm = rdm; 519 break; 520 case RESTRICTION: 521 testname = "Restriction"; 522 idm = rdm; 523 fdm = dm; 524 break; 525 case INJECTION: 526 testname = "Injection"; 527 idm = rdm; 528 fdm = dm; 529 break; 530 } 531 PetscCall(DMGetGlobalVector(idm, &iu)); 532 PetscCall(DMGetGlobalVector(fdm, &fu)); 533 PetscCall(DMSetApplicationContext(dm, user)); 534 PetscCall(DMSetApplicationContext(rdm, user)); 535 /* Project function into initial FE function space */ 536 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user)); 537 PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu)); 538 /* Interpolate function into final FE function space */ 539 switch (inType) { 540 case INTERPOLATION: 541 PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling)); 542 PetscCall(MatInterpolate(Interp, iu, fu)); 543 break; 544 case RESTRICTION: 545 PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling)); 546 PetscCall(MatRestrict(Interp, iu, fu)); 547 PetscCall(VecPointwiseMult(fu, scaling, fu)); 548 break; 549 case INJECTION: 550 PetscCall(DMCreateInjection(dm, rdm, &Interp)); 551 PetscCall(MatRestrict(Interp, iu, fu)); 552 break; 553 } 554 PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, testname, fu, user)); 555 if (user->K && (inType == INTERPOLATION)) { 556 KSP smoother; 557 Mat A, iVM, fVM; 558 Vec iV, fV; 559 PetscInt k, dim, d, im, fm; 560 561 PetscCall(PetscPrintf(comm, " Adapting interpolator using %s\n", user->usePoly ? "polynomials" : "harmonics")); 562 PetscCall(DMGetDimension(dm, &dim)); 563 /* Project coarse modes into initial and final FE function space */ 564 PetscCall(DMGetGlobalVector(idm, &iV)); 565 PetscCall(DMGetGlobalVector(fdm, &fV)); 566 PetscCall(VecGetLocalSize(iV, &im)); 567 PetscCall(VecGetLocalSize(fV, &fm)); 568 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)dm), im, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &iVM)); 569 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)dm), fm, PETSC_DECIDE, PETSC_DECIDE, user->K * dim, NULL, &fVM)); 570 PetscCall(DMRestoreGlobalVector(idm, &iV)); 571 PetscCall(DMRestoreGlobalVector(fdm, &fV)); 572 for (k = 0; k < user->K; ++k) { 573 for (d = 0; d < dim; ++d) { 574 PetscCall(MatDenseGetColumnVecWrite(iVM, k * dim + d, &iV)); 575 PetscCall(MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV)); 576 PetscCall(SetupFunctions(idm, user->usePoly, user->usePoly ? k : k + 1, d, exactFuncs, exactFuncDers, user)); 577 PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iV)); 578 PetscCall(DMProjectFunction(fdm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, fV)); 579 PetscCall(MatDenseRestoreColumnVecWrite(iVM, k * dim + d, &iV)); 580 PetscCall(MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV)); 581 } 582 } 583 584 /* Adapt interpolator */ 585 PetscCall(DMCreateMatrix(rdm, &A)); 586 PetscCall(MatShift(A, 1.0)); 587 PetscCall(KSPCreate(comm, &smoother)); 588 PetscCall(KSPSetFromOptions(smoother)); 589 PetscCall(KSPSetOperators(smoother, A, A)); 590 PetscCall(DMAdaptInterpolator(dm, rdm, Interp, smoother, fVM, iVM, &InterpAdapt, user)); 591 /* Interpolate function into final FE function space */ 592 PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s poly", testname)); 593 PetscCall(MatInterpolate(InterpAdapt, iu, fu)); 594 PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, checkname, fu, user)); 595 for (k = 0; k < user->K; ++k) { 596 for (d = 0; d < dim; ++d) { 597 PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s trig (%" PetscInt_FMT ", %" PetscInt_FMT ")", testname, k, d)); 598 PetscCall(MatDenseGetColumnVecRead(iVM, k * dim + d, &iV)); 599 PetscCall(MatDenseGetColumnVecWrite(fVM, k * dim + d, &fV)); 600 PetscCall(MatInterpolate(InterpAdapt, iV, fV)); 601 PetscCall(CheckTransferError(fdm, PETSC_FALSE, k + 1, d, checkname, fV, user)); 602 PetscCall(MatDenseRestoreColumnVecRead(iVM, k * dim + d, &iV)); 603 PetscCall(MatDenseRestoreColumnVecWrite(fVM, k * dim + d, &fV)); 604 } 605 } 606 /* Cleanup */ 607 PetscCall(KSPDestroy(&smoother)); 608 PetscCall(MatDestroy(&A)); 609 PetscCall(MatDestroy(&InterpAdapt)); 610 PetscCall(MatDestroy(&iVM)); 611 PetscCall(MatDestroy(&fVM)); 612 } 613 PetscCall(DMRestoreGlobalVector(idm, &iu)); 614 PetscCall(DMRestoreGlobalVector(fdm, &fu)); 615 PetscCall(MatDestroy(&Interp)); 616 PetscCall(VecDestroy(&scaling)); 617 PetscCall(DMDestroy(&rdm)); 618 PetscFunctionReturn(0); 619 } 620 621 int main(int argc, char **argv) { 622 DM dm; 623 PetscFE fe; 624 AppCtx user; 625 PetscInt dim; 626 PetscBool simplex; 627 628 PetscFunctionBeginUser; 629 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 630 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 631 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 632 633 PetscCall(DMGetDimension(dm, &dim)); 634 PetscCall(DMPlexIsSimplex(dm, &simplex)); 635 PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.Nc < 0 ? dim : user.Nc, simplex, NULL, user.qorder, &fe)); 636 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 637 PetscCall(PetscFEDestroy(&fe)); 638 PetscCall(DMCreateDS(dm)); 639 640 PetscCall(CheckFunctions(dm, user.porder, &user)); 641 PetscCall(CheckTransfer(dm, INTERPOLATION, user.porder, &user)); 642 PetscCall(CheckTransfer(dm, INJECTION, user.porder, &user)); 643 PetscCall(DMDestroy(&dm)); 644 PetscCall(PetscFinalize()); 645 return 0; 646 } 647 648 /*TEST 649 650 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 651 # 2D/3D P_1 on a simplex 652 test: 653 suffix: p1 654 requires: triangle ctetgen 655 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} 656 test: 657 suffix: p1_pragmatic 658 requires: triangle ctetgen pragmatic 659 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} 660 test: 661 suffix: p1_adapt 662 requires: triangle ctetgen 663 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} 664 665 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 666 # 2D/3D P_2 on a simplex 667 test: 668 suffix: p2 669 requires: triangle ctetgen 670 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} 671 test: 672 suffix: p2_pragmatic 673 requires: triangle ctetgen pragmatic 674 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} 675 676 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 677 # TODO This is broken. Check ex3 which worked 678 # 2D/3D P_3 on a simplex 679 test: 680 TODO: gll Lagrange nodes break this 681 suffix: p3 682 requires: triangle ctetgen !single 683 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} 684 test: 685 TODO: gll Lagrange nodes break this 686 suffix: p3_pragmatic 687 requires: triangle ctetgen pragmatic !single 688 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} 689 690 # 2D/3D Q_1 on a tensor cell 691 test: 692 suffix: q1 693 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} 694 695 # 2D/3D Q_2 on a tensor cell 696 test: 697 suffix: q2 698 requires: !single 699 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} 700 701 # 2D/3D Q_3 on a tensor cell 702 test: 703 TODO: gll Lagrange nodes break this 704 suffix: q3 705 requires: !single 706 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} 707 708 # 2D/3D P_1disc on a triangle/quadrilateral 709 # TODO Missing injection functional for simplices 710 test: 711 suffix: p1d 712 requires: triangle ctetgen 713 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} 714 715 TEST*/ 716