Preliminary support for mixed cell type flow (#501)

* working on multi-mesh type support

* working on multi-mesh type support

* tinny bug fix that took way too long

* test updates

* fixed bug in projection for FE

* formatting

* preliminary support for flow with mixed cell types

* updated tests for petsc update

* updated test for new petsc branch

* version bump

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.12)
+project(ablateLibrary VERSION 0.12.13)
 
 # Load the Required 3rd Party Libaries
 pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)

config/findXdmfGenerator.cmake

@@ -13,7 +13,7 @@ ELSE ()
     FetchContent_Declare(
             xdmfGeneratorLibrary
             GIT_REPOSITORY https://github.com/UBCHREST/XdmfGenerator.git
-            GIT_TAG v0.1.15
+            GIT_TAG v0.2.1
     )
     FetchContent_MakeAvailable(xdmfGeneratorLibrary)
     # Put the library into CHREST namespace

src/domain/boxMeshBoundaryCells.cpp

@@ -10,9 +10,9 @@
 
 ablate::domain::BoxMeshBoundaryCells::BoxMeshBoundaryCells(const std::string& name, const std::vector<std::shared_ptr<FieldDescriptor>>& fieldDescriptors,
                                                            std::vector<std::shared_ptr<modifiers::Modifier>> preModifiers, std::vector<std::shared_ptr<modifiers::Modifier>> postModifiers,
-                                                           std::vector<int> faces, const std::vector<double>& lower, const std::vector<double>& upper, bool simplex,
+                                                           std::vector<int> faces, const std::vector<double>& lower, const std::vector<double>& upper,
                                                            const std::shared_ptr<parameters::Parameters>& options)
