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