diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 443bd712e..c5d8c9442 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -39,11 +39,19 @@ #include "blackoilmicpmodules.hh" #include #include +#include + +#include +#include +#include + #include #include #include +#include + namespace Opm { /*! * \ingroup BlackOilModel @@ -111,6 +119,32 @@ class BlackOilIntensiveQuantities using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities; using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities; + struct DirectionalMobility { + using array_type = std::array; + DirectionalMobility(const DirectionalMobility& other) + : mobilityX_{other.mobilityX_}, mobilityY_{other.mobilityY_}, mobilityZ_{other.mobilityZ_} {} + DirectionalMobility(const array_type& mX, const array_type& mY, const array_type& mZ) + : mobilityX_{mX}, mobilityY_{mY}, mobilityZ_{mZ} {} + DirectionalMobility() : mobilityX_{}, mobilityY_{}, mobilityZ_{} {} + array_type& getArray(int index) { + switch(index) { + case 0: + return mobilityX_; + case 1: + return mobilityY_; + case 2: + return mobilityZ_; + default: + throw std::runtime_error("Unexpected mobility array index"); + } + } + array_type mobilityX_; + array_type mobilityY_; + array_type mobilityZ_; + }; + using DirectionalMobilityPtr = Opm::Utility::CopyablePtr; + + public: using FluidState = BlackOilFluidState pC; const auto& materialParams = problem.materialLawParams(globalSpaceIdx); MaterialLaw::capillaryPressures(pC, materialParams, fluidState_); + BlackOilIntensiveQuantities::updateRelperms(mobility_, dirMob_, fluidState_, problem, materialParams, globalSpaceIdx); // oil is the reference phase for pressure if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv || priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) { @@ -242,11 +275,6 @@ public: fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx])); } - // calculate relative permeabilities. note that we store the result into the - // mobility_ class attribute. the division by the phase viscosity happens later. - MaterialLaw::relativePermeabilities(mobility_, materialParams, fluidState_); - Valgrind::CheckDefined(mobility_); - // update the Saturation functions for the blackoil solvent module. asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx); @@ -395,20 +423,31 @@ public: paramCache.updateAll(fluidState_); // compute the phase densities and transform the phase permeabilities into mobilities + int nmobilities = 1; + std::vector*> mobilities = {&mobility_}; + if (dirMob_) { + for (int i=0; i<3; i++) { + nmobilities += 1; + mobilities.push_back(&(dirMob_->getArray(i))); + } + } for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) continue; - const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx); fluidState_.setInvB(phaseIdx, b); - const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx); - if (enableExtbo && phaseIdx == oilPhaseIdx) - mobility_[phaseIdx] /= asImp_().oilViscosity(); - else if (enableExtbo && phaseIdx == gasPhaseIdx) - mobility_[phaseIdx] /= asImp_().gasViscosity(); - else - mobility_[phaseIdx] /= mu; + for (int i = 0; imobilityX_[phaseIdx]; + case Dir::YPlus: + return dirMob_->mobilityY_[phaseIdx]; + case Dir::ZPlus: + return dirMob_->mobilityZ_[phaseIdx]; + default: + throw std::runtime_error("Unexpected face direction"); + } + } + else { + return mobility_[phaseIdx]; + } + + } + /*! * \copydoc ImmiscibleIntensiveQuantities::porosity */ @@ -588,6 +648,46 @@ private: friend BlackOilBrineIntensiveQuantities; friend BlackOilMICPIntensiveQuantities; + template + static void updateRelperms( + std::array &mobility, + DirectionalMobilityPtr &dirMob, + FluidState &fluidState, + const Problem& problem, + const MaterialLawParams& materialParams, + unsigned globalSpaceIdx) + { + // calculate relative permeabilities. note that we store the result into the + // mobility_ class attribute. the division by the phase viscosity happens later. + MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState); + Valgrind::CheckDefined(mobility); + const auto* materialLawManager = problem.materialLawManagerPtr(); + if (materialLawManager && materialLawManager->hasDirectionalRelperms()) { + auto satnumIdx = materialLawManager->satnumRegionIdx(globalSpaceIdx); + using Dir = FaceDir::DirEnum; + constexpr int ndim = 3; + dirMob = std::make_unique(); + Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus}; + for (int i = 0; igetKrnumSatIdx(globalSpaceIdx, facedirs[i]); + auto& mob_array = dirMob->getArray(i); + if (krnumSatIdx != satnumIdx) { + // This hack is also used by StandardWell_impl.hpp:getMobilityEval() to temporarily use a different + // satnum index for a cell + const auto& paramsCell = materialLawManager->connectionMaterialLawParams(krnumSatIdx, globalSpaceIdx); + MaterialLaw::relativePermeabilities(mob_array, paramsCell, fluidState); + // reset the cell's satnum index back to the original + materialLawManager->connectionMaterialLawParams(satnumIdx, globalSpaceIdx); + } + else { + // Copy the default (non-directional dependent) mobility + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + mob_array[phaseIdx] = mobility[phaseIdx]; + } + } + } + } + } Implementation& asImp_() { return *static_cast(this); } @@ -597,6 +697,23 @@ private: Evaluation porosity_; Evaluation rockCompTransMultiplier_; std::array mobility_; + + // Instead of writing a custom copy constructor and a custom assignment operator just to handle + // the dirMob_ unique ptr member variable when copying BlackOilIntensiveQuantites (see for example + // updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we write the below + // custom wrapper class CopyablePtr which wraps the unique ptr and makes it copyable. + // + // The advantage of this approach is that we avoid having to call all the base class copy constructors and + // assignment operators explicitly (which is needed when writing the custom copy constructor and assignment + // operators) which could become a maintenance burden. For example, when adding a new base class (if that should + // be needed sometime in the future) to BlackOilIntensiveQuantites we could forget to update the copy + // constructor and assignment operators. + // + // We want each copy of the BlackOilIntensiveQuantites to be unique, (TODO: why?) so we have to make a copy + // of the unique_ptr each time we copy construct or assign to it from another BlackOilIntensiveQuantites. + // (On the other hand, if a copy could share the ptr with the original, a shared_ptr could be used instead and the + // wrapper would not be needed) + DirectionalMobilityPtr dirMob_; }; } // namespace Opm diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 4b54c332a..e5fc3242b 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -38,6 +38,8 @@ #include "blackoildiffusionmodule.hh" #include "blackoilmicpmodules.hh" #include +#include + namespace Opm { /*! @@ -198,7 +200,8 @@ public: const IntensiveQuantities& intQuantsIn, const IntensiveQuantities& intQuantsEx, const Scalar trans, - const Scalar faceArea) + const Scalar faceArea, + const FaceDir::DirEnum facedir) { flux = 0.0; Scalar Vin = problem.model().dofTotalVolume(globalIndexIn); @@ -233,7 +236,8 @@ public: distZ * g, thpres, trans, - faceArea); + faceArea, + facedir); } // This function demonstrates compatibility with the ElementContext-based interface. @@ -264,6 +268,11 @@ public: const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx); Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); Scalar faceArea = scvf.area(); + const auto& materialLawManager = problem.materialLawManager(); + FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown; // Use an arbitrary + if (materialLawManager->hasDirectionalRelperms()) { + facedir = scvf.faceDirFromDirId(); + } Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); // estimate the gravity correction: for performance reasons we use a simplified @@ -295,7 +304,8 @@ public: distZ * g, thpres, trans, - faceArea); + faceArea, + facedir); } static void calculateFluxes_(RateVector& flux, @@ -308,7 +318,8 @@ public: const Scalar& distZg, const Scalar& thpres, const Scalar& trans, - const Scalar& faceArea) + const Scalar& faceArea, + const FaceDir::DirEnum facedir) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) @@ -346,9 +357,10 @@ public: darcyFlux = 0.0; // NB maybe we could drop calculations } else { if (globalUpIndex == globalIndexIn) - darcyFlux = pressureDifference * up.mobility(phaseIdx) * transMult * (-trans / faceArea); + darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult * (-trans / faceArea); else - darcyFlux = pressureDifference * (Toolbox::value(up.mobility(phaseIdx)) * Toolbox::value(transMult) * (-trans / faceArea)); + darcyFlux = pressureDifference * + (Toolbox::value(up.mobility(phaseIdx, facedir)) * Toolbox::value(transMult) * (-trans / faceArea)); } unsigned pvtRegionIdx = up.pvtRegionIndex(); diff --git a/opm/models/common/multiphasebaseproblem.hh b/opm/models/common/multiphasebaseproblem.hh index 85964a5c4..f0380b64a 100644 --- a/opm/models/common/multiphasebaseproblem.hh +++ b/opm/models/common/multiphasebaseproblem.hh @@ -40,7 +40,13 @@ #include namespace Opm { - +// TODO: This hack is used to be able compile blackoilintensitivequantities.hh (see the function updateRelperms()) when +// the problem is not an EclProblem. For example if the problem is a ReservoirBlackOilVcfvProblem, the problem will not +// have a materialLawManager pointer (as the EclProblem has). Since this class MuitPhaseBaseProblem (see below) is a parent +// class for both those problem types, we can solve this problem by forward declaring EclMaterialLawManager here +// and defining a method materialLawManagerPtr() here that returns a nullptr, but is overridden in EclProblem to +// return the real EclMaterialManager pointer. +template class EclMaterialLawManager; /*! * \ingroup Discretization * @@ -246,6 +252,13 @@ public: return dummy; } + // TODO: See the comment at the top of this file for the reason why we need this method + template + const ::Opm::EclMaterialLawManager* materialLawManagerPtr() const + { + return nullptr; + } + /*! * \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume. * diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index fbc33a5e1..809d37dee 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -35,6 +35,7 @@ #include #include +#include #include #include @@ -314,7 +315,8 @@ private: unsigned numCells = model.numTotalDof(); neighborInfo_.reserve(numCells, 6 * numCells); std::vector loc_nbinfo; - + const auto& materialLawManager = problem_().materialLawManager(); + using FaceDirection = FaceDir::DirEnum; ElementIterator elemIt = gridView_().template begin<0>(); const ElementIterator elemEndIt = gridView_().template end<0>(); for (; elemIt != elemEndIt; ++elemIt) { @@ -330,8 +332,14 @@ private: sparsityPattern[myIdx].insert(neighborIdx); if (dofIdx > 0) { const double trans = problem_().transmissibility(myIdx, neighborIdx); - const double area = stencil.interiorFace(dofIdx - 1).area(); - loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area}; + const auto scvfIdx = dofIdx - 1; + const auto& scvf = stencil.interiorFace(scvfIdx); + const double area = scvf.area(); + FaceDirection dirId = FaceDirection::Unknown; + if (materialLawManager->hasDirectionalRelperms()) { + dirId = scvf.faceDirFromDirId(); + } + loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area, dirId}; } } neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end()); @@ -410,7 +418,8 @@ private: } const IntensiveQuantities& intQuantsEx = *intQuantsExP; LocalResidual::computeFlux( - adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, nbInfo.trans, nbInfo.faceArea); + adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, + nbInfo.trans, nbInfo.faceArea, nbInfo.faceDirection); adres *= nbInfo.faceArea; setResAndJacobi(res, bMat, adres); residual_[globI] += res; @@ -467,6 +476,7 @@ private: unsigned int neighbor; double trans; double faceArea; + FaceDir::DirEnum faceDirection; }; SparseTable neighborInfo_; }; diff --git a/opm/models/discretization/ecfv/ecfvstencil.hh b/opm/models/discretization/ecfv/ecfvstencil.hh index 8e233a7ad..343692436 100644 --- a/opm/models/discretization/ecfv/ecfvstencil.hh +++ b/opm/models/discretization/ecfv/ecfvstencil.hh @@ -37,6 +37,8 @@ #include #include #include +#include +#include #include @@ -157,11 +159,15 @@ public: if (needNormal) (*normal_) = intersection.centerUnitOuterNormal(); - const auto& geometry = intersection.geometry(); if (needIntegrationPos) (*integrationPos_) = geometry.center(); area_ = geometry.volume(); + // TODO: facedir_ = intersection.boundaryId(); // This did not compile for some reason + // using indexInInside() as in ecltransmissibility.cc instead + // Note that indexInInside() returns -1 for NNC faces, that's the reason we + // cannot convert directly to a FaceDir::DirEnum (see FaceDir.hpp) + dirId_ = intersection.indexInInside(); // Legal values: -1, 0, 1, 2, 3, 4, 5 } /*! @@ -205,10 +211,43 @@ public: Scalar area() const { return area_; } + /*! + * \brief Returns the direction of the face + */ + int dirId() const + { + return dirId_; + } + + /*! + * \brief Returns the direction of the face + */ + FaceDir::DirEnum faceDirFromDirId() const + { + using Dir = FaceDir::DirEnum; + if (dirId_ == -1) { + OPM_THROW(std::runtime_error, "NNC faces does not have a face id"); + } + switch(dirId_) { + case 0: + case 1: + return Dir::XPlus; + case 2: + case 3: + return Dir::YPlus; + case 4: + case 5: + return Dir::ZPlus; + default: + OPM_THROW(std::runtime_error, "Unexpected face id" << dirId_); + } + } + private: ConditionalStorage integrationPos_; ConditionalStorage normal_; Scalar area_; + int dirId_; unsigned short exteriorIdx_; };