cellRangeWithoutGhost and Periodicity

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)

src/domain/boxMesh.cpp

@@ -36,7 +36,8 @@ DM ablate::domain::BoxMesh::CreateBoxDM(const std::string& name, std::vector<int
     // Make copy with PetscInt
     std::vector<PetscInt> 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;

src/domain/boxMeshBoundaryCells.cpp

@@ -66,7 +66,7 @@ DM ablate::domain::BoxMeshBoundaryCells::CreateBoxDM(const std::string& name, st
     // Make copy with PetscInt
     std::vector<PetscInt> 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<ablate::domain::modifiers::Modifier>, "postModifiers", "a list of domain modifiers to apply after ghost labeling"),
          ARG(std::vector<int>, "faces", "the number of faces in each direction"), ARG(std::vector<double>, "lower", "the lower bound of the mesh"),
          ARG(std::vector<double>, "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."));
+         OPT(ablate::parameters::Parameters, "options", "PETSc options specific to this dm.  Default value allows the dm to access global options."));

src/domain/domain.cpp

@@ -160,6 +160,12 @@ void ablate::domain::Domain::InitializeSubDomains(const std::vector<std::shared_
         subDomain->CreateSubDomainStructures();
     }
 
+    // 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;

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;
     }

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();

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

src/solver/solver.cpp

@@ -4,7 +4,7 @@
 #include "utilities/petscUtilities.hpp"
 
 ablate::solver::Solver::Solver(std::string solverId, std::shared_ptr<domain::Region> region, std::shared_ptr<parameters::Parameters> options)
-    : solverId(std::move(solverId)), region(std::move(region)), petscOptions(nullptr) {
+    : solverId(std::move(solverId)), region(std::move(region)), regionMinusGhost(std::make_shared<domain::Region>(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<ablate::domain::SubDomain> 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;

src/solver/solver.hpp

@@ -28,6 +28,9 @@ class Solver {
     // The region of this solver.
     const std::shared_ptr<domain::Region> region;
 
+    // The region excluding any FVM ghost cells at boundaries
+    const std::shared_ptr<domain::Region> 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