Simulation partially supported for CpGrid with LGRs

opm/simulators/flow/FlowProblem.hpp

@@ -667,7 +667,13 @@ public:
         }
         }
         bool isSubStep = !EWOMS_GET_PARAM(TypeTag, bool, EnableWriteAllSolutions) && !this->simulator().episodeWillBeOver();
-        eclWriter_->evalSummaryState(isSubStep);
+        // For CpGrid with LGRs, ecl/vtk output is not supported yet.
+        const auto& grid =  this->simulator().vanguard().gridView().grid();
+        using GridType =  std::remove_cv_t< typename std::remove_reference<decltype(grid)>::type>;
+        bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
+        if ( !isCpGrid || (this->simulator().vanguard().gridView().grid().maxLevel()==0)) {
+            eclWriter_->evalSummaryState(isSubStep);
+        }
 
         int episodeIdx = this->episodeIndex();
 
@@ -2074,7 +2080,7 @@ protected:
         this->referencePorosity_[/*timeIdx=*/0].resize(numDof);
 
         const auto& fp = eclState.fieldProps();
-        const std::vector<double> porvData = fp.porv(false);
+        const std::vector<double> porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV");
         for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
             Scalar poreVolume = porvData[dofIdx];
 

opm/simulators/flow/GenericCpGridVanguard.cpp

@@ -35,6 +35,7 @@
 
 #include <opm/input/eclipse/Schedule/Schedule.hpp>
 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/LgrCollection.hpp>
 
 #include <opm/simulators/utils/ParallelEclipseState.hpp>
 #include <opm/simulators/utils/ParallelSerialization.hpp>
@@ -370,6 +371,8 @@ void GenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(Eclips
         global_porv = eclState.fieldProps().porv(true);
         OpmLog::info("\nProcessing grid");
     }
+
+    // --- Create grid ---
     OPM_TIMEBLOCK(createGrids);
 #if HAVE_MPI
     this->grid_ = std::make_unique<Dune::CpGrid>(FlowGenericVanguard::comm());
@@ -408,6 +411,17 @@ void GenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(Eclips
 
     cartesianIndexMapper_ = std::make_unique<CartesianIndexMapper>(*grid_);
 
+    // --- Add LGRs and update Leaf Grid View ---
+    // Check if input file contains Lgrs.
+    const auto& lgrs = eclState.getLgrs();
+    const auto lgrsSize = lgrs.size();
+    // If there are lgrs, create the grid with them, and update the leaf grid view.
+    if (lgrsSize)
+    {
+        OpmLog::info("\nAdding LGRs to the grid and updating its leaf grid view");
+        this->addLgrsUpdateLeafView(lgrs, lgrsSize);
+    }
+
 #if HAVE_MPI
     {
         const bool has_numerical_aquifer = eclState.aquifer().hasNumericalAquifer();
@@ -428,6 +442,7 @@ void GenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(Eclips
     }
 #endif
 
+    // --- Copy grid with LGRs to equilGrid_ ---
     // We use separate grid objects: one for the calculation of the initial
     // condition via EQUIL and one for the actual simulation. The reason is
     // that the EQUIL code is allergic to distributed grids and the
@@ -461,6 +476,29 @@ void GenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(Eclips
     eclState.pruneDeactivatedAquiferConnections(removed_cells);
 }
 
+template<class ElementMapper, class GridView, class Scalar>
+void GenericCpGridVanguard<ElementMapper,GridView,Scalar>::addLgrsUpdateLeafView(const LgrCollection& lgrCollection, const int lgrsSize)
+{
+    std::vector<std::array<int,3>> cells_per_dim_vec;
+    std::vector<std::array<int,3>> startIJK_vec;
+    std::vector<std::array<int,3>> endIJK_vec;
+    std::vector<std::string> lgrName_vec;
+    cells_per_dim_vec.reserve(lgrsSize);
+    startIJK_vec.reserve(lgrsSize);
+    endIJK_vec.reserve(lgrsSize);
+    lgrName_vec.reserve(lgrsSize);
+    for (int lgr = 0; lgr < lgrsSize; ++lgr)
+    {
+        const auto lgrCarfin = lgrCollection.getLgr(lgr);
+        cells_per_dim_vec.push_back({lgrCarfin.NX()/(lgrCarfin.I2() +1 - lgrCarfin.I1()),
+                lgrCarfin.NY()/(lgrCarfin.J2() +1 - lgrCarfin.J1()), lgrCarfin.NZ()/(lgrCarfin.K2() +1 - lgrCarfin.K1())});
+        startIJK_vec.push_back({lgrCarfin.I1(), lgrCarfin.J1(), lgrCarfin.K1()});
+        endIJK_vec.push_back({lgrCarfin.I2()+1, lgrCarfin.J2()+1, lgrCarfin.K2()+1});
+        lgrName_vec.emplace_back(lgrCarfin.NAME());
+    }
+    this->grid_->addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgrName_vec);
+};
+
 template<class ElementMapper, class GridView, class Scalar>
 void GenericCpGridVanguard<ElementMapper,GridView,Scalar>::
 doFilterConnections_(Schedule& schedule)

opm/simulators/flow/GenericCpGridVanguard.hpp

@@ -45,6 +45,7 @@ namespace Opm {
     class Schedule;
     class Well;
     class ParallelEclipseState;
+    class LgrCollection;
 }
 
 namespace Opm {
@@ -193,6 +194,7 @@ protected:
     void allocCartMapper();
 
     void doCreateGrids_(EclipseState& eclState);
+    void addLgrsUpdateLeafView(const LgrCollection& lgrCollection, const int lgrsSize);
 
     virtual void allocTrans() = 0;
     virtual double getTransmissibility(unsigned I, unsigned J) const = 0;

opm/simulators/flow/Transmissibility_impl.hpp

@@ -308,7 +308,7 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
             }
 
             const auto& outsideElem = intersection.outside();
-            unsigned outsideElemIdx = elemMapper.index(outsideElem);
+            unsigned outsideElemIdx = elemMapper.index(outsideElem);  
 
             // Get the Cartesian indices of the origen cells (parent or equivalent cell on level zero), for CpGrid with LGRs.
             // For genral grids and no LGRs, get the usual Cartesian Index.
@@ -802,7 +802,7 @@ updateFromEclState_(bool global)
             if(grid_.maxLevel()>0) {
                 OPM_THROW(std::invalid_argument, "Calculations on TRANX/TRANY/TRANZ arrays are not support with LGRS, yet.");
             }
-            fp->apply_tran(*key, *it);
+        fp->apply_tran(*key, *it);
         }
     }
 
