diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh index 30b5d0cad..57f04391b 100644 --- a/ebos/eclbasevanguard.hh +++ b/ebos/eclbasevanguard.hh @@ -209,6 +209,7 @@ public: protected: static const int dimension = Grid::dimension; using Element = typename GridView::template Codim<0>::Entity; + using CartesianIndexMapper = Dune::CartesianIndexMapper; public: @@ -637,6 +638,10 @@ public: const std::string& caseName() const { return caseName_; } + // TODO: revising this function later + const CartesianIndexMapper& cartesianMapper() const + { return asImp_().cartesianIndexMapper(); } + /*! * \brief Returns the number of logically Cartesian cells in each direction */ diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index 5c940f5c2..a5aef79cf 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -101,6 +101,7 @@ public: EQUIL::DeckDependent::InitialStateComputer initialState(materialLawManager, eclState, vanguard.gridView(), + vanguard.cartesianMapper(), simulator.problem().gravity()[dimWorld - 1]); // copy the result into the array of initial fluid states diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh index 6a544ea17..0d80356f9 100644 --- a/ebos/equil/initstateequil.hh +++ b/ebos/equil/initstateequil.hh @@ -1567,12 +1567,14 @@ class InitialStateComputer using Element = typename GridView::template Codim<0>::Entity; using ElementMapper = GetPropType; using Grid = GetPropType; + using CartesianIndexMapper = Dune::CartesianIndexMapper; public: template InitialStateComputer(MaterialLawManager& materialLawManager, const Opm::EclipseState& eclipseState, const GridView& gridView, + const CartesianIndexMapper& cartMapper, const double grav = Opm::unit::gravity, const bool applySwatInit = true) : temperature_(gridView.size(/*codim=*/0)), @@ -1582,7 +1584,8 @@ public: sat_(FluidSystem::numPhases, std::vector(gridView.size(/*codim=*/0))), rs_(gridView.size(/*codim=*/0)), - rv_(gridView.size(/*codim=*/0)) + rv_(gridView.size(/*codim=*/0)), + cartesianIndexMapper_(cartMapper) { //Check for presence of kw SWATINIT if (applySwatInit) { @@ -1594,7 +1597,7 @@ public: //Querry cell depth, cell top-bottom. // aquifer should enter here const auto& num_aquifers = eclipseState.aquifer().numericalAquifers(); - updateCellProps_(gridView); + updateCellProps_(gridView, num_aquifers); // Get the equilibration records. const std::vector rec = getEquil(eclipseState); @@ -1709,7 +1712,7 @@ public: // Compute pressures, saturations, rs and rv factors. const auto& comm = gridView.comm(); - calcPressSatRsRv(eqlmap, rec, materialLawManager, num_aquifers, comm, grav); + calcPressSatRsRv(eqlmap, rec, materialLawManager, gridView, num_aquifers, comm, grav); // Modify oil pressure in no-oil regions so that the pressures of present phases can // be recovered from the oil pressure and capillary relations. @@ -1772,12 +1775,14 @@ private: PVec sat_; Vec rs_; Vec rv_; + const CartesianIndexMapper& cartesianIndexMapper_; Vec swatInit_; Vec cellCenterDepth_; std::vector> cellZSpan_; std::vector> cellZMinMax_; - void updateCellProps_(const GridView& gridView) + void updateCellProps_(const GridView& gridView, + const NumericalAquifers& aquifer) { ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); int numElements = gridView.size(/*codim=*/0); @@ -1791,6 +1796,12 @@ private: const Element& element = *elemIt; const unsigned int elemIdx = elemMapper.index(element); cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); + const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); + if (aquifer.hasCell(cartIx)) { + const auto& aqu_cells = aquifer.aquiferCells(); + const auto& aqu_cell = aqu_cells.at(cartIx); + cellCenterDepth_[elemIdx] = aqu_cell.depth; + } cellZSpan_[elemIdx] = Details::cellZSpan(element); cellZMinMax_[elemIdx] = Details::cellZMinMax(element); } @@ -1811,6 +1822,7 @@ private: void calcPressSatRsRv(const RMap& reg, const std::vector< Opm::EquilRecord >& rec, MaterialLawManager& materialLawManager, + const GridView& gridView, const NumericalAquifers& aquifer, const Comm& comm, const double grav) @@ -1848,7 +1860,7 @@ private: const auto acc = eqreg.equilibrationAccuracy(); if (acc == 0) { // Centre-point method - this->equilibrateCellCentres(cells, eqreg, ptable, aquifer, psat); + this->equilibrateCellCentres(cells, eqreg, ptable, gridView, aquifer, psat); } else if (acc < 0) { // Horizontal subdivision @@ -1912,14 +1924,17 @@ private: void equilibrateCellCentres(const CellRange& cells, const EquilReg& eqreg, const PressTable& ptable, + const GridView& gridView, const NumericalAquifers& aquifer, PhaseSat& psat) { using CellPos = typename PhaseSat::Position; using CellID = std::remove_cv_t().cell)>>; + // TODO: might not needed + ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - this->cellLoop(cells, [this, &eqreg, &ptable, &aquifer, &psat] + this->cellLoop(cells, [this, &eqreg, &ptable, &elemMapper, &aquifer, &psat] (const CellID cell, Details::PhaseQuantityValue& pressures, Details::PhaseQuantityValue& saturations, @@ -1930,18 +1945,11 @@ private: cell, cellCenterDepth_[cell] }; -/* const size_t global_index = UgGridHelpers::globalCell(grid)[cell]; - - if (aquifer.hasCell(global_index)) { - const auto& aqu_cells = aquifer.aquiferCells(); - const auto& aqu_cell = aqu_cells.at(global_index); - pos = CellPos { - cell, aqu_cell.depth - }; - } */ - saturations = psat.deriveSaturations(pos, eqreg, ptable); -/* if (aquifer.hasCell(global_index)) { + + // TODO: not totally sure this is the cartesian Index + const auto global_index = cartesianIndexMapper_.cartesianIndex(cell); + if (aquifer.hasCell(global_index)) { saturations = {0.0, 0.0, 1.0}; const auto &aqu_cell = aquifer.getCell(global_index); std::ostringstream ss; @@ -1949,16 +1957,16 @@ private: << aqu_cell.K + 1 << " } OF NUMERICAL AQUIFER " << aqu_cell.aquifer_id << " with global_index " << global_index << ", " << "WATER SATURATION IS SET TO BE UNITY"; OpmLog::info(ss.str()); - } */ + } pressures = psat.correctedPhasePressures(); - /* if (aquifer.hasCell(global_index)) { + if (aquifer.hasCell(global_index)) { const auto &aqu_cell = aquifer.getCell(global_index); - // TODO: NOT totally sure what we should do here to empoly the pressure specified by AQUNUM + // TODO: NOT totally sure what we should do here to employ the pressure specified by AQUNUM if (aqu_cell.init_pressure > 0.) { const double pres = aqu_cell.init_pressure; pressures = {pres, pres, pres}; } - } */ + } const auto temp = this->temperature_[cell];