From 5e1a37e49733a6b68db71bc43164c4c0e11e109f Mon Sep 17 00:00:00 2001 From: Antonella Ritorto Date: Mon, 29 Jan 2024 15:05:52 +0100 Subject: [PATCH] Simulation partially supported for CpGrid with LGRs --- opm/simulators/flow/FlowProblem.hpp | 10 +++- opm/simulators/flow/GenericCpGridVanguard.cpp | 38 +++++++++++++ opm/simulators/flow/GenericCpGridVanguard.hpp | 2 + opm/simulators/flow/Transmissibility_impl.hpp | 57 +++++++++++++++++-- 4 files changed, 100 insertions(+), 7 deletions(-) diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index 89ec8fae3..e622203ae 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/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::type>; + bool isCpGrid = std::is_same_v; + 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 porvData = fp.porv(false); + const std::vector porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV"); for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { Scalar poreVolume = porvData[dofIdx]; diff --git a/opm/simulators/flow/GenericCpGridVanguard.cpp b/opm/simulators/flow/GenericCpGridVanguard.cpp index b823d8363..8202849f4 100644 --- a/opm/simulators/flow/GenericCpGridVanguard.cpp +++ b/opm/simulators/flow/GenericCpGridVanguard.cpp @@ -35,6 +35,7 @@ #include #include +#include #include #include @@ -370,6 +371,8 @@ void GenericCpGridVanguard::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(FlowGenericVanguard::comm()); @@ -408,6 +411,17 @@ void GenericCpGridVanguard::doCreateGrids_(Eclips cartesianIndexMapper_ = std::make_unique(*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::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::doCreateGrids_(Eclips eclState.pruneDeactivatedAquiferConnections(removed_cells); } +template +void GenericCpGridVanguard::addLgrsUpdateLeafView(const LgrCollection& lgrCollection, const int lgrsSize) +{ + std::vector> cells_per_dim_vec; + std::vector> startIJK_vec; + std::vector> endIJK_vec; + std::vector 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 void GenericCpGridVanguard:: doFilterConnections_(Schedule& schedule) diff --git a/opm/simulators/flow/GenericCpGridVanguard.hpp b/opm/simulators/flow/GenericCpGridVanguard.hpp index 2c0b2e172..76db2eaf9 100644 --- a/opm/simulators/flow/GenericCpGridVanguard.hpp +++ b/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; diff --git a/opm/simulators/flow/Transmissibility_impl.hpp b/opm/simulators/flow/Transmissibility_impl.hpp index 0d48b0032..9cdf891e0 100644 --- a/opm/simulators/flow/Transmissibility_impl.hpp +++ b/opm/simulators/flow/Transmissibility_impl.hpp @@ -308,7 +308,7 @@ update(bool global, const std::function& 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