Change PORV computation for mix CpGrid

OPM · Dec 20, 2024 · e7a2b28 · e7a2b28
1 parent e714fe4
commit e7a2b28
Show file tree

Hide file tree

Showing 4 changed files with 64 additions and 5 deletions.
diff --git a/opm/grid/LookUpData.hh b/opm/grid/LookUpData.hh
@@ -93,7 +93,15 @@ public:
     operator()(const EntityType& elem, const std::vector<FieldPropType>& fieldProp) const;
 
     /// \brief: Get field property of type double from field properties manager by name.
-    std::vector<double> assignFieldPropsDoubleOnLeaf(const FieldPropsManager& fieldPropsManager,
+    template<typename GridType>
+    typename std::enable_if_t<!std::is_same_v<GridType,Dune::CpGrid>, std::vector<double>>
+    assignFieldPropsDoubleOnLeaf(const FieldPropsManager& fieldPropsManager,
+                                                     const std::string& propString) const;
+
+    /// \brief: Get field property of type double from field properties manager by name.
+    template<typename GridType>
+    typename std::enable_if_t<std::is_same_v<GridType,Dune::CpGrid>, std::vector<double>>
+    assignFieldPropsDoubleOnLeaf(const FieldPropsManager& fieldPropsManager,
                                                      const std::string& propString) const;
 
     /// \brief: Get field property of type int from field properties manager by name.
@@ -303,7 +311,9 @@ Opm::LookUpData<Grid,GridView>::operator()(const EntityType& elem,
 }
 
 template<typename Grid, typename GridView>
-std::vector<double> Opm::LookUpData<Grid,GridView>::assignFieldPropsDoubleOnLeaf(const FieldPropsManager& fieldPropsManager,
+template<typename GridType>
+typename std::enable_if_t<!std::is_same_v<GridType,Dune::CpGrid>, std::vector<double>>
+Opm::LookUpData<Grid,GridView>::assignFieldPropsDoubleOnLeaf(const FieldPropsManager& fieldPropsManager,
                                                                                  const std::string& propString) const
 {
     std::vector<double> fieldPropOnLeaf;
@@ -338,6 +348,50 @@ std::vector<double> Opm::LookUpData<Grid,GridView>::assignFieldPropsDoubleOnLeaf
     return fieldPropOnLeaf;
 }
 
+
+template<typename Grid, typename GridView>
+template<typename GridType>
+typename std::enable_if_t<std::is_same_v<GridType,Dune::CpGrid>, std::vector<double>>
+Opm::LookUpData<Grid,GridView>::assignFieldPropsDoubleOnLeaf(const FieldPropsManager& fieldPropsManager,
+                                                                                 const std::string& propString) const
+{
+    std::vector<double> fieldPropOnLeaf;
+    unsigned int numElements = gridView_.size(0);
+    fieldPropOnLeaf.resize(numElements);
+    const auto& fieldProp = fieldPropsManager.get_double(propString);
+    if ( (propString == "PORV") && (gridView_.grid().maxLevel() > 0)) {
+        // PORV poreVolume. LGRs supported (so far) only for CpGrid.
+        // For CpGrid with LGRs, poreVolume of a cell on the leaf grid view which has a parent cell on level 0,
+        // is computed as  porv[parent] * leafCellVolume / parentCellVolume. In this way, the sum of the pore
+        // volume of a parent cell coincides with the sum of the pore volume of its children.
+        for (const auto& element : elements(gridView_)) {
+            const auto& elemIdx = this-> elemMapper_.index(element);
+            const auto& fieldPropIdx = this->getFieldPropIdx<Grid>(elemIdx); // gets parentIdx (or (lgr)levelIdx) for CpGrid with LGRs
+            if (element.hasFather()) {
+                // const auto fatherVolume = element.father().geometry().volume();
+                // Compute total amount of siblings (total amount of refined child cells of father cell)
+                const auto& cells_per_dim = gridView_.grid().currentData()[element.level()]-> getCellsPerDim();
+                const auto& total_children = cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2];
+                // const auto& elemVolume = element.geometry().volume();
+                fieldPropOnLeaf[elemIdx] = fieldProp[fieldPropIdx] / total_children;
+                    //fieldProp[fieldPropIdx] * elemVolume / fatherVolume;
+            }
+            else {
+                fieldPropOnLeaf[elemIdx] = fieldProp[fieldPropIdx];
+            }
+        }
+    }
+    else {
+        for (const auto& element : elements(gridView_)) {
+            const auto& elemIdx = this-> elemMapper_.index(element);
+            const auto& fieldPropIdx = this->getFieldPropIdx<Grid>(elemIdx); // gets parentIdx (or (lgr)levelIdx) for CpGrid with LGRs
+            fieldPropOnLeaf[elemIdx] = fieldProp[fieldPropIdx];
+        }
+    }
+    return fieldPropOnLeaf;
+}
+
+
 template<typename Grid, typename GridView>
 template<typename IntType>
 std::vector<IntType> Opm::LookUpData<Grid,GridView>::assignFieldPropsIntOnLeaf(const FieldPropsManager& fieldPropsManager,

diff --git a/opm/grid/cpgrid/CpGridData.hpp b/opm/grid/cpgrid/CpGridData.hpp
@@ -536,6 +536,11 @@ class CpGridData
         return *level_data_ptr_;
     }
 
+    std::array<int,3> getCellsPerDim() const
+    {
+        return cells_per_dim_;
+    }
+
 
     /// @brief Refine a single cell and return a shared pointer of CpGridData type.
     ///

diff --git a/tests/cpgrid/lookupdataCpGrid_test.cpp b/tests/cpgrid/lookupdataCpGrid_test.cpp
@@ -303,13 +303,13 @@ void fieldProp_check(const Dune::CpGrid& grid, Opm::EclipseGrid eclGrid, const s
     const Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>
         mapper(leaf_view, Dune::mcmgElementLayout());
 
-    const auto& poroOnLeaf = lookUpData.assignFieldPropsDoubleOnLeaf(fpm, "PORO");
+    const auto& poroOnLeaf = lookUpData.template assignFieldPropsDoubleOnLeaf<Dune::CpGrid>(fpm, "PORO");
     const auto& poroOnLeafCart = lookUpCartesianData.assignFieldPropsDoubleOnLeaf(fpm, "PORO");
 
     const auto& eqlnumOnLeaf = lookUpData.assignFieldPropsIntOnLeaf<int>(fpm, "EQLNUM", true);
     const auto& eqlnumOnLeafCart = lookUpCartesianData.assignFieldPropsIntOnLeaf<int>(fpm, "EQLNUM", true);
 
-    const auto& porvOnLeaf = lookUpData.assignFieldPropsDoubleOnLeaf(fpm, "PORV");
+    const auto& porvOnLeaf = lookUpData.template assignFieldPropsDoubleOnLeaf<Dune::CpGrid>(fpm, "PORV");
 
     for (const auto& elem : elements(leaf_view))
     {

diff --git a/tests/test_lookupdata_polyhedral.cpp b/tests/test_lookupdata_polyhedral.cpp
@@ -191,7 +191,7 @@ void fieldProp_check(const Dune::PolyhedralGrid<3,3>& grid, Opm::EclipseGrid ecl
     // Element mapper
     const Dune::MultipleCodimMultipleGeomTypeMapper<GridView> mapper(leaf_view, Dune::mcmgElementLayout());
 
-    const auto& poroOnLeaf = lookUpData.assignFieldPropsDoubleOnLeaf(fpm, "PORO");
+    const auto& poroOnLeaf = lookUpData.template assignFieldPropsDoubleOnLeaf<Dune::PolyhedralGrid<3,3>>(fpm, "PORO");
     const auto& poroOnLeafCart = lookUpCartesianData.assignFieldPropsDoubleOnLeaf(fpm, "PORO");
 
     const auto& eqlnumOnLeaf = lookUpData.assignFieldPropsIntOnLeaf<int>(fpm, "EQLNUM", true);