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 hexhedra", "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] = {user, user, user}; 539 DM rdm, idm, fdm; 540 Mat Interp, InterpAdapt = NULL; 541 Vec iu, fu, scaling; 542 MPI_Comm comm; 543 const char *testname; 544 char checkname[PETSC_MAX_PATH_LEN]; 545 PetscErrorCode ierr; 546 547 PetscFunctionBeginUser; 548 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 549 ierr = DMRefine(dm, comm, &rdm);CHKERRQ(ierr); 550 ierr = DMSetCoarseDM(rdm, dm);CHKERRQ(ierr); 551 ierr = DMCopyDisc(dm, rdm);CHKERRQ(ierr); 552 switch (inType) { 553 case INTERPOLATION: 554 testname = "Interpolation"; 555 idm = dm; 556 fdm = rdm; 557 break; 558 case RESTRICTION: 559 testname = "Restriction"; 560 idm = rdm; 561 fdm = dm; 562 break; 563 case INJECTION: 564 testname = "Injection"; 565 idm = rdm; 566 fdm = dm; 567 break; 568 } 569 ierr = DMGetGlobalVector(idm, &iu);CHKERRQ(ierr); 570 ierr = DMGetGlobalVector(fdm, &fu);CHKERRQ(ierr); 571 ierr = DMSetApplicationContext(dm, user);CHKERRQ(ierr); 572 ierr = DMSetApplicationContext(rdm, user);CHKERRQ(ierr); 573 /* Project function into initial FE function space */ 574 ierr = SetupFunctions(dm, PETSC_TRUE, order, 0, exactFuncs, exactFuncDers, user);CHKERRQ(ierr); 575 ierr = DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);CHKERRQ(ierr); 576 /* Interpolate function into final FE function space */ 577 switch (inType) { 578 case INTERPOLATION: 579 ierr = DMCreateInterpolation(dm, rdm, &Interp, &scaling);CHKERRQ(ierr); 580 ierr = MatInterpolate(Interp, iu, fu);CHKERRQ(ierr); 581 break; 582 case RESTRICTION: 583 ierr = DMCreateInterpolation(dm, rdm, &Interp, &scaling);CHKERRQ(ierr); 584 ierr = MatRestrict(Interp, iu, fu);CHKERRQ(ierr); 585 ierr = VecPointwiseMult(fu, scaling, fu);CHKERRQ(ierr); 586 break; 587 case INJECTION: 588 ierr = DMCreateInjection(dm, rdm, &Interp);CHKERRQ(ierr); 589 ierr = MatRestrict(Interp, iu, fu);CHKERRQ(ierr); 590 break; 591 } 592 ierr = CheckTransferError(fdm, PETSC_TRUE, order, 0, testname, fu, user);CHKERRQ(ierr); 593 if (user->K && (inType == INTERPOLATION)) { 594 KSP smoother; 595 Mat A; 596 Vec *iV, *fV; 597 PetscInt k, dim, d; 598 599 ierr = PetscPrintf(comm, " Adapting interpolator using %s\n", user->usePoly ? "polynomials" : "harmonics");CHKERRQ(ierr); 600 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 601 ierr = PetscMalloc2(user->K*dim, &iV, user->K*dim, &fV);CHKERRQ(ierr); 602 /* Project coarse modes into initial and final FE function space */ 603 for (k = 0; k < user->K; ++k) { 604 for (d = 0; d < dim; ++d) { 605 ierr = DMGetGlobalVector(idm, &iV[k*dim+d]);CHKERRQ(ierr); 606 ierr = DMGetGlobalVector(fdm, &fV[k*dim+d]);CHKERRQ(ierr); 607 ierr = SetupFunctions(idm, user->usePoly, user->usePoly ? k : k+1, d, exactFuncs, exactFuncDers, user);CHKERRQ(ierr); 608 ierr = DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iV[k*dim+d]);CHKERRQ(ierr); 609 ierr = DMProjectFunction(fdm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, fV[k*dim+d]);CHKERRQ(ierr); 610 } 611 } 612 /* Adapt interpolator */ 613 ierr = DMCreateMatrix(rdm, &A);CHKERRQ(ierr); 614 ierr = MatShift(A, 1.0);CHKERRQ(ierr); 615 ierr = KSPCreate(comm, &smoother);CHKERRQ(ierr); 616 ierr = KSPSetFromOptions(smoother);CHKERRQ(ierr); 617 ierr = KSPSetOperators(smoother, A, A);CHKERRQ(ierr); 618 ierr = DMAdaptInterpolator(dm, rdm, Interp, smoother, user->K*dim, fV, iV, &InterpAdapt, user);CHKERRQ(ierr); 619 /* Interpolate function into final FE function space */ 620 ierr = PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s poly", testname);CHKERRQ(ierr); 621 ierr = MatInterpolate(InterpAdapt, iu, fu);CHKERRQ(ierr); 622 ierr = CheckTransferError(fdm, PETSC_TRUE, order, 0, checkname, fu, user);CHKERRQ(ierr); 623 for (k = 0; k < user->K; ++k) { 624 for (d = 0; d < dim; ++d) { 625 ierr = PetscSNPrintf(checkname, PETSC_MAX_PATH_LEN, " %s trig (%D, %D)", testname, k, d);CHKERRQ(ierr); 626 ierr = MatInterpolate(InterpAdapt, iV[k*dim+d], fV[k*dim+d]);CHKERRQ(ierr); 627 ierr = CheckTransferError(fdm, PETSC_FALSE, k+1, d, checkname, fV[k*dim+d], user);CHKERRQ(ierr); 628 } 629 } 630 /* Cleanup */ 631 ierr = KSPDestroy(&smoother);CHKERRQ(ierr); 632 ierr = MatDestroy(&A);CHKERRQ(ierr); 633 for (k = 0; k < user->K; ++k) { 634 for (d = 0; d < dim; ++d) { 635 ierr = DMRestoreGlobalVector(idm, &iV[k*dim+d]);CHKERRQ(ierr); 636 ierr = DMRestoreGlobalVector(fdm, &fV[k*dim+d]);CHKERRQ(ierr); 637 } 638 } 639 ierr = PetscFree2(iV, fV);CHKERRQ(ierr); 640 ierr = MatDestroy(&InterpAdapt);CHKERRQ(ierr); 641 } 642 ierr = DMRestoreGlobalVector(idm, &iu);CHKERRQ(ierr); 643 ierr = DMRestoreGlobalVector(fdm, &fu);CHKERRQ(ierr); 644 ierr = MatDestroy(&Interp);CHKERRQ(ierr); 645 ierr = VecDestroy(&scaling);CHKERRQ(ierr); 646 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 647 PetscFunctionReturn(0); 648 } 649 650 int main(int argc, char **argv) 651 { 652 DM dm; 653 AppCtx user; /* user-defined work context */ 654 PetscErrorCode ierr; 655 656 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 657 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 658 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 659 ierr = PetscFECreateDefault(PETSC_COMM_WORLD, user.dim, user.Nc, user.simplex, NULL, user.qorder, &user.fe);CHKERRQ(ierr); 660 ierr = DMSetField(dm, 0, NULL, (PetscObject) user.fe);CHKERRQ(ierr); 661 ierr = DMCreateDS(dm);CHKERRQ(ierr); 662 ierr = CheckFunctions(dm, user.porder, &user);CHKERRQ(ierr); 663 ierr = CheckTransfer(dm, INTERPOLATION, user.porder, &user);CHKERRQ(ierr); 664 ierr = CheckTransfer(dm, INJECTION, user.porder, &user);CHKERRQ(ierr); 665 ierr = PetscFEDestroy(&user.fe);CHKERRQ(ierr); 666 ierr = DMDestroy(&dm);CHKERRQ(ierr); 667 ierr = PetscFinalize(); 668 return ierr; 669 } 670 671 /*TEST 672 673 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 674 # 2D/3D P_1 on a simplex 675 test: 676 suffix: p1 677 requires: triangle ctetgen 678 args: -dim {{2}separate output} -petscspace_degree 1 -num_comp 1 -qorder 1 -porder {{1}separate output} 679 test: 680 suffix: p1_pragmatic 681 requires: triangle ctetgen pragmatic 682 args: -dim {{2}separate output} -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder {{1 2}separate output} 683 test: 684 suffix: p1_adapt 685 requires: triangle ctetgen 686 args: -dim {{2}separate output} -dm_refine 3 -petscspace_degree 1 -qorder 1 -porder {{1 2}separate output} 687 688 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 689 # 2D/3D P_2 on a simplex 690 test: 691 suffix: p2 692 requires: triangle ctetgen 693 args: -dim {{2}separate output} -petscspace_degree 2 -qorder 2 -porder {{1 2 3}separate output} 694 test: 695 suffix: p2_pragmatic 696 requires: triangle ctetgen pragmatic 697 args: -dim {{2}separate output} -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder {{1 2 3}separate output} 698 699 # TODO dim 3 will not work until I get composite elements in 3D (see plexrefine.c:34) 700 # TODO This is broken. Check ex3 which worked 701 # 2D/3D P_3 on a simplex 702 test: 703 suffix: p3 704 requires: triangle ctetgen !single 705 args: -dim {{2}separate output} -petscspace_degree 3 -qorder 3 -porder {{1 2 3 4}separate output} 706 test: 707 suffix: p3_pragmatic 708 requires: triangle ctetgen pragmatic !single 709 args: -dim {{2}separate output} -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder {{1 2 3 4}separate output} 710 711 # 2D/3D Q_1 on a tensor cell 712 test: 713 suffix: q1 714 requires: mpi_type_get_envelope 715 args: -dim {{2 3}separate output} -simplex 0 -petscspace_degree 1 -qorder 1 -porder {{1 2}separate output} 716 717 # 2D/3D Q_2 on a tensor cell 718 test: 719 suffix: q2 720 requires: mpi_type_get_envelope 721 args: -dim {{2 3}separate output} -simplex 0 -petscspace_degree 2 -qorder 2 -porder {{1 2 3}separate output} 722 723 # 2D/3D Q_3 on a tensor cell 724 test: 725 suffix: q3 726 requires: mpi_type_get_envelope !single 727 args: -dim {{2 3}separate output} -simplex 0 -petscspace_degree 3 -qorder 3 -porder {{1 2 3 4}separate output} 728 729 # 2D/3D P_1disc on a triangle/quadrilateral 730 # TODO Missing injection functional for simplices 731 test: 732 suffix: p1d 733 requires: triangle ctetgen 734 args: -dim {{2}separate output} -simplex {{0}separate output} -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder {{1 2}separate output} 735 736 TEST*/ 737