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 PetscErrorCode ierr; 339 340 PetscFunctionBeginUser; 341 options->qorder = 0; 342 options->Nc = PETSC_DEFAULT; 343 options->porder = 0; 344 options->m = 1; 345 options->dir = 0; 346 options->K = 0; 347 options->usePoly = PETSC_TRUE; 348 349 ierr = PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");PetscCall(ierr); 350 PetscCall(PetscOptionsInt("-qorder", "The quadrature order", "ex8.c", options->qorder, &options->qorder, NULL)); 351 PetscCall(PetscOptionsInt("-num_comp", "The number of field components", "ex8.c", options->Nc, &options->Nc, NULL)); 352 PetscCall(PetscOptionsInt("-porder", "The order of polynomials to test", "ex8.c", options->porder, &options->porder, NULL)); 353 PetscCall(PetscOptionsInt("-K", "The number of coarse modes used in optimization", "ex8.c", options->K, &options->K, NULL)); 354 PetscCall(PetscOptionsBool("-use_poly", "Use polynomials (or harmonics) to adapt interpolator", "ex8.c", options->usePoly, &options->usePoly, NULL)); 355 ierr = PetscOptionsEnd();PetscCall(ierr); 356 PetscFunctionReturn(0); 357 } 358 359 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 360 { 361 PetscFunctionBeginUser; 362 PetscCall(DMCreate(comm, dm)); 363 PetscCall(DMSetType(*dm, DMPLEX)); 364 PetscCall(DMSetFromOptions(*dm)); 365 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 366 PetscFunctionReturn(0); 367 } 368 369 /* Setup functions to approximate */ 370 static PetscErrorCode SetupFunctions(DM dm, PetscBool usePoly, PetscInt order, PetscInt dir, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), 371 PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), AppCtx *user) 372 { 373 PetscInt dim; 374 375 PetscFunctionBeginUser; 376 user->dir = dir; 377 if (usePoly) { 378 switch (order) { 379 case 0: 380 exactFuncs[0] = constant; 381 exactFuncDers[0] = constantDer; 382 break; 383 case 1: 384 exactFuncs[0] = linear; 385 exactFuncDers[0] = linearDer; 386 break; 387 case 2: 388 exactFuncs[0] = quadratic; 389 exactFuncDers[0] = quadraticDer; 390 break; 391 case 3: 392 exactFuncs[0] = cubic; 393 exactFuncDers[0] = cubicDer; 394 break; 395 case 4: 396 exactFuncs[0] = quartic; 397 exactFuncDers[0] = quarticDer; 398 break; 399 default: 400 PetscCall(DMGetDimension(dm, &dim)); 401 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %d order %d", dim, order); 402 } 403 } else { 404 user->m = order; 405 exactFuncs[0] = trig; 406 exactFuncDers[0] = trigDer; 407 } 408 PetscFunctionReturn(0); 409 } 410 411 static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), 412 PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), 413 void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user) 414 { 415 Vec u; 416 PetscReal n[3] = {1.0, 1.0, 1.0}; 417 418 PetscFunctionBeginUser; 419 PetscCall(DMGetGlobalVector(dm, &u)); 420 /* Project function into FE function space */ 421 PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u)); 422 PetscCall(VecViewFromOptions(u, NULL, "-projection_view")); 423 /* Compare approximation to exact in L_2 */ 424 PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error)); 425 PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer)); 426 PetscCall(DMRestoreGlobalVector(dm, &u)); 427 PetscFunctionReturn(0); 428 } 429 430 static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user) 431 { 432 PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 433 PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx); 434 void *exactCtxs[3]; 435 MPI_Comm comm; 436 PetscReal error, errorDer, tol = PETSC_SMALL; 437 438 PetscFunctionBeginUser; 439 exactCtxs[0] = user; 440 exactCtxs[1] = user; 441 exactCtxs[2] = user; 442 user->constants[0] = 1.0; 443 user->constants[1] = 2.0; 444 user->constants[2] = 3.0; 445 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 446 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user)); 447 PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user)); 448 /* Report result */ 449 if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol,(double) error)); 450 else PetscCall(PetscPrintf(comm, "Function tests pass for order %D at tolerance %g\n", order, (double)tol)); 451 if (errorDer > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer)); 452 else PetscCall(PetscPrintf(comm, "Function tests pass for order %D derivatives at tolerance %g\n", order, (double)tol)); 453 PetscFunctionReturn(0); 454 } 455 456 /* Compare approximation to exact in L_2 */ 457 static PetscErrorCode CheckTransferError(DM fdm, PetscBool usePoly, PetscInt order, PetscInt dir, const char *testname, Vec fu, AppCtx *user) 458 { 459 PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 460 PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx); 461 PetscReal n[3] = {1.0, 1.0, 1.0}; 462 void *exactCtxs[3]; 463 MPI_Comm comm; 464 PetscReal error, errorDer, tol = PETSC_SMALL; 465 466 PetscFunctionBeginUser; 467 exactCtxs[0] = user; 468 exactCtxs[1] = user; 469 exactCtxs[2] = user; 470 user->constants[0] = 1.0; 471 user->constants[1] = 2.0; 472 user->constants[2] = 3.0; 473 PetscCall(PetscObjectGetComm((PetscObject) fdm, &comm)); 474 PetscCall(SetupFunctions(fdm, usePoly, order, dir, exactFuncs, exactFuncDers, user)); 475 PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error)); 476 PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer)); 477 /* Report result */ 478 if (error > tol) PetscCall(PetscPrintf(comm, "%s tests FAIL for order %D at tolerance %g error %g\n", testname, order, (double)tol, (double)error)); 479 else PetscCall(PetscPrintf(comm, "%s tests pass for order %D at tolerance %g\n", testname, order, (double)tol)); 480 if (errorDer > tol) PetscCall(PetscPrintf(comm, "%s tests FAIL for order %D derivatives at tolerance %g error %g\n", testname, order, (double)tol, (double)errorDer)); 481 else PetscCall(PetscPrintf(comm, "%s tests pass for order %D derivatives at tolerance %g\n", testname, order, (double)tol)); 482 PetscFunctionReturn(0); 483 } 484 485 static PetscErrorCode CheckTransfer(DM dm, InterpType inType, PetscInt order, AppCtx *user) 486 { 487 PetscErrorCode (*exactFuncs[1]) (PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 488 PetscErrorCode (*exactFuncDers[1]) (PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx); 489 void *exactCtxs[3]; 490 DM rdm = NULL, idm = NULL, fdm = NULL; 491 Mat Interp, InterpAdapt = NULL; 492 Vec iu, fu, scaling = NULL; 493 MPI_Comm comm; 494 const char *testname = "Unknown"; 495 char checkname[PETSC_MAX_PATH_LEN]; 496 497 PetscFunctionBeginUser; 498 exactCtxs[0] = exactCtxs[1] = exactCtxs[2] = user; 499 PetscCall(PetscObjectGetComm((PetscObject)dm,&comm)); 500 PetscCall(DMRefine(dm, comm, &rdm)); 501 PetscCall(DMViewFromOptions(rdm, NULL, "-ref_dm_view")); 502 PetscCall(DMSetCoarseDM(rdm, dm)); 503 PetscCall(DMCopyDisc(dm, rdm)); 504 switch (inType) { 505 case INTERPOLATION: 506 testname = "Interpolation"; 507 idm = dm; 508 fdm = rdm; 509 break; 510 case RESTRICTION: 511 testname = "Restriction"; 512 idm = rdm; 513 fdm = dm; 514 break; 515 case INJECTION: 516 testname = "Injection"; 517 idm = rdm; 518 fdm = dm; 519 break; 520 } 521 PetscCall(DMGetGlobalVector(idm, &iu)); 522 PetscCall(DMGetGlobalVector(fdm, &fu)); 523 PetscCall(DMSetApplicationContext(dm, user)); 524 PetscCall(DMSetApplicationContext(rdm, user)); 525 /* Project function into initial FE function space */ 526 PetscCall(SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user)); 527 PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu)); 528 /* Interpolate function into final FE function space */ 529 switch (inType) { 530 case INTERPOLATION: 531 PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling)); 532 PetscCall(MatInterpolate(Interp, iu, fu)); 533 break; 534 case RESTRICTION: 535 PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling)); 536 PetscCall(MatRestrict(Interp, iu, fu)); 537 PetscCall(VecPointwiseMult(fu, scaling, fu)); 538 break; 539 case INJECTION: 540 PetscCall(DMCreateInjection(dm, rdm, &Interp)); 541 PetscCall(MatRestrict(Interp, iu, fu)); 542 break; 543 } 544 PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, testname, fu, user)); 545 if (user->K && (inType == INTERPOLATION)) { 546 KSP smoother; 547 Mat A; 548 Vec *iV, *fV; 549 PetscInt k, dim, d; 550 551 PetscCall(PetscPrintf(comm, " Adapting interpolator using %s\n", user->usePoly ? "polynomials" : "harmonics")); 552 PetscCall(DMGetDimension(dm, &dim)); 553 PetscCall(PetscMalloc2(user->K*dim, &iV, user->K*dim, &fV)); 554 /* Project coarse modes into initial and final FE function space */ 555 for (k = 0; k < user->K; ++k) { 556 for (d = 0; d < dim; ++d) { 557 PetscCall(DMGetGlobalVector(idm, &iV[k*dim+d])); 558 PetscCall(DMGetGlobalVector(fdm, &fV[k*dim+d])); 559 PetscCall(SetupFunctions(idm, user->usePoly, user->usePoly ? k : k+1, d, exactFuncs, exactFuncDers, user)); 560 PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iV[k*dim+d])); 561 PetscCall(DMProjectFunction(fdm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, fV[k*dim+d])); 562 } 563 } 564 /* Adapt interpolator */ 565 PetscCall(DMCreateMatrix(rdm, &A)); 566 PetscCall(MatShift(A, 1.0)); 567 PetscCall(KSPCreate(comm, &smoother)); 568 PetscCall(KSPSetFromOptions(smoother)); 569 PetscCall(KSPSetOperators(smoother, A, A)); 570 PetscCall(DMAdaptInterpolator(dm, rdm, Interp, smoother, user->K*dim, fV, iV, &InterpAdapt, user)); 571 /* Interpolate function into final FE function space */ 572 PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s poly", testname)); 573 PetscCall(MatInterpolate(InterpAdapt, iu, fu)); 574 PetscCall(CheckTransferError(fdm, PETSC_TRUE, order, 0, checkname, fu, user)); 575 for (k = 0; k < user->K; ++k) { 576 for (d = 0; d < dim; ++d) { 577 PetscCall(PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s trig (%D, %D)", testname, k, d)); 578 PetscCall(MatInterpolate(InterpAdapt, iV[k*dim+d], fV[k*dim+d])); 579 PetscCall(CheckTransferError(fdm, PETSC_FALSE, k+1, d, checkname, fV[k*dim+d], user)); 580 } 581 } 582 /* Cleanup */ 583 PetscCall(KSPDestroy(&smoother)); 584 PetscCall(MatDestroy(&A)); 585 for (k = 0; k < user->K; ++k) { 586 for (d = 0; d < dim; ++d) { 587 PetscCall(DMRestoreGlobalVector(idm, &iV[k*dim+d])); 588 PetscCall(DMRestoreGlobalVector(fdm, &fV[k*dim+d])); 589 } 590 } 591 PetscCall(PetscFree2(iV, fV)); 592 PetscCall(MatDestroy(&InterpAdapt)); 593 } 594 PetscCall(DMRestoreGlobalVector(idm, &iu)); 595 PetscCall(DMRestoreGlobalVector(fdm, &fu)); 596 PetscCall(MatDestroy(&Interp)); 597 PetscCall(VecDestroy(&scaling)); 598 PetscCall(DMDestroy(&rdm)); 599 PetscFunctionReturn(0); 600 } 601 602 int main(int argc, char **argv) 603 { 604 DM dm; 605 PetscFE fe; 606 AppCtx user; 607 PetscInt dim; 608 PetscBool simplex; 609 610 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 611 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 612 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 613 614 PetscCall(DMGetDimension(dm, &dim)); 615 PetscCall(DMPlexIsSimplex(dm, &simplex)); 616 PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.Nc < 0 ? dim : user.Nc, simplex, NULL, user.qorder, &fe)); 617 PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe)); 618 PetscCall(PetscFEDestroy(&fe)); 619 PetscCall(DMCreateDS(dm)); 620 621 PetscCall(CheckFunctions(dm, user.porder, &user)); 622 PetscCall(CheckTransfer(dm, INTERPOLATION, user.porder, &user)); 623 PetscCall(CheckTransfer(dm, INJECTION, user.porder, &user)); 624 PetscCall(DMDestroy(&dm)); 625 PetscCall(PetscFinalize()); 626 return 0; 627 } 628 629 /*TEST 630 631 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 632 # 2D/3D P_1 on a simplex 633 test: 634 suffix: p1 635 requires: triangle ctetgen 636 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} 637 test: 638 suffix: p1_pragmatic 639 requires: triangle ctetgen pragmatic 640 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} 641 test: 642 suffix: p1_adapt 643 requires: triangle ctetgen 644 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} 645 646 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 647 # 2D/3D P_2 on a simplex 648 test: 649 suffix: p2 650 requires: triangle ctetgen 651 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} 652 test: 653 suffix: p2_pragmatic 654 requires: triangle ctetgen pragmatic 655 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} 656 657 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 658 # TODO This is broken. Check ex3 which worked 659 # 2D/3D P_3 on a simplex 660 test: 661 TODO: gll Lagrange nodes break this 662 suffix: p3 663 requires: triangle ctetgen !single 664 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} 665 test: 666 TODO: gll Lagrange nodes break this 667 suffix: p3_pragmatic 668 requires: triangle ctetgen pragmatic !single 669 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} 670 671 # 2D/3D Q_1 on a tensor cell 672 test: 673 suffix: q1 674 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} 675 676 # 2D/3D Q_2 on a tensor cell 677 test: 678 suffix: q2 679 requires: !single 680 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} 681 682 # 2D/3D Q_3 on a tensor cell 683 test: 684 TODO: gll Lagrange nodes break this 685 suffix: q3 686 requires: !single 687 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} 688 689 # 2D/3D P_1disc on a triangle/quadrilateral 690 # TODO Missing injection functional for simplices 691 test: 692 suffix: p1d 693 requires: triangle ctetgen 694 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} 695 696 TEST*/ 697