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 int ierr; 83 DMLabel label; 84 85 PetscFunctionBeginUser; 86 87 ierr = DMCreateLabel(dm, name); CHKERRQ(ierr); 88 ierr = DMGetLabel(dm, name, &label); CHKERRQ(ierr); 89 ierr = DMPlexMarkBoundaryFaces(dm, 1, label); CHKERRQ(ierr); 90 91 PetscFunctionReturn(0); 92 }; 93 94 // ----------------------------------------------------------------------------- 95 // This function sets up a DM for a given degree 96 // ----------------------------------------------------------------------------- 97 PetscErrorCode SetupDMByDegree(DM dm, PetscInt p_degree, PetscInt q_extra, 98 PetscInt num_comp_u, 99 PetscInt dim, bool enforce_bc, BCFunction bc_func) { 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 DM 129 //ierr = DMCreateDS(dm); CHKERRQ(ierr); 130 if (enforce_bc) { 131 PetscBool has_label; 132 DMHasLabel(dm, "marker", &has_label); 133 if (!has_label) {CreateBCLabel(dm, "marker");} 134 DMLabel label; 135 ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr); 136 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, 137 marker_ids, 0, 0, NULL, (void(*)(void))bc_func, 138 NULL, NULL, NULL); CHKERRQ(ierr); 139 } 140 141 if (!is_simplex) { 142 DM dm_coord; 143 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 144 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 145 CHKERRQ(ierr); 146 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 147 CHKERRQ(ierr); 148 } 149 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 150 151 PetscFunctionReturn(0); 152 }; 153 154 // ----------------------------------------------------------------------------- 155 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 156 // ----------------------------------------------------------------------------- 157 PetscInt Involute(PetscInt i) { 158 return i >= 0 ? i : -(i + 1); 159 }; 160 161 // ----------------------------------------------------------------------------- 162 // Get CEED restriction data from DMPlex 163 // ----------------------------------------------------------------------------- 164 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, 165 DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { 166 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; 167 PetscErrorCode ierr; 168 169 PetscFunctionBeginUser; 170 171 ierr = DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, 172 &elem_size, &num_comp, &num_dof, &elem_restr_offsets); 173 CHKERRQ(ierr); 174 175 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 176 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 177 elem_restr_offsets, elem_restr); 178 ierr = PetscFree(elem_restr_offsets); CHKERRQ(ierr); 179 180 PetscFunctionReturn(0); 181 }; 182 183 // ----------------------------------------------------------------------------- 184 // Utility function - convert from DMPolytopeType to CeedElemTopology 185 // ----------------------------------------------------------------------------- 186 CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { 187 switch (cell_type) { 188 case DM_POLYTOPE_TRIANGLE: return CEED_TOPOLOGY_TRIANGLE; 189 case DM_POLYTOPE_QUADRILATERAL: return CEED_TOPOLOGY_QUAD; 190 case DM_POLYTOPE_TETRAHEDRON: return CEED_TOPOLOGY_TET; 191 case DM_POLYTOPE_HEXAHEDRON: return CEED_TOPOLOGY_HEX; 192 default: return 0; 193 } 194 } 195 196 // ----------------------------------------------------------------------------- 197 // Get CEED Basis from DMPlex 198 // ----------------------------------------------------------------------------- 199 PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, 200 CeedInt label_value, CeedInt height, 201 CeedInt dm_field, CeedBasis *basis) { 202 PetscErrorCode ierr; 203 PetscDS ds; 204 PetscFE fe; 205 PetscQuadrature quadrature; 206 PetscBool is_simplex = PETSC_TRUE; 207 PetscInt dim, ds_field = -1, num_comp, P, Q; 208 209 PetscFunctionBeginUser; 210 211 // Get basis information 212 { 213 IS field_is; 214 const PetscInt *fields; 215 PetscInt num_fields; 216 217 ierr = DMGetRegionDS(dm, domain_label, &field_is, &ds); CHKERRQ(ierr); 218 // Translate dm_field to ds_field 219 ierr = ISGetIndices(field_is, &fields); CHKERRQ(ierr); 220 ierr = ISGetSize(field_is, &num_fields); CHKERRQ(ierr); 221 for (PetscInt i = 0; i < num_fields; i++) { 222 if (dm_field == fields[i]) { 223 ds_field = i; 224 break; 225 } 226 } 227 ierr = ISRestoreIndices(field_is, &fields); CHKERRQ(ierr); 228 } 229 if (ds_field == -1) { 230 // LCOV_EXCL_START 231 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, 232 "Could not find dm_field %D in DS", dm_field); 233 // LCOV_EXCL_STOP 234 } 235 236 // Get element information 237 { 238 PetscDualSpace dual_space; 239 PetscInt num_dual_basis_vectors; 240 241 ierr = PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe); 242 CHKERRQ(ierr); 243 ierr = PetscFEGetHeightSubspace(fe, height, &fe); CHKERRQ(ierr); 244 ierr = PetscFEGetSpatialDimension(fe, &dim); CHKERRQ(ierr); 245 ierr = PetscFEGetNumComponents(fe, &num_comp); CHKERRQ(ierr); 246 ierr = PetscFEGetDualSpace(fe, &dual_space); CHKERRQ(ierr); 247 ierr = PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors); 248 CHKERRQ(ierr); 249 P = num_dual_basis_vectors / num_comp; 250 ierr = PetscFEGetQuadrature(fe, &quadrature); CHKERRQ(ierr); 251 ierr = PetscQuadratureGetData(quadrature, NULL, NULL, &Q, NULL, NULL); 252 CHKERRQ(ierr); 253 } 254 255 // Check if simplex or tensor-product mesh 256 ierr = DMPlexIsSimplex(dm, &is_simplex); CHKERRQ(ierr); 257 // Build libCEED basis 258 if (is_simplex) { 259 PetscInt num_derivatives = 1, first_point; 260 PetscInt ids[1] = {label_value}; 261 PetscTabulation basis_tabulation; 262 const PetscScalar *q_points, *q_weights; 263 DMLabel depth_label; 264 DMPolytopeType cell_type; 265 CeedElemTopology elem_topo; 266 PetscScalar *interp, *grad; 267 268 // Use depth label if no domain label present 269 if (!domain_label) { 270 PetscInt depth; 271 272 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 273 ierr = DMPlexGetDepthLabel(dm, &depth_label); CHKERRQ(ierr); 274 ids[0] = depth - height; 275 } 276 // Get cell interp, grad, and quadrature data 277 ierr = PetscFEGetCellTabulation(fe, num_derivatives, &basis_tabulation); 278 CHKERRQ(ierr); 279 ierr = PetscQuadratureGetData(quadrature, NULL, NULL, NULL, &q_points, 280 &q_weights); CHKERRQ(ierr); 281 ierr = DMGetFirstLabeledPoint(dm, dm, domain_label ? domain_label : depth_label, 282 1, ids, height, &first_point, NULL); 283 CHKERRQ(ierr); 284 ierr = DMPlexGetCellType(dm, first_point, &cell_type); CHKERRQ(ierr); 285 elem_topo = ElemTopologyP2C(cell_type); 286 if (!elem_topo) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, 287 "DMPlex topology not supported"); 288 // Convert to libCEED orientation 289 ierr = PetscCalloc(P * Q * sizeof(PetscScalar), &interp); CHKERRQ(ierr); 290 ierr = PetscCalloc(P * Q * dim * sizeof(PetscScalar), &grad); CHKERRQ(ierr); 291 const CeedInt c = 0; 292 for (CeedInt q = 0; q < Q; q++) { 293 for (CeedInt p = 0; p < P; p++) { 294 interp[q*P + p] = basis_tabulation->T[0][(q*P + p)*num_comp*num_comp + c]; 295 for (CeedInt d = 0; d < dim; d++) { 296 grad[(d*Q + q)*P + p] = basis_tabulation->T[1][((q*P + p)*num_comp*num_comp + c) 297 *dim + d]; 298 } 299 } 300 } 301 // Finaly, create libCEED basis 302 ierr = CeedBasisCreateH1(ceed, elem_topo, num_comp, P, Q, interp, grad, 303 q_points, q_weights, basis); 304 CHKERRQ(ierr); 305 ierr = PetscFree(interp); CHKERRQ(ierr); 306 ierr = PetscFree(grad); CHKERRQ(ierr); 307 } else { 308 CeedInt P_1d = (CeedInt) round(pow(P, 1.0 / dim)); 309 CeedInt Q_1d = (CeedInt) round(pow(Q, 1.0 / dim)); 310 311 ierr = CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, 312 CEED_GAUSS, basis); 313 CHKERRQ(ierr); 314 } 315 316 PetscFunctionReturn(0); 317 }; 318 319 // ----------------------------------------------------------------------------- 320 // Utilities 321 // ----------------------------------------------------------------------------- 322 323 // Utility function, compute three factors of an integer 324 static void Split3(PetscInt size, PetscInt m[3], bool reverse) { 325 for (PetscInt d=0, size_left=size; d<3; d++) { 326 PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(size_left, 1./(3 - d))); 327 while (try * (size_left / try) != size_left) try++; 328 m[reverse ? 2-d : d] = try; 329 size_left /= try; 330 } 331 } 332 333 static int Max3(const PetscInt a[3]) { 334 return PetscMax(a[0], PetscMax(a[1], a[2])); 335 } 336 337 static int Min3(const PetscInt a[3]) { 338 return PetscMin(a[0], PetscMin(a[1], a[2])); 339 } 340 341 // ----------------------------------------------------------------------------- 342 // Create distribute dm 343 // ----------------------------------------------------------------------------- 344 PetscErrorCode CreateDistributedDM(RunParams rp, DM *dm) { 345 PetscErrorCode ierr; 346 347 PetscFunctionBeginUser; 348 // Setup DM 349 if (rp->read_mesh) { 350 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, rp->filename, NULL, PETSC_TRUE, 351 dm); 352 CHKERRQ(ierr); 353 } else { 354 if (rp->user_l_nodes) { 355 // Find a nicely composite number of elements no less than global nodes 356 PetscMPIInt size; 357 ierr = MPI_Comm_size(rp->comm, &size); CHKERRQ(ierr); 358 for (PetscInt g_elem = 359 PetscMax(1, size * rp->local_nodes / PetscPowInt(rp->degree, rp->dim)); 360 ; 361 g_elem++) { 362 Split3(g_elem, rp->mesh_elem, true); 363 if (Max3(rp->mesh_elem) / Min3(rp->mesh_elem) <= 2) break; 364 } 365 } 366 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, rp->dim, rp->simplex, 367 rp->mesh_elem, 368 NULL, NULL, NULL, PETSC_TRUE, dm); CHKERRQ(ierr); 369 } 370 371 ierr = DMSetFromOptions(*dm); CHKERRQ(ierr); 372 ierr = DMViewFromOptions(*dm, NULL, "-dm_view"); CHKERRQ(ierr); 373 374 PetscFunctionReturn(0); 375 } 376 377 // ----------------------------------------------------------------------------- 378