-    : Domain(CreateBoxDM(name, std::move(faces), lower, upper, simplex), name, fieldDescriptors,
+    : Domain(CreateBoxDM(name, std::move(faces), lower, upper, false), name, fieldDescriptors,
              // We need to get the optional dm_plex_scale to determine bounds
              AddBoundaryModifiers(lower, upper, options ? options->Get("dm_plex_scale", 1.0) : 1.0, std::move(preModifiers), std::move(postModifiers)), options) {
     // make sure that dm_refine was not set
@@ -163,5 +163,5 @@ REGISTER(ablate::domain::Domain, ablate::domain::BoxMeshBoundaryCells,
          OPT(std::vector<ablate::domain::modifiers::Modifier>, "preModifiers", "a list of domain modifiers to apply before ghost labeling"),
          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(bool, "simplex", "sets if the elements/cells are simplex"),
+         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."));

src/domain/boxMeshBoundaryCells.hpp

@@ -38,7 +38,7 @@ class BoxMeshBoundaryCells : public Domain {
    public:
     BoxMeshBoundaryCells(const std::string& name, const std::vector<std::shared_ptr<FieldDescriptor>>& fieldDescriptors, std::vector<std::shared_ptr<modifiers::Modifier>> preModifiers,
                          std::vector<std::shared_ptr<modifiers::Modifier>> postModifiers, std::vector<int> faces, const std::vector<double>& lower, const std::vector<double>& upper,
-                         bool simplex = true, const std::shared_ptr<parameters::Parameters>& options = {});
+                         const std::shared_ptr<parameters::Parameters>& options = {});
 
     ~BoxMeshBoundaryCells() override;
 };

src/domain/domain.cpp

@@ -1,11 +1,11 @@
 #include "domain.hpp"
 #include <set>
-#include <typeinfo>
 #include <utility>
 #include "solver/solver.hpp"
 #include "subDomain.hpp"
 #include "utilities/demangler.hpp"
 #include "utilities/mpiUtilities.hpp"
+#include "utilities/petscSupport.hpp"
 #include "utilities/petscUtilities.hpp"
 
 ablate::domain::Domain::Domain(DM dmIn, std::string name, std::vector<std::shared_ptr<FieldDescriptor>> fieldDescriptorsIn, std::vector<std::shared_ptr<modifiers::Modifier>> modifiersIn,
@@ -140,7 +140,7 @@ std::shared_ptr<ablate::domain::SubDomain> ablate::domain::Domain::GetSubDomain(
     }
 }
 
-void ablate::domain::Domain::InitializeSubDomains(const std::vector<std::shared_ptr<solver::Solver>>& solvers, std::shared_ptr<ablate::domain::Initializer> initializations,
+void ablate::domain::Domain::InitializeSubDomains(const std::vector<std::shared_ptr<solver::Solver>>& solvers, const std::shared_ptr<ablate::domain::Initializer>& initializations,
                                                   const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& exactSolutions) {
     // determine the number of fields
     for (auto& solver : solvers) {
@@ -240,7 +240,7 @@ void ablate::domain::Domain::ProjectFieldFunctions(const std::vector<std::shared
                 ISDestroy(&regionIS) >> utilities::PetscUtilities::checkError;
             }
         } else {
-            DMProjectFunctionLocal(dm, time, fieldFunctionsPts.data(), fieldContexts.data(), INSERT_VALUES, locVec) >> utilities::PetscUtilities::checkError;
+            DMProjectFunctionLocalMixedCells(dm, time, fieldFunctionsPts.data(), fieldContexts.data(), INSERT_VALUES, locVec) >> utilities::PetscUtilities::checkError;
         }
     }
 

src/domain/domain.hpp

@@ -93,7 +93,7 @@ class Domain : private utilities::Loggable<Domain>, private ablate::utilities::N
      * @param solvers
      * @param initializations
      */
-    void InitializeSubDomains(const std::vector<std::shared_ptr<solver::Solver>>& solvers = {}, std::shared_ptr<ablate::domain::Initializer> initializations = {},
+    void InitializeSubDomains(const std::vector<std::shared_ptr<solver::Solver>>& solvers = {}, const std::shared_ptr<ablate::domain::Initializer>& initializations = {},
                               const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& = {});
 
     /**

src/domain/subDomain.cpp

@@ -2,6 +2,7 @@
 #include <set>
 #include <sstream>
 #include "utilities/mpiUtilities.hpp"
+#include "utilities/petscSupport.hpp"
 #include "utilities/petscUtilities.hpp"
 
 ablate::domain::SubDomain::SubDomain(Domain& domainIn, PetscInt dsNumber, const std::vector<std::shared_ptr<FieldDescription>>& allAuxFields)
@@ -200,7 +201,7 @@ PetscErrorCode ablate::domain::SubDomain::ProjectFieldFunctionsToSubDM(const std
         fieldFunctions[fieldId.subId] = fieldInitialization->GetSolutionField().GetPetscFunction();
     }
 
-    PetscCall(DMProjectFunctionLocal(subDM, time, &fieldFunctions[0], &fieldContexts[0], INSERT_VALUES, locVec));
+    PetscCall(DMProjectFunctionLocalMixedCells(subDM, time, &fieldFunctions[0], &fieldContexts[0], INSERT_VALUES, locVec));
 
     // push the results back to the global vector
     PetscCall(DMLocalToGlobal(dm, locVec, INSERT_VALUES, globVec));
@@ -739,7 +740,7 @@ void ablate::domain::SubDomain::ProjectFieldFunctionsToLocalVector(const std::ve
                 ISDestroy(&regionIS) >> utilities::PetscUtilities::checkError;
             }
         } else {
-            DMProjectFunctionLocal(dm, time, fieldFunctionsPts.data(), fieldContexts.data(), INSERT_VALUES, locVec) >> utilities::PetscUtilities::checkError;
+            DMProjectFunctionLocalMixedCells(dm, time, fieldFunctionsPts.data(), fieldContexts.data(), INSERT_VALUES, locVec) >> utilities::PetscUtilities::checkError;
         }
     }
 }

src/finiteVolume/stencils/faceStencilGenerator.hpp

@@ -14,7 +14,7 @@ class FaceStencilGenerator {
    public:
     virtual ~FaceStencilGenerator() = default;
 
-    virtual void Generate(PetscInt face, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region> solverRegion, DM cellDM, const PetscScalar* cellGeomArray,
+    virtual void Generate(PetscInt face, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region>& solverRegion, DM cellDM, const PetscScalar* cellGeomArray,
                           DM faceDM, const PetscScalar* faceGeomArray) = 0;
 };
 

src/finiteVolume/stencils/leastSquares.cpp

@@ -3,7 +3,7 @@
 #include "utilities/petscUtilities.hpp"
 
 void ablate::finiteVolume::stencil::LeastSquares::Generate(PetscInt face, ablate::finiteVolume::stencil::Stencil& stencil, const domain::SubDomain& subDomain,
-                                                           const std::shared_ptr<domain::Region> solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
+                                                           const std::shared_ptr<domain::Region>& solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
                                                            const PetscScalar* faceGeomArray) {
     auto dm = subDomain.GetDM();
     auto dim = subDomain.GetDimensions();

src/finiteVolume/stencils/leastSquares.hpp

@@ -16,7 +16,7 @@ class LeastSquares : public FaceStencilGenerator {
     PetscFV gradientCalculator = nullptr;
 
    public:
-    void Generate(PetscInt face, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region> solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
+    void Generate(PetscInt face, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region>& solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
                   const PetscScalar* faceGeomArray) override;
     ~LeastSquares() override;
 };

src/finiteVolume/stencils/leastSquaresAverage.cpp

@@ -3,7 +3,7 @@
 #include "utilities/petscUtilities.hpp"
 
 void ablate::finiteVolume::stencil::LeastSquaresAverage::Generate(PetscInt face, ablate::finiteVolume::stencil::Stencil& stencil, const domain::SubDomain& subDomain,
-                                                                  const std::shared_ptr<domain::Region> solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
+                                                                  const std::shared_ptr<domain::Region>& solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
                                                                   const PetscScalar* faceGeomArray) {
     auto dm = subDomain.GetDM();
     auto dim = subDomain.GetDimensions();
@@ -36,7 +36,7 @@ void ablate::finiteVolume::stencil::LeastSquaresAverage::Generate(PetscInt face,
     stencil.stencil.insert(stencil.stencil.end(), rightStencil.stencil.begin(), rightStencil.stencil.end());
     std::sort(stencil.stencil.begin(), stencil.stencil.end());
     stencil.stencil.erase(std::unique(stencil.stencil.begin(), stencil.stencil.end()), stencil.stencil.end());
-    stencil.stencilSize = stencil.stencil.size();
+    stencil.stencilSize = (PetscInt)stencil.stencil.size();
 
     // resize the weight arrays and init to zero
     stencil.weights.resize(stencil.stencilSize, 0.0);
@@ -67,7 +67,7 @@ void ablate::finiteVolume::stencil::LeastSquaresAverage::Generate(PetscInt face,
 }
 
 void ablate::finiteVolume::stencil::LeastSquaresAverage::ComputeNeighborCellStencil(PetscInt cell, ablate::finiteVolume::stencil::Stencil& stencil, const ablate::domain::SubDomain& subDomain,
-                                                                                    const std::shared_ptr<domain::Region> solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
+                                                                                    const std::shared_ptr<domain::Region>& solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
                                                                                     const PetscScalar* faceGeomArray) {
     // only add the cells if they are in the ds
     if (!subDomain.InRegion(cell)) {

src/finiteVolume/stencils/leastSquaresAverage.hpp

@@ -16,14 +16,14 @@ class LeastSquaresAverage : public FaceStencilGenerator {
     PetscFV gradientCalculator = nullptr;
 
    public:
-    void Generate(PetscInt face, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region> solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
+    void Generate(PetscInt face, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region>& solverRegion, DM cellDM, const PetscScalar* cellGeomArray, DM faceDM,
                   const PetscScalar* faceGeomArray) override;
     ~LeastSquaresAverage() override;
 
     /**
      * helper function to compute left or right stencil
      */
-    void ComputeNeighborCellStencil(PetscInt cell, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region> solverRegion, DM cellDM,
+    void ComputeNeighborCellStencil(PetscInt cell, Stencil& stencil, const domain::SubDomain& subDomain, const std::shared_ptr<domain::Region>& solverRegion, DM cellDM,
                                     const PetscScalar* cellGeomArray, DM faceDM, const PetscScalar* faceGeomArray);
 };
 

src/io/hdf5MultiFileSerializer.cpp

@@ -206,7 +206,11 @@ void ablate::io::Hdf5MultiFileSerializer::SaveMetadata(TS ts) const {
         // keep a back of the restart file incase writing fails
         if (exists(restartFilePath)) {
             auto backupRestartFilePath = rootOutputDirectory / "restart.bak";
-            std::filesystem::copy(restartFilePath, backupRestartFilePath, std::filesystem::copy_options::overwrite_existing);
+            try {
+                std::filesystem::copy(restartFilePath, backupRestartFilePath, std::filesystem::copy_options::overwrite_existing);
+            } catch (const std::filesystem::filesystem_error& error) {
+                std::cout << "Cannot create restart.bak file: " << error.what() << " Continuing without backup file." << std::endl;
+            }
         }
         std::ofstream restartFile;
         restartFile.open(restartFilePath);

src/io/hdf5Serializer.cpp

@@ -151,11 +151,11 @@ ablate::io::Hdf5Serializer::Hdf5ObjectSerializer::~Hdf5ObjectSerializer() {
         // If this is the root process generate the xdmf file
         PetscMPIInt rank;
         MPI_Comm_rank(PetscObjectComm(PetscObject(petscViewer)), &rank);
+        PetscViewerDestroy(&petscViewer) >> utilities::PetscUtilities::checkError;
+
         if (rank == 0 && !filePath.empty() && std::filesystem::exists(filePath)) {
             xdmfGenerator::Generate(filePath);
         }
-
-        PetscViewerDestroy(&petscViewer) >> utilities::PetscUtilities::checkError;
     }
 }
 

src/utilities/petscSupport.cpp

@@ -1398,3 +1398,33 @@ PetscErrorCode DMPlexCellGradFromCell(DM dm, const PetscInt c, Vec data, PetscIn
 
     PetscFunctionReturn(PETSC_SUCCESS);
 }
+
+PetscErrorCode DMProjectFunctionLocalMixedCells(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx), void **ctxs,
+                                                InsertMode mode, Vec localX) {
+    PetscFunctionBegin;
+    // Call once to cover the FE fields.
+    PetscCall(DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX));
+
+    // define the cells we would like to project over
+    const PetscInt PetscCellTypeCount = 12;
+    const PetscInt types[12] = {DM_POLYTOPE_POINT,
+                                DM_POLYTOPE_SEGMENT,
+                                DM_POLYTOPE_POINT_PRISM_TENSOR,
+                                DM_POLYTOPE_TRIANGLE,
+                                DM_POLYTOPE_QUADRILATERAL,
+                                DM_POLYTOPE_SEG_PRISM_TENSOR,
+                                DM_POLYTOPE_TETRAHEDRON,
+                                DM_POLYTOPE_HEXAHEDRON,
+                                DM_POLYTOPE_TRI_PRISM,
+                                DM_POLYTOPE_TRI_PRISM_TENSOR,
+                                DM_POLYTOPE_QUAD_PRISM_TENSOR,
+                                DM_POLYTOPE_PYRAMID};
+
+    // extract the current cell type label
+    DMLabel ctLabel;
+    PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
+
+    // project on to this mesh over all types
+    PetscCall(DMProjectFunctionLabelLocal(dm, time, ctLabel, PetscCellTypeCount, types, -1, NULL, funcs, ctxs, mode, localX));
+    PetscFunctionReturn(PETSC_SUCCESS);
+}

src/utilities/petscSupport.hpp

@@ -169,3 +169,9 @@ PetscErrorCode PetscSortedArrayComplement(const PetscInt na, const PetscInt a[],
  *    b - All integers in a and b
  */
 PetscErrorCode PetscSortedArrayCommon(const PetscInt na, const PetscInt a[], PetscInt *nb, PetscInt b[]);
+
+/**
+ * This is a copy of DMProjectFunctionLocal (https://petsc.org/main/manualpages/DM/DMProjectFunctionLocal/) but projects across all cells even with different cell types
+ * @return
+ */
+PetscErrorCode DMProjectFunctionLocalMixedCells(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);

tests/integrationTests/inputs/compressibleFlow/compressibleFlowCadExample/compressibleFlowCadExample.txt

@@ -58,8 +58,6 @@ MinMaxAvg velocity for timestep 0010:
 ResultFiles:
 domain.hdf5
 domain.xmf
-domain_1.xmf
 restart.rst
 vortexFlowField.hdf5
 vortexFlowField.xmf
-vortexFlowField_1.xmf

tests/integrationTests/inputs/compressibleFlow/compressibleFlowPgsLodi.yaml

@@ -37,7 +37,6 @@ timestepper:
     faces: [ 25, 25 ]
     lower: [ 0.0, -.5 ]
     upper: [ 1.0, .5 ]
-    simplex: false
     options:
       # pass in these options to petsc when setting up the domain.  Using an option list here prevents command line arguments from being seen.
       dm_refine: 0 # must be zero when using the BoxMeshBoundaryCells

tests/integrationTests/inputs/compressibleFlow/compressibleFlowVortexLodi.yaml

@@ -60,7 +60,6 @@ timestepper:
     faces: [ 25, 25]
     lower: [ 0.0, -.5 ]
     upper: [ 1.0, .5 ]
-    simplex: false
     # pass in these options to petsc when setting up the domain.  Using an option list here prevents command line arguments from being seen.
     options:
       dm_refine: 0 # must be zero when using the BoxMeshBoundaryCells

tests/integrationTests/inputs/compressibleFlow/compressibleFlowVortexLodiRestart.yaml

@@ -68,7 +68,6 @@ timestepper:
     faces: [ 25, 25]
     lower: [ 0.0, -.5 ]
     upper: [ 1.0, .5 ]
-    simplex: false
     # pass in these options to petsc when setting up the domain.  Using an option list here prevents command line arguments from being seen.
     options:
       dm_refine: 0 # must be zero when using the BoxMeshBoundaryCells

tests/integrationTests/inputs/compressibleFlow/compressibleSublimationPipe/compressibleSublimationPipe.yaml

@@ -34,7 +34,6 @@ timestepper:
     faces: [ 50, 10 ]
     lower: [ 0.0, 0.0 ]
     upper: [ .5, .1 ]
-    simplex: false
     # pass in these options to petsc when setting up the domain.  Using an option list here prevents command line arguments from being seen.
     options:
       dm_refine: 0 # must be zero when using the BoxMeshBoundaryCells

tests/integrationTests/inputs/compressibleFlow/mixedCellTypeTest/inHexMesh.csv

@@ -0,0 +1,4 @@
+time,Yi_N2,Yi_H2O,Yi_O2,
+0,0,0,0,
+(.*),(.*),(.*),(.*),<expects> ~ >0.44 >.20 >0.34
+(.*),(.*),(.*),(.*),<expects> ~ >0.7 <0.15 <0.2

tests/integrationTests/inputs/compressibleFlow/mixedCellTypeTest/inTriangleMesh.csv

@@ -0,0 +1,4 @@
+time,Yi_N2,Yi_H2O,Yi_O2,
+0,0,0,0,
+(.*),(.*),(.*),(.*),<expects> ~ <0.3 >0.2 >0.4
+(.*),(.*),(.*),(.*),<expects> ~ >0.5 <0.2 <0.3

tests/integrationTests/inputs/compressibleFlow/mixedCellTypeTest/mixedCellType2D.txt

@@ -0,0 +1,13 @@
+Timestep: 0000 time = 0        dt = (.*)<expects> ~
+Timestep: 0500 time = 0.02909  dt = (.*)<expects> ~
+Timestep: 1000 time = 0.05802  dt = (.*)<expects> ~
+ResultFiles:
+domain.hdf5
+domain.xmf
+inHexMesh.csv
+inTriangleMesh.csv
+mixedCellLatex.tex
+pgsLog
+restart.rst
+vortexFlowField.hdf5
+vortexFlowField.xmf