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