1 #include "../include/petscutils.h" 2 3 // ----------------------------------------------------------------------------- 4 // Convert PETSc MemType to libCEED MemType 5 // ----------------------------------------------------------------------------- 6 CeedMemType MemTypeP2C(PetscMemType mem_type) { 7 return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 8 } 9 10 // ----------------------------------------------------------------------------- 11 // Apply 3D Kershaw mesh transformation 12 // ----------------------------------------------------------------------------- 13 // Transition from a value of "a" for x=0, to a value of "b" for x=1. Optionally 14 // smooth -- see the commented versions at the end. 15 static double step(const double a, const double b, double x) { 16 if (x <= 0) return a; 17 if (x >= 1) return b; 18 return a + (b-a) * (x); 19 } 20 21 // 1D transformation at the right boundary 22 static double right(const double eps, const double x) { 23 return (x <= 0.5) ? (2-eps) * x : 1 + eps*(x-1); 24 } 25 26 // 1D transformation at the left boundary 27 static double left(const double eps, const double x) { 28 return 1-right(eps,1-x); 29 } 30 31 // Apply 3D Kershaw mesh transformation 32 // The eps parameters are in (0, 1] 33 // Uniform mesh is recovered for eps=1 34 PetscErrorCode Kershaw(DM dm_orig, PetscScalar eps) { 35 PetscErrorCode ierr; 36 Vec coord; 37 PetscInt ncoord; 38 PetscScalar *c; 39 40 PetscFunctionBeginUser; 41 ierr = DMGetCoordinatesLocal(dm_orig, &coord); CHKERRQ(ierr); 42 ierr = VecGetLocalSize(coord, &ncoord); CHKERRQ(ierr); 43 ierr = VecGetArray(coord, &c); CHKERRQ(ierr); 44 45 for (PetscInt i = 0; i < ncoord; i += 3) { 46 PetscScalar x = c[i], y = c[i+1], z = c[i+2]; 47 PetscInt layer = x*6; 48 PetscScalar lambda = (x-layer/6.0)*6; 49 c[i] = x; 50 51 switch (layer) { 52 case 0: 53 c[i+1] = left(eps, y); 54 c[i+2] = left(eps, z); 55 break; 56 case 1: 57 case 4: 58 c[i+1] = step(left(eps, y), right(eps, y), lambda); 59 c[i+2] = step(left(eps, z), right(eps, z), lambda); 60 break; 61 case 2: 62 c[i+1] = step(right(eps, y), left(eps, y), lambda/2); 63 c[i+2] = step(right(eps, z), left(eps, z), lambda/2); 64 break; 65 case 3: 66 c[i+1] = step(right(eps, y), left(eps, y), (1+lambda)/2); 67 c[i+2] = step(right(eps, z), left(eps, z), (1+lambda)/2); 68 break; 69 default: 70 c[i+1] = right(eps, y); 71 c[i+2] = right(eps, z); 72 } 73 } 74 ierr = VecRestoreArray(coord, &c); CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 // ----------------------------------------------------------------------------- 79 // Create BC label 80 // ----------------------------------------------------------------------------- 81 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) { 82 83 DMLabel label; 84 85 PetscFunctionBeginUser; 86 87 PetscCall(DMCreateLabel(dm, name)); 88 PetscCall(DMGetLabel(dm, name, &label)); 89 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, label)); 90 PetscCall(DMPlexLabelComplete(dm, label)); 91 92 PetscFunctionReturn(0); 93 }; 94 95 // ----------------------------------------------------------------------------- 96 // This function sets up a DM for a given degree 97 // ----------------------------------------------------------------------------- 98 PetscErrorCode SetupDMByDegree(DM dm, PetscInt p_degree, PetscInt q_extra, 99 PetscInt num_comp_u, PetscInt dim, bool enforce_bc) { 100 PetscInt ierr, marker_ids[1] = {1}; 101 PetscInt q_degree = p_degree + q_extra; 102 PetscFE fe; 103 MPI_Comm comm; 104 PetscBool is_simplex = PETSC_TRUE; 105 106 PetscFunctionBeginUser; 107 108 // Check if simplex or tensor-product mesh 109 ierr = DMPlexIsSimplex(dm, &is_simplex); CHKERRQ(ierr); 110 // Setup FE 111 ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); 112 ierr = PetscFECreateLagrange(comm, dim, num_comp_u, is_simplex, p_degree, 113 q_degree, &fe); CHKERRQ(ierr); 114 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 115 ierr = DMCreateDS(dm); CHKERRQ(ierr); 116 117 { 118 // create FE field for coordinates 119 PetscFE fe_coords; 120 PetscInt num_comp_coord; 121 ierr = DMGetCoordinateDim(dm, &num_comp_coord); CHKERRQ(ierr); 122 ierr = PetscFECreateLagrange(comm, dim, num_comp_coord, is_simplex, 1, q_degree, 123 &fe_coords); CHKERRQ(ierr); 124 ierr = DMProjectCoordinates(dm, fe_coords); CHKERRQ(ierr); 125 ierr = PetscFEDestroy(&fe_coords); CHKERRQ(ierr); 126 } 127 128 // Setup Dirichlet BC 129 // Note bp1, bp2 are projection and we don't need to apply BC 130 // For bp3,bp4, the target function is zero on the boundaries 131 // So we pass bcFunc = NULL in DMAddBoundary function 132 if (enforce_bc) { 133 PetscBool has_label; 134 DMHasLabel(dm, "marker", &has_label); 135 if (!has_label) {CreateBCLabel(dm, "marker");} 136 DMLabel label; 137 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 138 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, 139 marker_ids, 0, 0, NULL, NULL, 140 NULL, NULL, NULL); CHKERRQ(ierr); 141 PetscCall(DMSetOptionsPrefix(dm, "final_")); 142 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 143 } 144 145 if (!is_simplex) { 146 DM dm_coord; 147 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 148 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 149 CHKERRQ(ierr); 150 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 151 CHKERRQ(ierr); 152 } 153 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 154 155 PetscFunctionReturn(0); 156 }; 157 158 // ----------------------------------------------------------------------------- 159 // Get CEED restriction data from DMPlex 160 // ----------------------------------------------------------------------------- 161 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, 162 DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { 163 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; 164 PetscErrorCode ierr; 165 166 PetscFunctionBeginUser; 167 168 ierr = DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, 169 &elem_size, &num_comp, &num_dof, &elem_restr_offsets); 170 CHKERRQ(ierr); 171 172 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 173 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 174 elem_restr_offsets, elem_restr); 175 ierr = PetscFree(elem_restr_offsets); CHKERRQ(ierr); 176 177 PetscFunctionReturn(0); 178 }; 179 180 // ----------------------------------------------------------------------------- 181 // Utility function - convert from DMPolytopeType to CeedElemTopology 182 // ----------------------------------------------------------------------------- 183 CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { 184 switch (cell_type) { 185 case DM_POLYTOPE_TRIANGLE: return CEED_TOPOLOGY_TRIANGLE; 186 case DM_POLYTOPE_QUADRILATERAL: return CEED_TOPOLOGY_QUAD; 187 case DM_POLYTOPE_TETRAHEDRON: return CEED_TOPOLOGY_TET; 188 case DM_POLYTOPE_HEXAHEDRON: return CEED_TOPOLOGY_HEX; 189 default: return 0; 190 } 191 } 192 193 // ----------------------------------------------------------------------------- 194 // Convert DM field to DS field 195 // ----------------------------------------------------------------------------- 196 PetscErrorCode DMFieldToDSField(DM dm, DMLabel domain_label, PetscInt dm_field, 197 PetscInt *ds_field) { 198 PetscDS ds; 199 IS field_is; 200 const PetscInt *fields; 201 PetscInt num_fields; 202 203 PetscFunctionBeginUser; 204 205 // Translate dm_field to ds_field 206 PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds)); 207 PetscCall(ISGetIndices(field_is, &fields)); 208 PetscCall(ISGetSize(field_is, &num_fields)); 209 for (PetscInt i = 0; i < num_fields; i++) { 210 if (dm_field == fields[i]) { 211 *ds_field = i; 212 break; 213 } 214 } 215 PetscCall(ISRestoreIndices(field_is, &fields)); 216 217 if (*ds_field == -1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, 218 "Could not find dm_field %" PetscInt_FMT " in DS", dm_field); 219 220 PetscFunctionReturn(0); 221 } 222 223 // ----------------------------------------------------------------------------- 224 // Create libCEED Basis from PetscTabulation 225 // ----------------------------------------------------------------------------- 226 PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, 227 PetscInt label_value, PetscInt height, PetscInt face, 228 PetscFE fe, PetscTabulation basis_tabulation, PetscQuadrature quadrature, 229 CeedBasis *basis) { 230 231 PetscInt first_point; 232 PetscInt ids[1] = {label_value}; 233 DMLabel depth_label; 234 DMPolytopeType cell_type; 235 CeedElemTopology elem_topo; 236 PetscScalar *q_points, *interp, *grad; 237 const PetscScalar *q_weights; 238 PetscDualSpace dual_space; 239 PetscInt num_dual_basis_vectors; 240 PetscInt dim, num_comp, P, Q; 241 242 PetscFunctionBeginUser; 243 244 // General basis information 245 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 246 PetscCall(PetscFEGetNumComponents(fe, &num_comp)); 247 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 248 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); 249 P = num_dual_basis_vectors / num_comp; 250 251 // Use depth label if no domain label present 252 if (!domain_label) { 253 PetscInt depth; 254 255 PetscCall(DMPlexGetDepth(dm, &depth)); 256 PetscCall(DMPlexGetDepthLabel(dm, &depth_label)); 257 ids[0] = depth - height; 258 } 259 260 // Get cell interp, grad, and quadrature data 261 PetscCall(DMGetFirstLabeledPoint(dm, dm, 262 domain_label ? domain_label : depth_label, 1, ids, height, &first_point, NULL)); 263 PetscCall(DMPlexGetCellType(dm, first_point, &cell_type)); 264 elem_topo = ElemTopologyP2C(cell_type); 265 if (!elem_topo) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, 266 "DMPlex topology not supported"); 267 { 268 size_t q_points_size; 269 const PetscScalar *q_points_petsc; 270 PetscInt q_dim; 271 272 PetscCall(PetscQuadratureGetData(quadrature, &q_dim, NULL, &Q, &q_points_petsc, 273 &q_weights)); 274 q_points_size = Q * dim * sizeof(CeedScalar); 275 PetscCall(PetscCalloc(q_points_size, &q_points)); 276 for (PetscInt q = 0; q < Q; q++) { 277 for (PetscInt d = 0; d < q_dim; 278 d++) q_points[q * dim + d] = q_points_petsc[q * q_dim + d]; 279 } 280 } 281 282 // Convert to libCEED orientation 283 { 284 PetscBool is_simplex = PETSC_FALSE; 285 IS permutation = NULL; 286 const PetscInt *permutation_indices; 287 288 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 289 if (!is_simplex) { 290 PetscSection section; 291 292 // -- Get permutation 293 PetscCall(DMGetLocalSection(dm, §ion)); 294 PetscCall(PetscSectionGetClosurePermutation(section, (PetscObject)dm, dim, 295 num_comp * P, &permutation)); 296 PetscCall(ISGetIndices(permutation, &permutation_indices)); 297 } 298 299 // -- Copy interp, grad matrices 300 PetscCall(PetscCalloc(P * Q * sizeof(CeedScalar), &interp)); 301 PetscCall(PetscCalloc(P * Q * dim * sizeof(CeedScalar), &grad)); 302 const CeedInt c = 0; 303 for (CeedInt q = 0; q < Q; q++) { 304 for (CeedInt p_ceed = 0; p_ceed < P; p_ceed++) { 305 CeedInt p_petsc = is_simplex ? (p_ceed * num_comp) : permutation_indices[p_ceed 306 * num_comp]; 307 308 interp[q * P + p_ceed] = basis_tabulation->T[0][((face * Q + q) * P * num_comp + 309 p_petsc) * num_comp + c]; 310 for (CeedInt d = 0; d < dim; d++) { 311 grad[(d * Q + q) * P + p_ceed] = basis_tabulation->T[1][((( 312 face * Q + q) * P * num_comp + p_petsc) * num_comp + c) * dim + d]; 313 } 314 } 315 } 316 317 // -- Cleanup 318 if (permutation) PetscCall(ISRestoreIndices(permutation, &permutation_indices)); 319 PetscCall(ISDestroy(&permutation)); 320 } 321 322 // Finally, create libCEED basis 323 CeedBasisCreateH1(ceed, elem_topo, num_comp, P, Q, interp, grad, q_points, 324 q_weights, basis); 325 PetscCall(PetscFree(q_points)); 326 PetscCall(PetscFree(interp)); 327 PetscCall(PetscFree(grad)); 328 329 PetscFunctionReturn(0); 330 } 331 332 // ----------------------------------------------------------------------------- 333 // Get CEED Basis from DMPlex 334 // ----------------------------------------------------------------------------- 335 PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, 336 CeedInt label_value, CeedInt height, 337 CeedInt dm_field, BPData bp_data, CeedBasis *basis) { 338 PetscDS ds; 339 PetscFE fe; 340 PetscQuadrature quadrature; 341 PetscBool is_simplex = PETSC_TRUE; 342 PetscInt ds_field = -1; 343 344 PetscFunctionBeginUser; 345 346 // Get element information 347 PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds)); 348 PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field)); 349 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe)); 350 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe)); 351 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 352 353 // Check if simplex or tensor-product mesh 354 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 355 356 // Build libCEED basis 357 if (is_simplex) { 358 PetscTabulation basis_tabulation; 359 PetscInt num_derivatives = 1, face = 0; 360 361 PetscCall(PetscFEGetCellTabulation(fe, num_derivatives, &basis_tabulation)); 362 PetscCall(BasisCreateFromTabulation(ceed, dm, domain_label, label_value, height, 363 face, fe, basis_tabulation, quadrature, basis)); 364 } else { 365 PetscDualSpace dual_space; 366 PetscInt num_dual_basis_vectors; 367 PetscInt dim, num_comp, P, Q; 368 369 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 370 PetscCall(PetscFEGetNumComponents(fe, &num_comp)); 371 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 372 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); 373 P = num_dual_basis_vectors / num_comp; 374 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &Q, NULL, NULL)); 375 376 CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim)); 377 CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim)); 378 379 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, 380 bp_data.q_mode, basis); 381 } 382 383 PetscFunctionReturn(0); 384 } 385 386 // ----------------------------------------------------------------------------- 387 // Utilities 388 // ----------------------------------------------------------------------------- 389 390 // Utility function, compute three factors of an integer 391 static void Split3(PetscInt size, PetscInt m[3], bool reverse) { 392 for (PetscInt d=0, size_left=size; d<3; d++) { 393 PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(size_left, 1./(3 - d))); 394 while (try * (size_left / try) != size_left) try++; 395 m[reverse ? 2-d : d] = try; 396 size_left /= try; 397 } 398 } 399 400 static int Max3(const PetscInt a[3]) { 401 return PetscMax(a[0], PetscMax(a[1], a[2])); 402 } 403 404 static int Min3(const PetscInt a[3]) { 405 return PetscMin(a[0], PetscMin(a[1], a[2])); 406 } 407 408 // ----------------------------------------------------------------------------- 409 // Create distribute dm 410 // ----------------------------------------------------------------------------- 411 PetscErrorCode CreateDistributedDM(RunParams rp, DM *dm) { 412 PetscErrorCode ierr; 413 414 PetscFunctionBeginUser; 415 // Setup DM 416 if (rp->read_mesh) { 417 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, rp->filename, NULL, PETSC_TRUE, 418 dm); 419 CHKERRQ(ierr); 420 } else { 421 if (rp->user_l_nodes) { 422 // Find a nicely composite number of elements no less than global nodes 423 PetscMPIInt size; 424 ierr = MPI_Comm_size(rp->comm, &size); CHKERRQ(ierr); 425 for (PetscInt g_elem = 426 PetscMax(1, size * rp->local_nodes / PetscPowInt(rp->degree, rp->dim)); 427 ; 428 g_elem++) { 429 Split3(g_elem, rp->mesh_elem, true); 430 if (Max3(rp->mesh_elem) / Min3(rp->mesh_elem) <= 2) break; 431 } 432 } 433 434 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, rp->dim, rp->simplex, 435 rp->mesh_elem, 436 NULL, NULL, NULL, PETSC_TRUE, dm); CHKERRQ(ierr); 437 } 438 439 ierr = DMSetFromOptions(*dm); CHKERRQ(ierr); 440 ierr = DMViewFromOptions(*dm, NULL, "-dm_view"); CHKERRQ(ierr); 441 442 PetscFunctionReturn(0); 443 } 444 445 // ----------------------------------------------------------------------------- 446