From b65f5e9808ba034f59b7b0fa7ac0d913288704bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 14 Sep 2022 10:33:47 +0200 Subject: [PATCH 1/7] Add new overload for transmissibilityBoundary() without element context. --- ebos/eclproblem.hh | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index e5812607c..058632ef4 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -1340,6 +1340,15 @@ public: return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx); } + /*! + * \brief Direct access to a boundary transmissibility. + */ + Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, + const unsigned boundaryFaceIdx) const + { + return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx); + } + /*! * \copydoc EclTransmissiblity::thermalHalfTransmissibility */ From 68e1479cafa3a4d5c870a4cf3cdfac3b8d1cb175 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 14 Sep 2022 10:36:14 +0200 Subject: [PATCH 2/7] [WIP] Refactor calculateBoundaryGradients_() Not addressing solvent yet. --- ebos/eclfluxmodule.hh | 91 ++++++++++++++++++++++++++++--------------- 1 file changed, 59 insertions(+), 32 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 138673586..8f1c013cb 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -418,34 +418,61 @@ protected: unsigned timeIdx, const FluidState& exFluidState) { + const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx); + Scalar faceArea = scvf.area(); + Scalar zEx = scvf.integrationPos()[dimWorld - 1]; const auto& problem = elemCtx.problem(); + calculateBoundaryGradients_(problem, + elemCtx.globalSpaceIndex(0, timeIdx), + elemCtx.intensiveQuantities(0, timeIdx), + scvfIdx, + timeIdx, + faceArea, + zEx, + exFluidState, + upIdx_, + dnIdx_, + volumeFlux_, + pressureDifference_); + } + +public: + /*! + * \brief Update the required gradients for boundary faces + */ + template + static void calculateBoundaryGradients_(const Problem& problem, + const unsigned globalSpaceIdx, + const IntensiveQuantities& intQuantsIn, + const unsigned bfIdx, + const unsigned timeIdx, + const double faceArea, + const double zEx, + const FluidState& exFluidState, + short (&upIdx)[numPhases], + short (&dnIdx)[numPhases], + EvaluationContainer& volumeFlux, + EvaluationContainer& pressureDifference) + { bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions(); if (!enableBoundaryMassFlux) return; - const auto& stencil = elemCtx.stencil(timeIdx); - const auto& scvf = stencil.boundaryFace(scvfIdx); - unsigned interiorDofIdx = scvf.interiorIndex(); - - Scalar trans = problem.transmissibilityBoundary(elemCtx, scvfIdx); - Scalar faceArea = scvf.area(); + Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx); // estimate the gravity correction: for performance reasons we use a simplified // approach for this flux module that assumes that gravity is constant and always // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) - Scalar g = elemCtx.problem().gravity()[dimWorld - 1]; - - const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx); + Scalar g = problem.gravity()[dimWorld - 1]; // this is quite hacky because the dune grid interface does not provide a // cellCenterDepth() method (so we ask the problem to provide it). The "good" // solution would be to take the Z coordinate of the element centroids, but since // ECL seems to like to be inconsistent on that front, it needs to be done like // here... - Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx); - Scalar zEx = scvf.integrationPos()[dimWorld - 1]; + Scalar zIn = problem.dofCenterDepth(globalSpaceIdx); // the distances from the DOF's depths. (i.e., the additional depth of the // exterior DOF) @@ -465,62 +492,62 @@ protected: Evaluation pressureExterior = exFluidState.pressure(phaseIdx); pressureExterior += rhoAvg*(distZ*g); - pressureDifference_[phaseIdx] = pressureExterior - pressureInterior; + pressureDifference[phaseIdx] = pressureExterior - pressureInterior; // decide the upstream index for the phase. for this we make sure that the // degree of freedom which is regarded upstream if both pressures are equal // is always the same: if the pressure is equal, the DOF with the lower // global index is regarded to be the upstream one. - if (pressureDifference_[phaseIdx] > 0.0) { - upIdx_[phaseIdx] = -1; - dnIdx_[phaseIdx] = interiorDofIdx; + const unsigned interiorDofIdx = 0; // Valid only for cell-centered FV. + if (pressureDifference[phaseIdx] > 0.0) { + upIdx[phaseIdx] = -1; + dnIdx[phaseIdx] = interiorDofIdx; } else { - upIdx_[phaseIdx] = interiorDofIdx; - dnIdx_[phaseIdx] = -1; + upIdx[phaseIdx] = interiorDofIdx; + dnIdx[phaseIdx] = -1; } Evaluation transModified = trans; - unsigned upstreamIdx = upstreamIndex_(phaseIdx); - if (upstreamIdx == interiorDofIdx) { + if (upIdx[phaseIdx] == interiorDofIdx) { // this is slightly hacky because in the automatic differentiation case, it // only works for the element centered finite volume method. for ebos this // does not matter, though. - const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); + const auto& up = intQuantsIn; // deal with water induced rock compaction const double transMult = Toolbox::value(up.rockCompTransMultiplier()); transModified *= transMult; - volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); + volumeFlux[phaseIdx] = + pressureDifference[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); - if (enableSolvent && phaseIdx == gasPhaseIdx) - asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-transModified/faceArea)); + // TODO: Figure out if this did have any effect. It should? + // if (enableSolvent && phaseIdx == gasPhaseIdx) + // asImp_().setSolventVolumeFlux( pressureDifference[phaseIdx]*up.solventMobility()*(-transModified/faceArea)); } else { // compute the phase mobility using the material law parameters of the // interior element. TODO: this could probably be done more efficiently - const auto& matParams = - elemCtx.problem().materialLawParams(elemCtx, - interiorDofIdx, - /*timeIdx=*/0); + const auto& matParams = problem.materialLawParams(globalSpaceIdx); std::array kr; MaterialLaw::relativePermeabilities(kr, matParams, exFluidState); const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx); - volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*mob*(-transModified/faceArea); + volumeFlux[phaseIdx] = + pressureDifference[phaseIdx]*mob*(-transModified/faceArea); // Solvent inflow is not yet supported - if (enableSolvent && phaseIdx == gasPhaseIdx) - asImp_().setSolventVolumeFlux(0.0); + // if (enableSolvent && phaseIdx == gasPhaseIdx) + // asImp_().setSolventVolumeFlux(0.0); } } } +protected: + /*! * \brief Update the volumetric fluxes for all fluid phases on the interior faces of the context */ From 5e08a9fc12591a5b7a1b61b8255bddf8589e553b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 16 Sep 2022 16:17:37 +0200 Subject: [PATCH 3/7] Make boundary condition data accessor function. --- ebos/eclproblem.hh | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 058632ef4..25018499e 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -2061,6 +2061,43 @@ public: return this->rockCompTransMultWc_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); } + std::pair boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) + { + if (!nonTrivialBoundaryConditions_ || directionId == -1) { + return std::make_pair(false, RateVector(0.0)); + } + const std::vector* flag = nullptr; + const std::vector* massrate = nullptr; + switch (directionId) { + case 0: + flag = &freebcXMinus_; + massrate = &massratebcXMinus_; + break; + case 1: + flag = &freebcX_; + massrate = &massratebcX_; + break; + case 2: + flag = &freebcYMinus_; + massrate = &massratebcYMinus_; + break; + case 3: + flag = &freebcY_; + massrate = &massratebcY_; + break; + case 4: + flag = &freebcZMinus_; + massrate = &massratebcZMinus_; + break; + case 5: + flag = &freebcZ_; + massrate = &massratebcZ_; + break; + } + const bool free = (*flag)[globalSpaceIdx]; + return { free, (*massrate)[globalSpaceIdx] }; + } + private: // update the parameters needed for DRSDT and DRVDT void updateCompositionChangeLimits_() From bc13f57e95ed5f29ba4d7962997913f2ce54722c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 22 Sep 2022 09:31:54 +0200 Subject: [PATCH 4/7] Simplify the boundaryCondition() method. --- ebos/eclproblem.hh | 34 ++++++++++------------------------ 1 file changed, 10 insertions(+), 24 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 25018499e..294a9419a 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -2063,39 +2063,25 @@ public: std::pair boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) { - if (!nonTrivialBoundaryConditions_ || directionId == -1) { - return std::make_pair(false, RateVector(0.0)); + if (!nonTrivialBoundaryConditions_) { + return { false, RateVector(0.0) }; } - const std::vector* flag = nullptr; - const std::vector* massrate = nullptr; switch (directionId) { case 0: - flag = &freebcXMinus_; - massrate = &massratebcXMinus_; - break; + return { freebcXMinus_[globalSpaceIdx], massratebcXMinus_[globalSpaceIdx] }; case 1: - flag = &freebcX_; - massrate = &massratebcX_; - break; + return { freebcX_[globalSpaceIdx], massratebcX_[globalSpaceIdx] }; case 2: - flag = &freebcYMinus_; - massrate = &massratebcYMinus_; - break; + return { freebcYMinus_[globalSpaceIdx], massratebcYMinus_[globalSpaceIdx] }; case 3: - flag = &freebcY_; - massrate = &massratebcY_; - break; + return { freebcY_[globalSpaceIdx], massratebcY_[globalSpaceIdx] }; case 4: - flag = &freebcZMinus_; - massrate = &massratebcZMinus_; - break; + return { freebcZMinus_[globalSpaceIdx], massratebcZMinus_[globalSpaceIdx] }; case 5: - flag = &freebcZ_; - massrate = &massratebcZ_; - break; + return { freebcZ_[globalSpaceIdx], massratebcZ_[globalSpaceIdx] }; + default: + return { false, RateVector(0.0) }; } - const bool free = (*flag)[globalSpaceIdx]; - return { free, (*massrate)[globalSpaceIdx] }; } private: From b7758aa7061dca05337a6bf6dab3b99d3f4c3215 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 22 Sep 2022 09:32:17 +0200 Subject: [PATCH 5/7] Use std::array rather than C arrays. --- ebos/eclfluxmodule.hh | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 8f1c013cb..fe8908ed1 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -43,6 +43,8 @@ #include #include +#include + namespace Opm { template @@ -212,8 +214,8 @@ protected: public: - static void volumeAndPhasePressureDifferences(short (&upIdx)[numPhases], - short (&dnIdx)[numPhases], + static void volumeAndPhasePressureDifferences(std::array& upIdx, + std::array& dnIdx, Evaluation (&volumeFlux)[numPhases], Evaluation (&pressureDifferences)[numPhases], const ElementContext& elemCtx, @@ -449,8 +451,8 @@ public: const double faceArea, const double zEx, const FluidState& exFluidState, - short (&upIdx)[numPhases], - short (&dnIdx)[numPhases], + std::array& upIdx, + std::array& dnIdx, EvaluationContainer& volumeFlux, EvaluationContainer& pressureDifference) { @@ -572,9 +574,9 @@ private: Evaluation pressureDifference_[numPhases]; // the local indices of the interior and exterior degrees of freedom - short upIdx_[numPhases]; - short dnIdx_[numPhases]; -}; + std::array upIdx_; + std::array dnIdx_; + }; } // namespace Opm From 04183dea312b79d695052fe7b7b93c32a39e17c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 23 Sep 2022 16:03:18 +0200 Subject: [PATCH 6/7] Ensure solvent boundary fluxes work. This does the job for the element-context-based linearizer, but not for the TpfaLinearizer, which when combined with the solvent model will require a bit more work for bcs. --- ebos/eclfluxmodule.hh | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index fe8908ed1..ed0bcf9a2 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -421,12 +421,15 @@ protected: const FluidState& exFluidState) { const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx); - Scalar faceArea = scvf.area(); - Scalar zEx = scvf.integrationPos()[dimWorld - 1]; + const Scalar faceArea = scvf.area(); + const Scalar zEx = scvf.integrationPos()[dimWorld - 1]; const auto& problem = elemCtx.problem(); + const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx); + const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx); + calculateBoundaryGradients_(problem, - elemCtx.globalSpaceIndex(0, timeIdx), - elemCtx.intensiveQuantities(0, timeIdx), + globalSpaceIdx, + intQuantsIn, scvfIdx, timeIdx, faceArea, @@ -436,6 +439,20 @@ protected: dnIdx_, volumeFlux_, pressureDifference_); + + // Treating solvent here and not in the static method, since that would require more + // extensive refactoring. It means that the TpfaLinearizer will not support bcs for solvent until this is + // addressed. + if constexpr (enableSolvent) { + if (upIdx_[gasPhaseIdx] == 0) { + const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx); + const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier()); + const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea); + asImp_().setSolventVolumeFlux(solventFlux); + } else { + asImp_().setSolventVolumeFlux(0.0); + } + } } public: @@ -461,7 +478,6 @@ public: if (!enableBoundaryMassFlux) return; - Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx); // estimate the gravity correction: for performance reasons we use a simplified @@ -520,15 +536,11 @@ public: const auto& up = intQuantsIn; // deal with water induced rock compaction - const double transMult = Toolbox::value(up.rockCompTransMultiplier()); + const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier()); transModified *= transMult; volumeFlux[phaseIdx] = pressureDifference[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); - - // TODO: Figure out if this did have any effect. It should? - // if (enableSolvent && phaseIdx == gasPhaseIdx) - // asImp_().setSolventVolumeFlux( pressureDifference[phaseIdx]*up.solventMobility()*(-transModified/faceArea)); } else { // compute the phase mobility using the material law parameters of the @@ -540,10 +552,6 @@ public: const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx); volumeFlux[phaseIdx] = pressureDifference[phaseIdx]*mob*(-transModified/faceArea); - - // Solvent inflow is not yet supported - // if (enableSolvent && phaseIdx == gasPhaseIdx) - // asImp_().setSolventVolumeFlux(0.0); } } } From 54284fad7f3711551c96f14c3e11408a37084b82 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 29 Sep 2022 14:13:49 +0200 Subject: [PATCH 7/7] Remove unneeded timeIdx argument. --- ebos/eclfluxmodule.hh | 2 -- 1 file changed, 2 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index ed0bcf9a2..51d50cb61 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -431,7 +431,6 @@ protected: globalSpaceIdx, intQuantsIn, scvfIdx, - timeIdx, faceArea, zEx, exFluidState, @@ -464,7 +463,6 @@ public: const unsigned globalSpaceIdx, const IntensiveQuantities& intQuantsIn, const unsigned bfIdx, - const unsigned timeIdx, const double faceArea, const double zEx, const FluidState& exFluidState,