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 { 151 AppCtx *user = (AppCtx *)ctx; 152 PetscInt d = user->dir; 153 154 if (Nc > 1) { 155 for (d = 0; d < Nc; ++d) u[d] = user->constants[d]; 156 } else { 157 u[0] = user->constants[d]; 158 } 159 return 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; 165 166 if (Nc > 1) { 167 for (d = 0; d < Nc; ++d) u[d] = 0.0; 168 } else { 169 u[0] = user->constants[d]; 170 } 171 return 0; 172 } 173 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; 179 180 if (Nc > 1) { 181 for (d = 0; d < Nc; ++d) u[d] = coords[d]; 182 } else { 183 u[0] = coords[d]; 184 } 185 return 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; 191 192 if (Nc > 1) { 193 PetscInt e; 194 for (d = 0; d < Nc; ++d) { 195 u[d] = 0.0; 196 for (e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e]; 197 } 198 } else { 199 u[0] = n[d]; 200 } 201 return 0; 202 } 203 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; 209 210 if (Nc > 1) { 211 if (Nc > 2) { 212 u[0] = coords[0] * coords[1]; 213 u[1] = coords[1] * coords[2]; 214 u[2] = coords[2] * coords[0]; 215 } else { 216 u[0] = coords[0] * coords[0]; 217 u[1] = coords[0] * coords[1]; 218 } 219 } else { 220 u[0] = coords[d] * coords[d]; 221 } 222 return 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; 228 229 if (Nc > 1) { 230 if (Nc > 2) { 231 u[0] = coords[1] * n[0] + coords[0] * n[1]; 232 u[1] = coords[2] * n[1] + coords[1] * n[2]; 233 u[2] = coords[2] * n[0] + coords[0] * n[2]; 234 } else { 235 u[0] = 2.0 * coords[0] * n[0]; 236 u[1] = coords[1] * n[0] + coords[0] * n[1]; 237 } 238 } else { 239 u[0] = 2.0 * coords[d] * n[d]; 240 } 241 return 0; 242 } 243 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; 249 250 if (Nc > 1) { 251 if (Nc > 2) { 252 u[0] = coords[0] * coords[0] * coords[1]; 253 u[1] = coords[1] * coords[1] * coords[2]; 254 u[2] = coords[2] * coords[2] * coords[0]; 255 } else { 256 u[0] = coords[0] * coords[0] * coords[0]; 257 u[1] = coords[0] * coords[0] * coords[1]; 258 } 259 } else { 260 u[0] = coords[d] * coords[d] * coords[d]; 261 } 262 return 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; 268 269 if (Nc > 1) { 270 if (Nc > 2) { 271 u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1]; 272 u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2]; 273 u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0]; 274 } else { 275 u[0] = 3.0 * coords[0] * coords[0] * n[0]; 276 u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1]; 277 } 278 } else { 279 u[0] = 3.0 * coords[d] * coords[d] * n[d]; 280 } 281 return 0; 282 } 283 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; 289 290 if (Nc > 1) { 291 if (Nc > 2) { 292 u[0] = coords[0] * coords[0] * coords[1] * coords[1]; 293 u[1] = coords[1] * coords[1] * coords[2] * coords[2]; 294 u[2] = coords[2] * coords[2] * coords[0] * coords[0]; 295 } else { 296 u[0] = coords[0] * coords[0] * coords[0] * coords[0]; 297 u[1] = coords[0] * coords[0] * coords[1] * coords[1]; 298 } 299 } else { 300 u[0] = coords[d] * coords[d] * coords[d] * coords[d]; 301 } 302 return 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; 308 309 if (Nc > 1) { 310 if (Nc > 2) { 311 u[0] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1]; 312 u[1] = 2.0 * coords[1] * coords[2] * coords[2] * n[1] + 2.0 * coords[1] * coords[1] * coords[2] * n[2]; 313 u[2] = 2.0 * coords[2] * coords[0] * coords[0] * n[2] + 2.0 * coords[2] * coords[2] * coords[0] * n[0]; 314 } else { 315 u[0] = 4.0 * coords[0] * coords[0] * coords[0] * n[0]; 316 u[1] = 2.0 * coords[0] * coords[1] * coords[1] * n[0] + 2.0 * coords[0] * coords[0] * coords[1] * n[1]; 317 } 318 } else { 319 u[0] = 4.0 * coords[d] * coords[d] * coords[d] * n[d]; 320 } 321 return 0; 322 } 323 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; 328 329 if (Nc > 1) { 330 for (d = 0; d < Nc; ++d) u[d] = PetscTanhReal(coords[d] - 0.5); 331 } else { 332 u[0] = PetscTanhReal(coords[d] - 0.5); 333 } 334 return 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; 340 341 if (Nc > 1) { 342 for (d = 0; d < Nc; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d]; 343 } else { 344 u[0] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d]; 345 } 346 return 0; 347 } 348 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; 353 354 if (Nc > 1) { 355 for (d = 0; d < Nc; ++d) u[d] = PetscSinReal(PETSC_PI * m * coords[d]); 356 } else { 357 u[0] = PetscSinReal(PETSC_PI * m * coords[d]); 358 } 359 return 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; 365 366 if (Nc > 1) { 367 for (d = 0; d < Nc; ++d) u[d] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d]; 368 } else { 369 u[0] = PETSC_PI * m * PetscCosReal(PETSC_PI * m * coords[d]) * n[d]; 370 } 371 return 0; 372 } 373 374 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 375 { 376 PetscFunctionBeginUser; 377 options->qorder = 0; 378 options->Nc = PETSC_DEFAULT; 379 options->porder = 0; 380 options->m = 1; 381 options->dir = 0; 382 options->K = 0; 383 options->usePoly = PETSC_TRUE; 384 385 PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex"); 386 PetscCall(PetscOptionsInt("-qorder", "The quadrature order", "ex8.c", options->qorder, &options->qorder, NULL)); 387 PetscCall(PetscOptionsInt("-num_comp", "The number of field components", "ex8.c", options->Nc, &options->Nc, NULL)); 388 PetscCall(PetscOptionsInt("-porder", "The order of polynomials to test", "ex8.c", options->porder, &options->porder, NULL)); 389 PetscCall(PetscOptionsInt("-K", "The number of coarse modes used in optimization", "ex8.c", options->K, &options->K, NULL)); 390 PetscCall(PetscOptionsBool("-use_poly", "Use polynomials (or harmonics) to adapt interpolator", "ex8.c", options->usePoly, &options->usePoly, NULL)); 391 PetscOptionsEnd(); 392 PetscFunctionReturn(0); 393 } 394 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(0); 403 } 404 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; 409 410 PetscFunctionBeginUser; 411 user->dir = dir; 412 if (usePoly) { 413 switch (order) { 414 case 0: 415 exactFuncs[0] = constant; 416 exactFuncDers[0] = constantDer; 417 break; 418 case 1: 419 exactFuncs[0] = linear; 420 exactFuncDers[0] = linearDer; 421 break; 422 case 2: 423 exactFuncs[0] = quadratic; 424 exactFuncDers[0] = quadraticDer; 425 break; 426 case 3: 427 exactFuncs[0] = cubic; 428 exactFuncDers[0] = cubicDer; 429 break; 430 case 4: 431 exactFuncs[0] = quartic; 432 exactFuncDers[0] = quarticDer; 433 break; 434 default: 435 PetscCall(DMGetDimension(dm, &dim)); 436 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order); 437 } 438 } else { 439 user->m = order; 440 exactFuncs[0] = trig; 441 exactFuncDers[0] = trigDer; 442 } 443 PetscFunctionReturn(0); 444 } 445 446 static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user) 447 { 448 Vec u; 449 PetscReal n[3] = {1.0, 1.0, 1.0}; 450 451 PetscFunctionBeginUser; 452 PetscCall(DMGetGlobalVector(dm, &u)); 453 /* Project function into FE function space */ 454 PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u)); 455 PetscCall(VecViewFromOptions(u, NULL, "-projection_view")); 456 /* Compare approximation to exact in L_2 */ 457 PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error)); 458 PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer)); 459 PetscCall(DMRestoreGlobalVector(dm, &u)); 460 PetscFunctionReturn(0); 461 } 462 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; 470 471 PetscFunctionBeginUser; 472 exactCtxs[0] = user; 473 exactCtxs[1] = user; 474 exactCtxs[2] = user; 475 user->constants[0] = 1.0; 476 user->constants[1] = 2.0; 477 user->constants[2] = 3.0; 478 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 479 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user)); 480 PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user)); 481 /* Report result */ 482 if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error)); 483 else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol)); 484 if (errorDer > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer)); 485 else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol)); 486 PetscFunctionReturn(0); 487 } 488 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; 498 499 PetscFunctionBeginUser; 500 exactCtxs[0] = user; 501 exactCtxs[1] = user; 502 exactCtxs[2] = user; 503 user->constants[0] = 1.0; 504 user->constants[1] = 2.0; 505 user->constants[2] = 3.0; 506 PetscCall(PetscObjectGetComm((PetscObject)fdm, &comm)); 507 PetscCall(SetupFunctions(fdm, usePoly, order, dir, exactFuncs, exactFuncDers, user)); 508 PetscCall(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(0); 516 } 517 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]; 529 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; 583 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 } 606 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(0); 642 } 643 644 int main(int argc, char **argv) 645 { 646 DM dm; 647 PetscFE fe; 648 AppCtx user; 649 PetscInt dim; 650 PetscBool simplex; 651 652 PetscFunctionBeginUser; 653 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 654 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 655 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 656 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)); 663 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 } 671 672 /*TEST 673 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} 688 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} 699 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} 713 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} 718 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} 724 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} 731 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} 738 739 TEST*/ 740