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