diff --git a/CMakeLists.txt b/CMakeLists.txt index 06536d057..de92e4a9e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.18.4) include(config/petscCompilers.cmake) # Set the project details -project(ablateLibrary VERSION 0.12.34) +project(ablateLibrary VERSION 0.12.35) # Load the Required 3rd Party Libaries pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc) diff --git a/src/domain/boxMesh.cpp b/src/domain/boxMesh.cpp index 006e2d537..071bd6230 100644 --- a/src/domain/boxMesh.cpp +++ b/src/domain/boxMesh.cpp @@ -36,7 +36,8 @@ DM ablate::domain::BoxMesh::CreateBoxDM(const std::string& name, std::vector facesPetsc(faces.begin(), faces.end()); DM dm; - DMPlexCreateBoxMesh(PETSC_COMM_WORLD, (PetscInt)dimensions, simplex ? PETSC_TRUE : PETSC_FALSE, &facesPetsc[0], &lower[0], &upper[0], &boundaryTypes[0], PETSC_TRUE, 0, PETSC_TRUE, &dm) >> + DMPlexCreateBoxMesh( + PETSC_COMM_WORLD, (PetscInt)dimensions, simplex ? PETSC_TRUE : PETSC_FALSE, &facesPetsc[0], &lower[0], &upper[0], &boundaryTypes[0], PETSC_TRUE, (PetscInt)dimensions, PETSC_TRUE, &dm) >> utilities::PetscUtilities::checkError; PetscObjectSetName((PetscObject)dm, name.c_str()) >> utilities::PetscUtilities::checkError; return dm; diff --git a/src/domain/boxMeshBoundaryCells.cpp b/src/domain/boxMeshBoundaryCells.cpp index 04fb2bdb9..a51e5f0aa 100644 --- a/src/domain/boxMeshBoundaryCells.cpp +++ b/src/domain/boxMeshBoundaryCells.cpp @@ -66,7 +66,7 @@ DM ablate::domain::BoxMeshBoundaryCells::CreateBoxDM(const std::string& name, st // Make copy with PetscInt std::vector facesPetsc(faces.begin(), faces.end()); DM dm; - DMPlexCreateBoxMesh(PETSC_COMM_WORLD, (PetscInt)dimensions, simplex ? PETSC_TRUE : PETSC_FALSE, &facesPetsc[0], &lower[0], &upper[0], nullptr, PETSC_TRUE, 0, PETSC_TRUE, &dm) >> + DMPlexCreateBoxMesh(PETSC_COMM_WORLD, (PetscInt)dimensions, simplex ? PETSC_TRUE : PETSC_FALSE, &facesPetsc[0], &lower[0], &upper[0], nullptr, PETSC_TRUE, (PetscInt)dimensions, PETSC_TRUE, &dm) >> utilities::PetscUtilities::checkError; PetscObjectSetName((PetscObject)dm, name.c_str()) >> utilities::PetscUtilities::checkError; return dm; @@ -164,4 +164,4 @@ REGISTER(ablate::domain::Domain, ablate::domain::BoxMeshBoundaryCells, OPT(std::vector, "postModifiers", "a list of domain modifiers to apply after ghost labeling"), ARG(std::vector, "faces", "the number of faces in each direction"), ARG(std::vector, "lower", "the lower bound of the mesh"), ARG(std::vector, "upper", "the upper bound of the mesh"), - OPT(ablate::parameters::Parameters, "options", "PETSc options specific to this dm. Default value allows the dm to access global options.")); \ No newline at end of file + OPT(ablate::parameters::Parameters, "options", "PETSc options specific to this dm. Default value allows the dm to access global options.")); diff --git a/src/domain/domain.cpp b/src/domain/domain.cpp index ff5dca434..a03a25c81 100644 --- a/src/domain/domain.cpp +++ b/src/domain/domain.cpp @@ -160,6 +160,12 @@ void ablate::domain::Domain::InitializeSubDomains(const std::vectorCreateSubDomainStructures(); } + // Setup the label to extract cell ranges without FVM ghosts. + // This needs to be done after CreateSubDomainStructures so that the local section has been created + for (auto& solver : solvers) { + solver->SetupCellRangeWithoutGhost(); + } + // set all values to nan to allow for a output check if (initializations) { VecSet(solGlobalField, NAN) >> utilities::PetscUtilities::checkError; diff --git a/src/domain/meshGenerator.cpp b/src/domain/meshGenerator.cpp index de521b5a9..5ff40af5f 100644 --- a/src/domain/meshGenerator.cpp +++ b/src/domain/meshGenerator.cpp @@ -29,6 +29,19 @@ void ablate::domain::MeshGenerator::ReplaceDm(DM& originalDm, DM& replaceDm) { PetscObjectSetOptions((PetscObject)replaceDm, options) >> utilities::PetscUtilities::checkError; ((DM_Plex*)(replaceDm)->data)->useHashLocation = ((DM_Plex*)originalDm->data)->useHashLocation; + // Fix the periodic localization here until it's fixed in PETSc + DM cdm; + PetscInt localizationHeight; + PetscBool sparseLocalize; + DMGetCoordinateDM(originalDm, &cdm) >> utilities::PetscUtilities::checkError; + DMPlexGetMaxProjectionHeight(cdm, &localizationHeight) >> utilities::PetscUtilities::checkError; + DMGetSparseLocalize(cdm, &sparseLocalize) >> utilities::PetscUtilities::checkError; + + DMGetCoordinateDM(replaceDm, &cdm) >> utilities::PetscUtilities::checkError; + DMPlexSetMaxProjectionHeight(cdm, localizationHeight) >> utilities::PetscUtilities::checkError; + DMSetSparseLocalize(cdm, sparseLocalize) >> utilities::PetscUtilities::checkError; + DMLocalizeCoordinates(replaceDm) >> utilities::PetscUtilities::checkError; + DMDestroy(&originalDm) >> utilities::PetscUtilities::checkError; originalDm = replaceDm; } diff --git a/src/finiteVolume/finiteVolumeSolver.cpp b/src/finiteVolume/finiteVolumeSolver.cpp index 80d7d74e4..9f8d03e0d 100644 --- a/src/finiteVolume/finiteVolumeSolver.cpp +++ b/src/finiteVolume/finiteVolumeSolver.cpp @@ -436,24 +436,6 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::Restore(PetscViewer vie PetscFunctionReturn(0); } -void ablate::finiteVolume::FiniteVolumeSolver::GetCellRangeWithoutGhost(ablate::domain::Range& faceRange) const { - // Get the point range - DMLabel solverRegionMinusGhostLabel; - PetscInt solverRegionMinusGhostValue; - domain::Region::GetLabel(solverRegionMinusGhost, GetSubDomain().GetDM(), solverRegionMinusGhostLabel, solverRegionMinusGhostValue); - - DMLabelGetStratumIS(solverRegionMinusGhostLabel, solverRegionMinusGhostValue, &faceRange.is) >> utilities::PetscUtilities::checkError; - if (faceRange.is == nullptr) { - // There are no points in this region, so skip - faceRange.start = 0; - faceRange.end = 0; - faceRange.points = nullptr; - } else { - // Get the range - ISGetPointRange(faceRange.is, &faceRange.start, &faceRange.end, &faceRange.points) >> utilities::PetscUtilities::checkError; - } -} - PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeBoundary(PetscReal time, Vec locX, Vec locX_t) { PetscFunctionBeginUser; auto dm = subDomain->GetDM(); diff --git a/src/finiteVolume/finiteVolumeSolver.hpp b/src/finiteVolume/finiteVolumeSolver.hpp index bfa46ab3d..0c49cd824 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -205,15 +205,6 @@ class FiniteVolumeSolver : public solver::CellSolver, */ PetscErrorCode Restore(PetscViewer viewer, PetscInt sequenceNumber, PetscReal time) override; - /** - * Get the cellIS and range over valid cells in this region without ghost cells (boundary or mpi) - * @param cellIS - * @param pStart - * @param pEnd - * @param points - */ - void GetCellRangeWithoutGhost(ablate::domain::Range& faceRange) const; - /** * Returns first instance of process of type specifed * @tparam T diff --git a/src/solver/solver.cpp b/src/solver/solver.cpp index 6660559a1..1d010891e 100644 --- a/src/solver/solver.cpp +++ b/src/solver/solver.cpp @@ -4,7 +4,7 @@ #include "utilities/petscUtilities.hpp" ablate::solver::Solver::Solver(std::string solverId, std::shared_ptr region, std::shared_ptr options) - : solverId(std::move(solverId)), region(std::move(region)), petscOptions(nullptr) { + : solverId(std::move(solverId)), region(std::move(region)), regionMinusGhost(std::make_shared(solverId + "_minusGhost")), petscOptions(nullptr) { // Set the options if (options) { PetscOptionsCreate(&petscOptions) >> utilities::PetscUtilities::checkError; @@ -20,6 +20,45 @@ ablate::solver::Solver::~Solver() { void ablate::solver::Solver::Register(std::shared_ptr subDomainIn) { subDomain = std::move(subDomainIn); } +void ablate::solver::Solver::SetupCellRangeWithoutGhost() { + // Get the original range + ablate::domain::Range cellRange; + GetCellRange(cellRange); + + // create a new label + auto dm = subDomain->GetDM(); + DMCreateLabel(dm, regionMinusGhost->GetName().c_str()) >> utilities::PetscUtilities::checkError; + DMLabel regionMinusGhostLabel; + PetscInt regionMinusGhostValue; + domain::Region::GetLabel(regionMinusGhost, dm, regionMinusGhostLabel, regionMinusGhostValue); + + // Get the ghost cell label + DMLabel ghostLabel; + DMGetLabel(dm, "ghost", &ghostLabel) >> utilities::PetscUtilities::checkError; + + // check if it is an exterior boundary cell ghost + PetscInt boundaryCellStart; + DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &boundaryCellStart, nullptr) >> utilities::PetscUtilities::checkError; + + // march over every cell + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + PetscInt cell = cellRange.points ? cellRange.points[c] : c; + + // check if it is boundary ghost + PetscInt isGhost = -1; + if (ghostLabel) { + DMLabelGetValue(ghostLabel, cell, &isGhost) >> utilities::PetscUtilities::checkError; + } + + PetscInt owned; + DMPlexGetPointGlobal(dm, cell, &owned, nullptr) >> utilities::PetscUtilities::checkError; + if (owned >= 0 && isGhost < 0 && (boundaryCellStart < 0 || cell < boundaryCellStart)) { + DMLabelSetValue(regionMinusGhostLabel, cell, regionMinusGhostValue); + } + } + RestoreRange(cellRange); +} + void ablate::solver::Solver::PreStage(TS ts, PetscReal stagetime) { for (auto &function : preStageFunctions) { function(ts, *this, stagetime); @@ -47,6 +86,24 @@ static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], Pe return 0; } +void ablate::solver::Solver::GetCellRangeWithoutGhost(ablate::domain::Range &cellRange) const { + // Get the point range + DMLabel regionMinusGhostLabel; + PetscInt regionMinusGhostValue; + domain::Region::GetLabel(regionMinusGhost, GetSubDomain().GetDM(), regionMinusGhostLabel, regionMinusGhostValue); + + DMLabelGetStratumIS(regionMinusGhostLabel, regionMinusGhostValue, &cellRange.is) >> utilities::PetscUtilities::checkError; + if (cellRange.is == nullptr) { + // There are no points in this region, so skip + cellRange.start = 0; + cellRange.end = 0; + cellRange.points = nullptr; + } else { + // Get the range + ISGetPointRange(cellRange.is, &cellRange.start, &cellRange.end, &cellRange.points) >> utilities::PetscUtilities::checkError; + } +} + PetscErrorCode ablate::solver::Solver::DMPlexInsertBoundaryValues_Plex(DM dm, PetscDS prob, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) { PetscObject isZero; PetscInt numBd, b; diff --git a/src/solver/solver.hpp b/src/solver/solver.hpp index 8a55297a7..ae580a209 100644 --- a/src/solver/solver.hpp +++ b/src/solver/solver.hpp @@ -28,6 +28,9 @@ class Solver { // The region of this solver. const std::shared_ptr region; + // The region excluding any FVM ghost cells at boundaries + const std::shared_ptr regionMinusGhost; + protected: // an optional petscOptions that is used for this solver PetscOptions petscOptions; @@ -54,6 +57,9 @@ class Solver { /*** Set up mesh dependent initialization, this may be called multiple times if the mesh changes **/ virtual void Initialize() = 0; + /** Label all ghost cells so that the range can be returned without them. **/ + void SetupCellRangeWithoutGhost(); + /** string id for this solver **/ [[nodiscard]] inline const std::string& GetSolverId() const { return solverId; } @@ -124,6 +130,13 @@ class Solver { */ void GetRange(PetscInt depth, ablate::domain::Range& range) const { this->subDomain->GetRange(this->GetRegion(), depth, range); } + /** + * Get the range of cells defined for this solver without any FVM ghost (boundary) cells + * @param depth + * @param range + */ + void GetCellRangeWithoutGhost(ablate::domain::Range& cellRange) const; + /** * Restores the is and range - This needs to be removed and replaced with subDomain->RestoreRange * @param cellIS