From 81cf436630233defcc5fcbecbf5952a7877050cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Thu, 18 Aug 2022 09:55:39 +0200 Subject: [PATCH 01/15] Add support for KRNUM --- .../blackoil/blackoilintensivequantities.hh | 86 ++++++++++++++++--- opm/models/discretization/ecfv/ecfvstencil.hh | 30 ++++++- 2 files changed, 101 insertions(+), 15 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 443bd712e..0674cba19 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -39,11 +39,18 @@ #include "blackoilmicpmodules.hh" #include #include +#include + +#include +#include + #include #include #include +#include + namespace Opm { /*! * \ingroup BlackOilModel @@ -77,6 +84,8 @@ class BlackOilIntensiveQuantities using Indices = GetPropType; using GridView = GetPropType; using FluxModule = GetPropType; + using EclMaterialLawManager = typename GetProp::EclMaterialLawManager; + using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams; enum { numEq = getPropValue() }; enum { enableSolvent = getPropValue() }; @@ -144,7 +153,6 @@ public: const auto& problem = elemCtx.problem(); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); - const auto& linearizationType = problem.model().linearizer().getLinearizationType(); unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx); Scalar RvMax = FluidSystem::enableVaporizedOil() @@ -226,6 +234,7 @@ public: std::array pC; const auto& materialParams = problem.materialLawParams(globalSpaceIdx); MaterialLaw::capillaryPressures(pC, materialParams, fluidState_); + updateRelperms(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 +251,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 +399,26 @@ public: paramCache.updateAll(fluidState_); // compute the phase densities and transform the phase permeabilities into mobilities + constexpr int nmobilities = 4; + std::array*,nmobilities> mobilities + = {&mobility_, &mobilityX_, &mobilityY_, &mobilityZ_}; 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; i; friend BlackOilMICPIntensiveQuantities; + void updateRelperms(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.materialLawManager(); + auto satnumIdx = materialLawManager->satnumRegionIdx(globalSpaceIdx); + using Dir = FaceDir::DirEnum; + constexpr int ndim = 3; + std::array*,ndim> relperms = {&mobilityX_, &mobilityY_, &mobilityZ_}; + Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus}; + for (int i = 0; igetKrnumSatIdx(globalSpaceIdx, facedirs[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(*relperms[i], 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) { + (*relperms[i])[phaseIdx] = mobility_[phaseIdx]; + } + } + } + } Implementation& asImp_() { return *static_cast(this); } @@ -597,6 +652,9 @@ private: Evaluation porosity_; Evaluation rockCompTransMultiplier_; std::array mobility_; + std::array mobilityX_; + std::array mobilityY_; + std::array mobilityZ_; }; } // namespace Opm diff --git a/opm/models/discretization/ecfv/ecfvstencil.hh b/opm/models/discretization/ecfv/ecfvstencil.hh index 8e233a7ad..3b691365d 100644 --- a/opm/models/discretization/ecfv/ecfvstencil.hh +++ b/opm/models/discretization/ecfv/ecfvstencil.hh @@ -157,11 +157,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 +209,34 @@ public: Scalar area() const { return area_; } + /*! + * \brief Returns the direction of the face + */ + FaceDir::DirEnum dirId() const + { + if (dirId_ == -1) { + OPM_THROW(std::runtime_error, "NNC faces does not have a face id"); + } + switch(dirId_) { + case 0: + case 1: + return FaceDir::XPlus; + case 2: + case 3: + return FaceDir::YPlus; + case 4: + case 5: + return FaceDir::ZPlus; + default: + OPM_THROW(std::runtime_error, "Unexpected face id" << dirId_); + } + } + private: ConditionalStorage integrationPos_; ConditionalStorage normal_; Scalar area_; + int dirId_; unsigned short exteriorIdx_; }; From ca78d271ce335c646e345cfb1d6611bb0a7a0c5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Tue, 6 Sep 2022 19:22:23 +0200 Subject: [PATCH 02/15] Only store directional mobilities if needed --- .../blackoil/blackoilintensivequantities.hh | 102 ++++++++++++------ 1 file changed, 68 insertions(+), 34 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 0674cba19..8eab7576e 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -120,6 +120,28 @@ class BlackOilIntensiveQuantities using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities; using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities; + struct DirectionalMobility { + using array_type = std::array; + 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 index00"); + } + } + array_type mobilityX_; + array_type mobilityY_; + array_type mobilityZ_; + }; + public: using FluidState = BlackOilFluidState*,nmobilities> mobilities - = {&mobility_, &mobilityX_, &mobilityY_, &mobilityZ_}; + 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; @@ -545,16 +572,22 @@ public: const Evaluation& mobility(unsigned phaseIdx, FaceDir::DirEnum facedir) const { using Dir = FaceDir::DirEnum; - switch(facedir) { - case Dir::XPlus: - return mobilityX_[phaseIdx]; - case Dir::YPlus: - return mobilityY_[phaseIdx]; - case Dir::ZPlus: - return mobilityZ_[phaseIdx]; - default: - throw std::runtime_error("Unexpected face direction"); + if (dirMob_) { + switch(facedir) { + case Dir::XPlus: + return dirMob_->mobilityX_[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]; + } + } /*! @@ -620,25 +653,28 @@ private: MaterialLaw::relativePermeabilities(mobility_, materialParams, fluidState_); Valgrind::CheckDefined(mobility_); const auto& materialLawManager = problem.materialLawManager(); - auto satnumIdx = materialLawManager->satnumRegionIdx(globalSpaceIdx); - using Dir = FaceDir::DirEnum; - constexpr int ndim = 3; - std::array*,ndim> relperms = {&mobilityX_, &mobilityY_, &mobilityZ_}; - Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus}; - for (int i = 0; igetKrnumSatIdx(globalSpaceIdx, facedirs[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(*relperms[i], 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) { - (*relperms[i])[phaseIdx] = mobility_[phaseIdx]; + if (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]; + } } } } @@ -652,9 +688,7 @@ private: Evaluation porosity_; Evaluation rockCompTransMultiplier_; std::array mobility_; - std::array mobilityX_; - std::array mobilityY_; - std::array mobilityZ_; + std::shared_ptr dirMob_; }; } // namespace Opm From fe7b4153505a2f812ad0d25fd6a88f9604c4c5f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Thu, 8 Sep 2022 10:32:35 +0200 Subject: [PATCH 03/15] Fix typo --- opm/models/blackoil/blackoilintensivequantities.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 8eab7576e..49e1685df 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -134,7 +134,7 @@ class BlackOilIntensiveQuantities case 2: return mobilityZ_; default: - throw std::runtime_error("Unexpected mobility array index00"); + throw std::runtime_error("Unexpected mobility array index"); } } array_type mobilityX_; From 7eeca2801899323a0906a57406786ba29102fc4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Wed, 14 Sep 2022 15:33:36 +0200 Subject: [PATCH 04/15] Add more support for directional relperms Adds support for directional relperms for the tpfa linearizer. --- .../blackoil/blackoillocalresidualtpfa.hh | 24 ++++++++++++++----- .../discretization/common/tpfalinearizer.hh | 18 ++++++++++---- 2 files changed, 32 insertions(+), 10 deletions(-) diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 4b54c332a..5353e5874 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.dirId(); + } 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/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index fbc33a5e1..705943545 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.dirId(); + } + 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_; }; From ee8e5651375b818d16d202ee194a4f7471bd61b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 16 Sep 2022 09:17:39 +0200 Subject: [PATCH 05/15] Use unique_ptr for dirMob_ instead of shared_ptr To make each copy of the IQ unique, we change dirMob_ to be a shared_ptr instead of a unique_ptr. --- .../blackoil/blackoilintensivequantities.hh | 43 +++++++++++++++++-- 1 file changed, 40 insertions(+), 3 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 49e1685df..7c2031aa6 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -162,9 +162,46 @@ public: } } - BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities& other) = default; + BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities& other) { + if (other.dirMob_) { + dirMob_ = std::make_unique( + other.dirMob_->mobilityX_, + other.dirMob_->mobilityY_, + other.dirMob_->mobilityY_ + ); + } + }; - BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default; + BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) { + fluidState_ = other.fluidState_; + referencePorosity_ = other.referencePorosity_; + porosity_ = other.porosity_; + rockCompTransMultiplier_ = other.rockCompTransMultiplier_; + mobility_ = other.mobility_; + if (other.dirMob_) { + dirMob_ = std::make_unique( + other.dirMob_->mobilityX_, + other.dirMob_->mobilityY_, + other.dirMob_->mobilityY_ + ); + } + else { + dirMob_.release(); // release any memory, and assign nullptr + } + // call assignment operators in the parent classes + GetPropType::operator=(other); + GetPropType::FluxIntensiveQuantities::operator=(other); + BlackOilDiffusionIntensiveQuantities() >::operator=(other); + BlackOilSolventIntensiveQuantities::operator=(other); + BlackOilExtboIntensiveQuantities::operator=(other); + BlackOilPolymerIntensiveQuantities::operator=(other); + BlackOilFoamIntensiveQuantities::operator=(other); + BlackOilBrineIntensiveQuantities::operator=(other); + BlackOilEnergyIntensiveQuantities::operator=(other); + BlackOilMICPIntensiveQuantities::operator=(other); + + return *this; + }; /*! * \copydoc IntensiveQuantities::update @@ -688,7 +725,7 @@ private: Evaluation porosity_; Evaluation rockCompTransMultiplier_; std::array mobility_; - std::shared_ptr dirMob_; + std::unique_ptr dirMob_; }; } // namespace Opm From dc6141af13c6822da658a3da546b0c1bc1a1289d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 16 Sep 2022 09:20:48 +0200 Subject: [PATCH 06/15] Make updateRelperms() a static function To aid readability, we make updateRelperms() a static function. --- .../blackoil/blackoilintensivequantities.hh | 22 ++++++++++++------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 7c2031aa6..1cdb6c9e1 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -293,7 +293,7 @@ public: std::array pC; const auto& materialParams = problem.materialLawParams(globalSpaceIdx); MaterialLaw::capillaryPressures(pC, materialParams, fluidState_); - updateRelperms(problem, materialParams, globalSpaceIdx); + 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) { @@ -683,34 +683,40 @@ private: friend BlackOilBrineIntensiveQuantities; friend BlackOilMICPIntensiveQuantities; - void updateRelperms(const Problem& problem, const MaterialLawParams& materialParams, unsigned globalSpaceIdx) + static void updateRelperms( + std::array &mobility, + std::unique_ptr &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_); + MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState); + Valgrind::CheckDefined(mobility); const auto& materialLawManager = problem.materialLawManager(); if (materialLawManager->hasDirectionalRelperms()) { auto satnumIdx = materialLawManager->satnumRegionIdx(globalSpaceIdx); using Dir = FaceDir::DirEnum; constexpr int ndim = 3; - dirMob_ = std::make_unique(); + 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); + 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_); + 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]; + mob_array[phaseIdx] = mobility[phaseIdx]; } } } From 38add2ed66b60aecbba8c45be35f0b3d39517a8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 16 Sep 2022 09:22:17 +0200 Subject: [PATCH 07/15] Rename dirId() to faceDirFromDirId() --- opm/models/blackoil/blackoillocalresidualtpfa.hh | 2 +- opm/models/discretization/common/tpfalinearizer.hh | 2 +- opm/models/discretization/ecfv/ecfvstencil.hh | 10 +++++++++- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 5353e5874..e5fc3242b 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -271,7 +271,7 @@ public: const auto& materialLawManager = problem.materialLawManager(); FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown; // Use an arbitrary if (materialLawManager->hasDirectionalRelperms()) { - facedir = scvf.dirId(); + facedir = scvf.faceDirFromDirId(); } Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 705943545..809d37dee 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -337,7 +337,7 @@ private: const double area = scvf.area(); FaceDirection dirId = FaceDirection::Unknown; if (materialLawManager->hasDirectionalRelperms()) { - dirId = scvf.dirId(); + dirId = scvf.faceDirFromDirId(); } loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area, dirId}; } diff --git a/opm/models/discretization/ecfv/ecfvstencil.hh b/opm/models/discretization/ecfv/ecfvstencil.hh index 3b691365d..94322fc4e 100644 --- a/opm/models/discretization/ecfv/ecfvstencil.hh +++ b/opm/models/discretization/ecfv/ecfvstencil.hh @@ -212,7 +212,15 @@ public: /*! * \brief Returns the direction of the face */ - FaceDir::DirEnum dirId() const + int dirId() const + { + return dirId_; + } + + /*! + * \brief Returns the direction of the face + */ + FaceDir::DirEnum faceDirFromDirId() const { if (dirId_ == -1) { OPM_THROW(std::runtime_error, "NNC faces does not have a face id"); From d15c5b51c8b96e4d0116264541a6b16c947cc403 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Mon, 19 Sep 2022 01:42:09 +0200 Subject: [PATCH 08/15] Use a copyable unique pointer Use a copyable unique pointer instead of writing custom copy constructor and assignment operator. --- .../blackoil/blackoilintensivequantities.hh | 103 ++++++++++-------- 1 file changed, 60 insertions(+), 43 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 1cdb6c9e1..5ad66322a 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -122,7 +122,9 @@ class BlackOilIntensiveQuantities struct DirectionalMobility { using array_type = std::array; - DirectionalMobility(const array_type& mX, const array_type& mY, const array_type& mZ) + 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) { @@ -142,6 +144,59 @@ class BlackOilIntensiveQuantities array_type mobilityZ_; }; + // 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 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) + template + class CopyablePtr { + public: + CopyablePtr() : ptr_(nullptr) {} + CopyablePtr(const CopyablePtr& other) { + if (other) { // other does not contain a nullptr + ptr_ = std::make_unique(other.get()); + } + // else {ptr_ = nullptr;} // this is the default construction value + } + // assignment operator + CopyablePtr& operator=(const CopyablePtr& other) { + if (other) { + ptr_ = std::make_unique(other.get()); + } + else { + ptr_ = nullptr; + } + return *this; + } + // assign directly from a unique_ptr + CopyablePtr& operator=(std::unique_ptr&& uptr) { + ptr_ = std::move(uptr); + return *this; + } + // member access operator + T* operator->() const {return ptr_.get(); } + // boolean context operator + explicit operator bool() const noexcept { + return ptr_ ? true : false; + } + // get a reference to the stored value + T& get() const {return *ptr_;} + T* release() const {return ptr_.release();} + private: + std::unique_ptr ptr_; + }; + public: using FluidState = BlackOilFluidState( - other.dirMob_->mobilityX_, - other.dirMob_->mobilityY_, - other.dirMob_->mobilityY_ - ); - } - }; - - BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) { - fluidState_ = other.fluidState_; - referencePorosity_ = other.referencePorosity_; - porosity_ = other.porosity_; - rockCompTransMultiplier_ = other.rockCompTransMultiplier_; - mobility_ = other.mobility_; - if (other.dirMob_) { - dirMob_ = std::make_unique( - other.dirMob_->mobilityX_, - other.dirMob_->mobilityY_, - other.dirMob_->mobilityY_ - ); - } - else { - dirMob_.release(); // release any memory, and assign nullptr - } - // call assignment operators in the parent classes - GetPropType::operator=(other); - GetPropType::FluxIntensiveQuantities::operator=(other); - BlackOilDiffusionIntensiveQuantities() >::operator=(other); - BlackOilSolventIntensiveQuantities::operator=(other); - BlackOilExtboIntensiveQuantities::operator=(other); - BlackOilPolymerIntensiveQuantities::operator=(other); - BlackOilFoamIntensiveQuantities::operator=(other); - BlackOilBrineIntensiveQuantities::operator=(other); - BlackOilEnergyIntensiveQuantities::operator=(other); - BlackOilMICPIntensiveQuantities::operator=(other); - - return *this; - }; + BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default; /*! * \copydoc IntensiveQuantities::update @@ -685,7 +702,7 @@ private: static void updateRelperms( std::array &mobility, - std::unique_ptr &dirMob, + CopyablePtr &dirMob, FluidState &fluidState, const Problem& problem, const MaterialLawParams& materialParams, @@ -731,7 +748,7 @@ private: Evaluation porosity_; Evaluation rockCompTransMultiplier_; std::array mobility_; - std::unique_ptr dirMob_; + CopyablePtr dirMob_; }; } // namespace Opm From d9c6ecf79dc3b8effa414d66f78d283c6c9d0304 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Mon, 19 Sep 2022 12:03:22 +0200 Subject: [PATCH 09/15] Add missing include file --- opm/models/discretization/ecfv/ecfvstencil.hh | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/models/discretization/ecfv/ecfvstencil.hh b/opm/models/discretization/ecfv/ecfvstencil.hh index 94322fc4e..37538b5d3 100644 --- a/opm/models/discretization/ecfv/ecfvstencil.hh +++ b/opm/models/discretization/ecfv/ecfvstencil.hh @@ -37,6 +37,7 @@ #include #include #include +#include #include @@ -222,19 +223,20 @@ public: */ 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 FaceDir::XPlus; + return Dir::XPlus; case 2: case 3: - return FaceDir::YPlus; + return Dir::YPlus; case 4: case 5: - return FaceDir::ZPlus; + return Dir::ZPlus; default: OPM_THROW(std::runtime_error, "Unexpected face id" << dirId_); } From cedb92c6d241ae70cef6aa44bb7aa3919ea6e8be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Mon, 19 Sep 2022 18:15:04 +0200 Subject: [PATCH 10/15] Make get() return a pointer instead of a reference This is more in line with how std::unique_ptr works. --- opm/models/blackoil/blackoilintensivequantities.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 5ad66322a..1a8692c2f 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -165,14 +165,14 @@ class BlackOilIntensiveQuantities CopyablePtr() : ptr_(nullptr) {} CopyablePtr(const CopyablePtr& other) { if (other) { // other does not contain a nullptr - ptr_ = std::make_unique(other.get()); + ptr_ = std::make_unique(*other.get()); } // else {ptr_ = nullptr;} // this is the default construction value } // assignment operator CopyablePtr& operator=(const CopyablePtr& other) { if (other) { - ptr_ = std::make_unique(other.get()); + ptr_ = std::make_unique(*other.get()); } else { ptr_ = nullptr; @@ -191,7 +191,7 @@ class BlackOilIntensiveQuantities return ptr_ ? true : false; } // get a reference to the stored value - T& get() const {return *ptr_;} + T* get() const {return ptr_.get();} T* release() const {return ptr_.release();} private: std::unique_ptr ptr_; From 769577c8aa7f3104cf3a907bea92aad6f0535937 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Mon, 19 Sep 2022 18:48:33 +0200 Subject: [PATCH 11/15] Explicitly initialize member in copy constructor For clarity, explicitly default construct the unique_ptr in the copy constructor --- opm/models/blackoil/blackoilintensivequantities.hh | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 1a8692c2f..293b6a3a8 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -167,7 +167,9 @@ class BlackOilIntensiveQuantities if (other) { // other does not contain a nullptr ptr_ = std::make_unique(*other.get()); } - // else {ptr_ = nullptr;} // this is the default construction value + else { + ptr_ = nullptr; + } } // assignment operator CopyablePtr& operator=(const CopyablePtr& other) { From e143c9ed13c7d52f6c05aa1001bd6317f5a15444 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Tue, 20 Sep 2022 01:27:59 +0200 Subject: [PATCH 12/15] Move the CopyablePtr class to a separate file Moved the CopyablePtr template class to opm-common, see file opm/utility/CopyablePtr.hpp in opm-common --- .../blackoil/blackoilintensivequantities.hh | 76 +++++-------------- 1 file changed, 20 insertions(+), 56 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 293b6a3a8..f8947585b 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -43,6 +43,7 @@ #include #include +#include #include @@ -143,61 +144,8 @@ class BlackOilIntensiveQuantities array_type mobilityY_; array_type mobilityZ_; }; + using DirectionalMobilityPtr = Opm::Utility::CopyablePtr; - // 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 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) - template - class CopyablePtr { - public: - CopyablePtr() : ptr_(nullptr) {} - CopyablePtr(const CopyablePtr& other) { - if (other) { // other does not contain a nullptr - ptr_ = std::make_unique(*other.get()); - } - else { - ptr_ = nullptr; - } - } - // assignment operator - CopyablePtr& operator=(const CopyablePtr& other) { - if (other) { - ptr_ = std::make_unique(*other.get()); - } - else { - ptr_ = nullptr; - } - return *this; - } - // assign directly from a unique_ptr - CopyablePtr& operator=(std::unique_ptr&& uptr) { - ptr_ = std::move(uptr); - return *this; - } - // member access operator - T* operator->() const {return ptr_.get(); } - // boolean context operator - explicit operator bool() const noexcept { - return ptr_ ? true : false; - } - // get a reference to the stored value - T* get() const {return ptr_.get();} - T* release() const {return ptr_.release();} - private: - std::unique_ptr ptr_; - }; public: using FluidState = BlackOilFluidState &mobility, - CopyablePtr &dirMob, + DirectionalMobilityPtr &dirMob, FluidState &fluidState, const Problem& problem, const MaterialLawParams& materialParams, @@ -750,7 +698,23 @@ private: Evaluation porosity_; Evaluation rockCompTransMultiplier_; std::array mobility_; - CopyablePtr dirMob_; + + // 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 From 90a3ae8940cd064fa07741126bcf2e97f7e69357 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Tue, 20 Sep 2022 21:59:00 +0200 Subject: [PATCH 13/15] Include missing header --- opm/models/discretization/ecfv/ecfvstencil.hh | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/models/discretization/ecfv/ecfvstencil.hh b/opm/models/discretization/ecfv/ecfvstencil.hh index 37538b5d3..343692436 100644 --- a/opm/models/discretization/ecfv/ecfvstencil.hh +++ b/opm/models/discretization/ecfv/ecfvstencil.hh @@ -37,6 +37,7 @@ #include #include #include +#include #include #include From a941ba49c09c7b813571d949c97f2fed1639e597 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Tue, 20 Sep 2022 22:09:04 +0200 Subject: [PATCH 14/15] Try to fix: no type named EclMaterialLawManager Try to fix jenkins build error: no type named EclMaterialLawManager --- opm/models/blackoil/blackoilintensivequantities.hh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index f8947585b..e1bb0d226 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -85,8 +85,6 @@ class BlackOilIntensiveQuantities using Indices = GetPropType; using GridView = GetPropType; using FluxModule = GetPropType; - using EclMaterialLawManager = typename GetProp::EclMaterialLawManager; - using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams; enum { numEq = getPropValue() }; enum { enableSolvent = getPropValue() }; @@ -650,6 +648,7 @@ private: friend BlackOilBrineIntensiveQuantities; friend BlackOilMICPIntensiveQuantities; + template static void updateRelperms( std::array &mobility, DirectionalMobilityPtr &dirMob, From f36aa67804522132c4ce56d4349f6527c60e11b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 23 Sep 2022 00:05:09 +0200 Subject: [PATCH 15/15] Return a nullptr to the EclMaterialLawManager Define a method materialLawManagerPtr() that returns a nullpointer to EclMaterialLawManager, but that can be overridden in derived classes e.g. EclProblem --- .../blackoil/blackoilintensivequantities.hh | 4 ++-- opm/models/common/multiphasebaseproblem.hh | 15 ++++++++++++++- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index e1bb0d226..c5d8c9442 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -661,8 +661,8 @@ private: // mobility_ class attribute. the division by the phase viscosity happens later. MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState); Valgrind::CheckDefined(mobility); - const auto& materialLawManager = problem.materialLawManager(); - if (materialLawManager->hasDirectionalRelperms()) { + const auto* materialLawManager = problem.materialLawManagerPtr(); + if (materialLawManager && materialLawManager->hasDirectionalRelperms()) { auto satnumIdx = materialLawManager->satnumRegionIdx(globalSpaceIdx); using Dir = FaceDir::DirEnum; constexpr int ndim = 3; 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. *