From d1604b4279e25f6a9da29e29adb7a39e7865166c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 2 Feb 2026 12:02:52 +0100 Subject: [PATCH] Store rock fraction in problem --- opm/simulators/flow/FlowGenericProblem.hpp | 5 ++- .../flow/FlowGenericProblem_impl.hpp | 18 -------- opm/simulators/flow/FlowProblem.hpp | 42 +++++++++++++++++++ 3 files changed, 46 insertions(+), 19 deletions(-) diff --git a/opm/simulators/flow/FlowGenericProblem.hpp b/opm/simulators/flow/FlowGenericProblem.hpp index 5365634333..a739545066 100644 --- a/opm/simulators/flow/FlowGenericProblem.hpp +++ b/opm/simulators/flow/FlowGenericProblem.hpp @@ -146,7 +146,9 @@ public: * i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume * multiplier. This is to avoid negative rock volume for pvmult*porosity > 1 */ - Scalar rockFraction(unsigned elementIdx, unsigned timeIdx) const; + Scalar rockFraction(unsigned elementIdx, unsigned timeIdx) const { + return rockFraction_[timeIdx][elementIdx]; + }; /*! * \brief Sets the porosity of an element @@ -318,6 +320,7 @@ protected: static inline std::string briefDescription_; std::array, 2> referencePorosity_; + std::array, 2> rockFraction_; std::vector pvtnum_; std::vector satnum_; diff --git a/opm/simulators/flow/FlowGenericProblem_impl.hpp b/opm/simulators/flow/FlowGenericProblem_impl.hpp index d57c5cfd2c..0a51e581c9 100644 --- a/opm/simulators/flow/FlowGenericProblem_impl.hpp +++ b/opm/simulators/flow/FlowGenericProblem_impl.hpp @@ -338,24 +338,6 @@ porosity(unsigned globalSpaceIdx, unsigned timeIdx) const return this->referencePorosity_[timeIdx][globalSpaceIdx]; } -template -typename FlowGenericProblem::Scalar -FlowGenericProblem:: -rockFraction(unsigned elementIdx, unsigned timeIdx) const -{ - // For the energy equation, we need the volume of the rock. - // The volume of the rock is computed by rockFraction * geometric volume of the element. - // The reference porosity is defined as porosity * ntg * pore-volume-multiplier. - // A common practice in reservoir simulation is to use large pore-volume-multipliers in boundary cells - // to model boundary conditions other than no-flow. This may result in reference porosities that are larger than 1. - // A simple (1-reference porosity) * geometric volume of the element may give unphysical results. - // We therefore instead consider the pore-volume-multiplier as a volume multiplier. The rock fraction is thus given by - // (1 - porosity * ntg) * pore-volume-multiplier = (1 - porosity * ntg) * reference porosity / (porosity * ntg) - const auto ntg = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(), "NTG", elementIdx); - const auto poro_eff = ntg * this->lookUpData_.fieldPropDouble(eclState_.fieldProps(), "PORO", elementIdx); - return (1 - poro_eff) * referencePorosity(elementIdx, timeIdx) / poro_eff; -} - template template void FlowGenericProblem:: diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index bfc4f4579c..38d8f3c303 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -331,6 +331,8 @@ public: transmissibilities_.update(true, TransUpdateQuantities::All, equilGridToGrid); this->referencePorosity_[1] = this->referencePorosity_[0]; updateReferencePorosity_(); + this->rockFraction_[1] = this->rockFraction_[0]; + updateRockFraction_(); updatePffDofData_(); this->model().linearizer().updateDiscretizationParameters(); } @@ -1400,6 +1402,12 @@ protected: this->referencePorosity_[1] = this->referencePorosity_[0]; //////////////////////////////// + //////////////////////////////// + // rock fraction + updateRockFraction_(); + this->rockFraction_[1] = this->rockFraction_[0]; + //////////////////////////////// + //////////////////////////////// // fluid-matrix interactions (saturation functions; relperm/capillary pressure) materialLawManager_ = std::make_shared(); @@ -1452,6 +1460,40 @@ protected: } } + void updateRockFraction_() { + + const bool solveEnergyEquation = (energyModuleType == EnergyModules::FullyImplicitThermal || + energyModuleType == EnergyModules::SequentialImplicitThermal); + if (!solveEnergyEquation) + return; + + const auto& simulator = this->simulator(); + const auto& vanguard = simulator.vanguard(); + const auto& eclState = vanguard.eclState(); + + std::size_t numDof = this->model().numGridDof(); + this->rockFraction_[/*timeIdx=*/0].resize(numDof); + // For the energy equation, we need the volume of the rock. + // The volume of the rock is computed by rockFraction * geometric volume of the element. + // The reference porosity is defined as porosity * ntg * pore-volume-multiplier. + // A common practice in reservoir simulation is to use large pore-volume-multipliers in boundary cells + // to model boundary conditions other than no-flow. This may result in reference porosities that are larger than 1. + // A simple (1-reference porosity) * geometric volume of the element may give unphysical results. + // We therefore instead consider the pore-volume-multiplier as a volume multiplier. The rock fraction is thus given by + // (1 - porosity * ntg) * pore-volume-multiplier = (1 - porosity * ntg) * reference porosity / (porosity * ntg) + const auto& fp = eclState.fieldProps(); + const std::vector poroData = this->fieldPropDoubleOnLeafAssigner_()(fp, "PORO"); + const std::vector ntgData = this->fieldPropDoubleOnLeafAssigner_()(fp, "NTG"); + + for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { + const auto ntg = ntgData[dofIdx]; + const auto poro_eff = ntg * poroData[dofIdx]; + const int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx); + const auto rock_fraction = (1 - poro_eff) * this->referencePorosity_[/*timeIdx=*/0][sfcdofIdx] / poro_eff; + this->rockFraction_[/*timeIdx=*/0][sfcdofIdx] = rock_fraction; + } + } + virtual void readInitialCondition_() { // TODO: whether we should move this to FlowProblemBlackoil