@@ -974,9 +974,56 @@ computeFaceProperties(const Intersection& intersection,
                       /*isCpGrid=*/std::true_type) const
 {
     int faceIdx = intersection.id();
-    faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx);
-    faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx);
-    faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
+
+    if(grid_.maxLevel() == 0) {
+        faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
+        faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
+        faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
+    }
+    else {
+        if ((intersection.inside().level() != intersection.outside().level())) {
+
+            // For CpGrid with LGRs, intersection laying on the boundary of an LGR, having two neighboring cells:
+            // one coarse neighboring cell and one refined neighboring cell, we get the corresponding parent
+            // intersection (from level 0), and use the center of the parent intersection for the coarse
+            // neighboring cell.
+
+            // Get parent intersection and its geometry
+            const auto& parentIntersection = grid_.getParentIntersectionFromLgrBoundaryFace(intersection);
+            const auto& parentIntersectionGeometry = parentIntersection.geometry();
+
+            // For the coarse neighboring cell, take the center of the parent intersection.
+            // For the refined neighboring cell, take the 'usual' center. 
+            faceCenterInside =  (intersection.inside().level() == 0) ? parentIntersectionGeometry.center() :
+                grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
+            faceCenterOutside = (intersection.outside().level() == 0) ?  parentIntersectionGeometry.center() :
+                grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
+
+            // For some computations, it seems to be benefitial to replace the actual area of the refined face, by
+            // the area of its parent face. 
+            // faceAreaNormal = parentIntersection.centerUnitOuterNormal();
+            // faceAreaNormal *= parentIntersectionGeometry.volume();
+
+            /// Alternatively, the actual area of the refined face can be computed as follows:
+            faceAreaNormal = intersection.centerUnitOuterNormal();
+            faceAreaNormal *= intersection.geometry().volume();
+        }
+        else {
+            assert(intersection.inside().level() == intersection.outside().level());
+        
+            faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection);
+            faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection);
+        
+            // When the CpGrid has LGRs, we compute the face area normal differently.
+            if (intersection.inside().level() > 0) {  // remove intersection.inside().level() > 0
+                faceAreaNormal = intersection.centerUnitOuterNormal();
+                faceAreaNormal *= intersection.geometry().volume();
+            }
+            else {
+                faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
+            }
+        }
+    }
 }
 
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>