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