Skip to content

Commit

Permalink
updated tests for new axisymmetric.
Browse files Browse the repository at this point in the history
  • Loading branch information
mmcgurn committed Jan 28, 2024
1 parent 5244a90 commit 9ab1a1b
Show file tree
Hide file tree
Showing 8 changed files with 126 additions and 78 deletions.
7 changes: 6 additions & 1 deletion src/domain/descriptions/axisymmetric.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@ ablate::domain::descriptions::Axisymmetric::Axisymmetric(const std::vector<Petsc

void ablate::domain::descriptions::Axisymmetric::BuildTopology(PetscInt cell, PetscInt *cellNodes) const {
// Note that the cell/vertex ordering grows in shells so that tri prism elements are specified first and together.
// flip the cell order so that hexes appear before tri prisms
cell = CellReverser(cell);

// determine the type of element
if (cell >= numberTriPrismCells) {
Expand Down Expand Up @@ -113,7 +115,10 @@ void ablate::domain::descriptions::Axisymmetric::SetCoordinate(PetscInt node, Pe
coordinate[1] += radius * radiusFactor * PetscSinReal(2.0 * nodeRotationIndex * PETSC_PI / numberWedges);
}

DMPolytopeType ablate::domain::descriptions::Axisymmetric::GetCellType(PetscInt cell) const { return cell < numberTriPrismCells ? DM_POLYTOPE_TRI_PRISM : DM_POLYTOPE_HEXAHEDRON; }
DMPolytopeType ablate::domain::descriptions::Axisymmetric::GetCellType(PetscInt cell) const {
// flip the cell order so that hexes appear before tri prisms
return CellReverser(cell) < numberTriPrismCells ? DM_POLYTOPE_TRI_PRISM : DM_POLYTOPE_HEXAHEDRON;
}

std::shared_ptr<ablate::domain::Region> ablate::domain::descriptions::Axisymmetric::GetRegion(const std::set<PetscInt> &face) const {
// check to see if each node in on the surface
Expand Down
7 changes: 7 additions & 0 deletions src/domain/descriptions/axisymmetric.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,13 @@ class Axisymmetric : public ablate::domain::descriptions::MeshDescription {
const static inline std::shared_ptr<ablate::domain::Region> lowerCapBoundary = std::make_shared<ablate::domain::Region>("lowerCap");
const static inline std::shared_ptr<ablate::domain::Region> upperCapBoundary = std::make_shared<ablate::domain::Region>("upperCap");

/**
* Reverse the cell order to make sure hexes appear before tri prism
* @param in
* @return
*/
inline PetscInt CellReverser(PetscInt in) const { return numberCells - in - 1; }

public:
/**
* generate and precompute a bunch of the required parameters
Expand Down
55 changes: 44 additions & 11 deletions src/domain/modifiers/extrudeLabel.cpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "extrudeLabel.hpp"
#include <petsc/private/dmimpl.h>
#include <petsc/private/petscfeimpl.h>
#include <petscdmplextransform.h>
#include <utility>
#include "tagLabelInterface.hpp"
Expand Down Expand Up @@ -84,18 +85,49 @@ void ablate::domain::modifiers::ExtrudeLabel::Modify(DM &dm) {
originalRegion->CreateLabel(dmAdapt, originalRegionLabel, originalRegionValue);
extrudedRegion->CreateLabel(dmAdapt, extrudedRegionLabel, extrudedRegionValue);

// Determine the current max cell int
PetscInt originalMaxCell;
DMPlexGetHeightStratum(dm, 0, nullptr, &originalMaxCell) >> utilities::PetscUtilities::checkError;

// March over each cell in this rank and determine if it is original or not
PetscInt cStart, cEnd;
DMPlexGetHeightStratum(dmAdapt, 0, &cStart, &cEnd) >> utilities::PetscUtilities::checkError;
for (PetscInt c = cStart; c < cEnd; ++c) {
if (c < originalMaxCell) {
DMLabelSetValue(originalRegionLabel, c, originalRegionValue) >> utilities::PetscUtilities::checkError;
} else {
DMLabelSetValue(extrudedRegionLabel, c, extrudedRegionValue) >> utilities::PetscUtilities::checkError;
PetscInt cAdaptStart, cAdaptEnd;
DMPlexGetHeightStratum(dmAdapt, 0, &cAdaptStart, &cAdaptEnd) >> utilities::PetscUtilities::checkError;

// cell the depths of the cell layer
PetscInt cellDepth;
DMPlexGetDepth(dm, &cellDepth) >> utilities::PetscUtilities::checkError;
DMLabel depthAdaptLabel, ctAdaptLabel, ctLabel;
DMPlexGetDepthLabel(dmAdapt, &depthAdaptLabel) >> utilities::PetscUtilities::checkError;
DMPlexGetCellTypeLabel(dmAdapt, &ctAdaptLabel) >> utilities::PetscUtilities::checkError;
DMPlexGetCellTypeLabel(dm, &ctLabel) >> utilities::PetscUtilities::checkError;

// because the new cells can be intertwined with the old cells for mixed use we need to do this cell type by cell type
for (PetscInt cellType = 0; cellType < DM_NUM_POLYTOPES; ++cellType) {
auto ict = (DMPolytopeType)cellType;

// get the new range for this cell type
PetscInt tAdaptStart, tAdaptEnd;
DMLabelGetStratumBounds(ctAdaptLabel, ict, &tAdaptStart, &tAdaptEnd) >> utilities::PetscUtilities::checkError;

// only check if there are cell of this type
if (tAdaptStart < 0) {
continue;
}
// determine the depth of this cell type
PetscInt cellTypeDepth;
DMLabelGetValue(depthAdaptLabel, tAdaptStart, &cellTypeDepth) >> utilities::PetscUtilities::checkError;
if (cellTypeDepth != cellDepth) {
continue;
}

// Get the original range for this cell type
PetscInt tStart, tEnd;
DMLabelGetStratumBounds(ctLabel, ict, &tStart, &tEnd) >> utilities::PetscUtilities::checkError;
PetscInt numberOldCells = tStart >= 0 ? tEnd - tStart : 0;

// march over each new cell
for (PetscInt c = tAdaptStart; c < tAdaptEnd; ++c) {
if ((c - tAdaptStart) < numberOldCells) {
DMLabelSetValue(originalRegionLabel, c, originalRegionValue) >> utilities::PetscUtilities::checkError;
} else {
DMLabelSetValue(extrudedRegionLabel, c, extrudedRegionValue) >> utilities::PetscUtilities::checkError;
}
}
}

Expand Down Expand Up @@ -134,6 +166,7 @@ PetscErrorCode ablate::domain::modifiers::ExtrudeLabel::DMPlexTransformAdaptLabe
PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
PetscCall(DMCopyDisc(dm, *rdm));
PetscCall(DMPlexTransformDestroy(&tr));
((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
PetscFunctionReturn(0);
}

Expand Down
6 changes: 3 additions & 3 deletions src/domain/range.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
#include <petsc/private/dmpleximpl.h> // For ISIntersect_Caching_Internal, used in ablate::domain::SubDomain::GetRange
#include "utilities/petscUtilities.hpp"

void ablate::domain::GetRange(DM dm, const std::shared_ptr<ablate::domain::Region>& region, PetscInt depth, ablate::domain::Range &range) {
void ablate::domain::GetRange(DM dm, const std::shared_ptr<ablate::domain::Region> &region, PetscInt depth, ablate::domain::Range &range) {
// Start out getting all the points
IS allPointIS;
DMGetStratumIS(dm, "dim", depth, &allPointIS) >> utilities::PetscUtilities::checkError;
Expand Down Expand Up @@ -39,14 +39,14 @@ void ablate::domain::GetRange(DM dm, const std::shared_ptr<ablate::domain::Regio
ISDestroy(&allPointIS) >> utilities::PetscUtilities::checkError;
}

void ablate::domain::GetCellRange(DM dm, const std::shared_ptr<ablate::domain::Region>& region, ablate::domain::Range &cellRange) {
void ablate::domain::GetCellRange(DM dm, const std::shared_ptr<ablate::domain::Region> &region, ablate::domain::Range &cellRange) {
// Start out getting all the cells
PetscInt depth;
DMPlexGetDepth(dm, &depth) >> utilities::PetscUtilities::checkError;
ablate::domain::GetRange(dm, region, depth, cellRange);
}

void ablate::domain::GetFaceRange(DM dm, const std::shared_ptr<ablate::domain::Region>& region, ablate::domain::Range &faceRange) {
void ablate::domain::GetFaceRange(DM dm, const std::shared_ptr<ablate::domain::Region> &region, ablate::domain::Range &faceRange) {
// Start out getting all the faces
PetscInt depth;
DMPlexGetDepth(dm, &depth) >> utilities::PetscUtilities::checkError;
Expand Down
10 changes: 5 additions & 5 deletions src/domain/range.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@ namespace ablate::domain {
*/
struct Range {
IS is = nullptr;
PetscInt start;
PetscInt end;
PetscInt start{};
PetscInt end{};
const PetscInt *points = nullptr;

/**
Expand All @@ -29,23 +29,23 @@ struct Range {
* @param depth
* @param range
*/
void GetRange(DM dm, const std::shared_ptr<Region>& region, PetscInt depth, Range &range);
void GetRange(DM dm, const std::shared_ptr<Region> &region, PetscInt depth, Range &range);

/**
* Get the range of cells defined over the region for this solver.
* @param dm
* @param region
* @param cellRange
*/
void GetCellRange(DM dm, const std::shared_ptr<Region>& region, Range &cellRange);
void GetCellRange(DM dm, const std::shared_ptr<Region> &region, Range &cellRange);

/**
* Get the range of faces/edges defined over the region for this solver.
* @param dm
* @param region
* @param faceRange
*/
void GetFaceRange(DM dm, const std::shared_ptr<Region>& region, Range &faceRange);
void GetFaceRange(DM dm, const std::shared_ptr<Region> &region, Range &faceRange);

/**
* Restores the range
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -29,43 +29,43 @@ timestepper:
length: 0.5
radius: ".1 - .1*z*z"
numberWedges: 8
numberSlices: 2
numberShells: 2
numberSlices: 4
numberShells: 4
# output the result label for interior cells
options:
dm_label_view: boundaryCells
dm_label_view: flowRegion
# force the new mesh to check everything
dm_plex_check_all: ""

# specify any modifications to be performed to the mesh/domain
modifiers:
# - # use the newly labels to extrude the boundary. Do not extrude the cell
# !ablate::domain::modifiers::ExtrudeLabel
# regions:
# - name: boundary
# # mark all the resulting boundary faces with boundaryFaces label
# boundaryRegion:
# name: boundaryFaces
# # tag the original mesh as the flow region
# originalRegion:
# name: flowRegion
# # tag the new boundary cells for easy boundary condition specifications
# extrudedRegion:
# name: boundaryCells
- # use the newly labels to extrude the boundary. Do not extrude the cell
!ablate::domain::modifiers::ExtrudeLabel
regions:
- name: boundary
# mark all the resulting boundary faces with boundaryFaces label
boundaryRegion:
name: boundaryFaces
# tag the original mesh as the flow region
originalRegion:
name: flowRegion
# tag the new boundary cells for easy boundary condition specifications
extrudedRegion:
name: boundaryCells

- !ablate::domain::modifiers::DistributeWithGhostCells
ghostCellDepth: 2

# setup some dummy fields
fields:
- name: exampleFVField
components: [ "xx", "yy", "zz" ]
components: [ "xx"]
type: FVM

# initialize the dummy field
initialization:
- fieldName: "exampleFVField"
field: "x, y, z"
field: "x"

# this is a test input file with no solvers
solvers: [ ]
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,14 @@ exampleAxisymmetricMesh in 3 dimensions:
Number of 0-cells per rank: 8241
Number of 1-cells per rank: 24440
Number of 2-cells per rank: 24200 (23380)
Number of 3-cells per rank: 8000 (7200)
Number of 3-cells per rank: 8000 (800)
Labels:
celltype: 6 strata with value/size (0 (8241), 8 (800), 7 (7200), 3 (820), 1 (24440), 4 (23380))
celltype: 6 strata with value/size (0 (8241), 7 (7200), 8 (800), 3 (820), 1 (24440), 4 (23380))
depth: 4 strata with value/size (0 (8241), 1 (24440), 2 (24200), 3 (8000))
boundary: 1 strata with value/size (1 (4762))
upperCap: 1 strata with value/size (1 (801))
lowerCap: 1 strata with value/size (1 (801))
outerShell: 1 strata with value/size (1 (3240))
ResultFiles:
domain.hdf5
domain.xmf
Expand Down
Loading

0 comments on commit 9ab1a1b

Please sign in to comment.