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