From 2afa58f0d6a43ef1407be78d9969178b92beab1c Mon Sep 17 00:00:00 2001 From: hanchuan Date: Tue, 23 Aug 2022 01:21:11 +0200 Subject: [PATCH 1/6] [pnm] Regularization for conductivity-pc/sw relation in pore-network model --- dumux/flux/porenetwork/advection.hh | 76 +++++-- .../porenetwork/throat/transmissibility2p.hh | 191 ++++++++++++++++++ dumux/porenetwork/2p/fluxvariablescache.hh | 171 +++++++++++++--- test/porenetwork/2p/params.input | 3 + 4 files changed, 404 insertions(+), 37 deletions(-) diff --git a/dumux/flux/porenetwork/advection.hh b/dumux/flux/porenetwork/advection.hh index ae81bbcc5f..d0d9a88cbf 100644 --- a/dumux/flux/porenetwork/advection.hh +++ b/dumux/flux/porenetwork/advection.hh @@ -27,6 +27,7 @@ #include #include +#include namespace Dumux::PoreNetwork::Detail { @@ -125,42 +126,91 @@ public: const int wPhaseIdx = spatialParams.template wettingPhase(element, elemVolVars); const bool invaded = fluxVarsCache.invaded(); const Scalar pc = fluxVarsCache.pc(); - const Scalar pcEntry = fluxVarsCache.pcEntry(); - const Scalar pcSnapoff = fluxVarsCache.pcSnapoff(); + // 1p wetting phase transmissibility const Scalar Kw1p = Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx); + // regulariztaion interval for invasion event + const Scalar invasionLeft = fluxVarsCache.regInvasionInterval(0); + const Scalar invasionRight = fluxVarsCache.regInvasionInterval(1); + // regularization interval for snapoff event + const Scalar snapoffLeft = fluxVarsCache.regSnapoffInterval(0); + const Scalar snapoffRight = fluxVarsCache.regSnapoffInterval(1); if (phaseIdx == wPhaseIdx) { - if (!invaded) // not invaded in last time step, drainage curve sohuld be used + if (!invaded) // not invaded in last time step { - if (pc < pcEntry) + if ( pc < invasionLeft ) return Kw1p; - else // invasion happens - return Transmissibility::wettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); //2p + // the regularization interval is [pce, pce + reg] + else if ( pc < invasionRight ) + { + auto entryKw = Transmissibility::entryWettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + auto slopeEntry = Transmissibility::dKwdPcEntry(element, fvGeometry, scvf, fluxVarsCache); + // TODO: only the slope for circlular throat is 0.0! + const auto slopes = std::array{0.0, slopeEntry}; + auto optionalKnSpline_ = Spline(invasionLeft, invasionRight,// x0, x1 + Kw1p, entryKw, // y0, y1 + slopes[0], slopes[1]); // m0, m1 + return optionalKnSpline_.eval(pc); + } + else + return Transmissibility::wettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); } - else // invaded in last time step, imbibition curve should be used + else // invaded in last time step { - if (pc < pcSnapoff) // snapoff happens - return Kw1p; + // the regularization interval is [pcs - reg, pcs] + if (pc < snapoffLeft) + return Kw1p; // snapoff occurs + else if (pc < snapoffRight) // reg interval + { + auto snapoffKw = Transmissibility::snapoffWettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + auto slopeSnapoff = Transmissibility::dKwdPcSnapoff(element, fvGeometry, scvf, fluxVarsCache); + const auto slopes = std::array{0.0, slopeSnapoff}; + auto optionalKnSpline_ = Spline(snapoffLeft, snapoffRight , // x0, x1 + Kw1p, snapoffKw, // y0, y1 + slopes[0], slopes[1]); // m0, m1 + return optionalKnSpline_.eval(pc); + + } else return Transmissibility::wettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); //2p } } else // non-wetting phase { - if (!invaded) // last time step not invaded, drainage curve is used + if (!invaded) { // the regularization interval is [pce, pce + reg] - if (pc < pcEntry) + if (pc < invasionLeft) return 0.0; + else if (pc < invasionRight) + { + auto entryKn = Transmissibility::entryNonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + auto slopeEntry = Transmissibility::dKndPcEntry(element, fvGeometry, scvf, fluxVarsCache); + const auto slopes = std::array{0.0, slopeEntry}; + auto optionalKnSpline_ = Spline(invasionLeft, invasionRight, // x0, x1 + 0.0, entryKn, // y0, y1 + slopes[0], slopes[1]); // m0, m1 + return optionalKnSpline_.eval(pc); + } else return Transmissibility::nonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache); } - else // last time step is invaded, imbibition curve is used + else { // the regularzazion interval is [pcs - reg, pcs] - if (pc < pcSnapoff) + if (pc < snapoffLeft) return 0.0; + else if (pc < snapoffRight) + { + auto snapoffKn = Transmissibility::snapoffNonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + auto slopeSnapoff = Transmissibility::dKndPcSnapoff(element, fvGeometry, scvf, fluxVarsCache); + const auto slopes = std::array{slopeSnapoff, 0.0}; + auto optionalKnSpline_ = Spline(snapoffLeft, snapoffRight, // x0, x1 + 0.0, snapoffKn, // y0, y1 + slopes[1], slopes[0]); // m0, m1 + return optionalKnSpline_.eval(pc); + } else return Transmissibility::nonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache); } diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/throat/transmissibility2p.hh b/dumux/material/fluidmatrixinteractions/porenetwork/throat/transmissibility2p.hh index df931fa59c..9624b35539 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/throat/transmissibility2p.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/throat/transmissibility2p.hh @@ -119,6 +119,111 @@ struct RansohoffRadke return result; } + + /*! + * \brief Returns the integral conductivity of all wetting layers at critical status entry pressure + */ + template + static Scalar entryWettingLayerTransmissibility(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar rC = fluxVarsCache.curvatureRadiusInvasion(0); + const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); + const auto shape = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); + const auto numCorners = Throat::numCorners(shape); + + // treat the wetting film layer in each corner of the throat individually (might have different corner half-angle and beta) + Scalar result = 0.0; + for (int i = 0; i < numCorners; ++i) + result += fluxVarsCache.entryWettingLayerArea(i) * rC * rC / (throatLength * fluxVarsCache.wettingLayerFlowVariables().creviceResistanceFactor(i)); + + return result; + } + + /*! + * \brief Returns the integral conductivity of all wetting layers at critical status snapoff + */ + template + static Scalar snapoffWettingLayerTransmissibility(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar rC = fluxVarsCache.curvatureRadiusSnapoff(0); + + const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); + const auto shape = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); + const auto numCorners = Throat::numCorners(shape); + + // treat the wetting film layer in each corner of the throat individually (might have different corner half-angle and beta) + Scalar result = 0.0; + for (int i = 0; i < numCorners; ++i) + result += fluxVarsCache.snapoffWettingLayerArea(i) * rC * rC / (throatLength * fluxVarsCache.wettingLayerFlowVariables().creviceResistanceFactor(i)); + + return result; + } + + /*! + * \brief Returns derivative of wetting layer transmissibility at entry pressure + * For convinience now we calculate the derivative numerically + * TODO: Implement the analytical formular of this derivative + */ + template + static Scalar dKwdPcEntry(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + const Scalar Kw = entryWettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar rCdelta = fluxVarsCache.curvatureRadiusInvasion(1); + + const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); + const auto shape = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); + const auto numCorners = Throat::numCorners(shape); + + // treat the wetting film layer in each corner of the throat individually (might have different corner half-angle and beta) + Scalar deltaKw = 0.0; + for (int i = 0; i < numCorners; ++i) + deltaKw += fluxVarsCache.epsilonEntryWettingLayerArea(i) * rCdelta * rCdelta / (throatLength * fluxVarsCache.wettingLayerFlowVariables().creviceResistanceFactor(i)); + + const auto epsilonPc = fluxVarsCache.epsilonPc(); + auto result = (deltaKw - Kw)/epsilonPc; + return result; + } + + /*! + * \brief Returns derivative of wetting layer transmissibility at entry pressure + * For convinience now we calculate the derivative numerically + * TODO: Implement the analytical formular of this derivative + */ + template + static Scalar dKwdPcSnapoff(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + const Scalar Kw = snapoffWettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar rCdelta = fluxVarsCache.curvatureRadiusSnapoff(1); + + const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); + const auto shape = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); + const auto numCorners = Throat::numCorners(shape); + + // treat the wetting film layer in each corner of the throat individually (might have different corner half-angle and beta) + Scalar deltaKw = 0.0; + for (int i = 0; i < numCorners; ++i) + deltaKw += fluxVarsCache.epsilonSnapoffWettingLayerArea(i) * rCdelta * rCdelta / (throatLength * fluxVarsCache.wettingLayerFlowVariables().creviceResistanceFactor(i)); + + const auto epsilonPc = fluxVarsCache.epsilonPc(); + auto result = (deltaKw - Kw)/epsilonPc; + return result; + } }; } // end namespace WettingLayerTransmissibility @@ -197,6 +302,92 @@ struct BakkeOren const Scalar result = rEff*rEff*aNw / (8.0*throatLength); return result; } + + /*! + * \brief Returns the conductivity of a throat at entry pressure. + * + * See Bakke & Oren (1997), eq. 9 + */ + template + static Scalar entryNonWettingPhaseTransmissibility(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + // Tora et al. (2012), quite close for single-phase value of square + using std::sqrt; + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar aNwCrit = fluxVarsCache.regBoundaryNonwettingThroatAreaInvasion(0); + const Scalar rEff = 0.5 * (sqrt(aNwCrit / M_PI) + fluxVarsCache.throatInscribedRadius()); + const Scalar result = rEff * rEff * aNwCrit / (8.0 * throatLength); + return result; + } + + /*! + * \brief Returns the conductivity of a throat at snap-off pressure. + * + * See Bakke & Oren (1997), eq. 9 + */ + template + static Scalar snapoffNonWettingPhaseTransmissibility(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + // Tora et al. (2012), quite close for single-phase value of square + using std::sqrt; + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar aNwCrit = fluxVarsCache.regBoundaryNonWettingThroatAreaSnapoff(0); + const Scalar rEff = 0.5*(sqrt(aNwCrit / M_PI) + fluxVarsCache.throatInscribedRadius()); + const Scalar result = rEff*rEff*aNwCrit / (8.0*throatLength); + return result; + } + + /*! + * \brief Returns the conductivity of a throat at entry pressure. + * + * See Bakke & Oren (1997), eq. 9 + */ + template + static Scalar dKndPcEntry(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + // Tora et al. (2012), quite close for single-phase value of square + const Scalar Kn = entryNonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + using std::sqrt; + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar aNwCritDelta = fluxVarsCache.regBoundaryNonwettingThroatAreaInvasion(1); + const Scalar rEffDelta = 0.5*(sqrt(aNwCritDelta / M_PI) + fluxVarsCache.throatInscribedRadius()); + const Scalar KnDelta = rEffDelta*rEffDelta*aNwCritDelta / (8.0*throatLength); + const auto epsilonPc = fluxVarsCache.epsilonPc(); + auto result = (KnDelta - Kn)/epsilonPc; + return result; + } + + /*! + * \brief Returns the conductivity of a throat at snap-off pressure. + * + * See Bakke & Oren (1997), eq. 9 + */ + template + static Scalar dKndPcSnapoff(const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + // Tora et al. (2012), quite close for single-phase value of square + const Scalar Kn = snapoffNonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache); + using std::sqrt; + const Scalar throatLength = fluxVarsCache.throatLength(); + const Scalar aNwCritDelta = fluxVarsCache.regBoundaryNonWettingThroatAreaSnapoff(1); + const Scalar rEffDelta = 0.5*(sqrt(aNwCritDelta / M_PI) + fluxVarsCache.throatInscribedRadius()); + const Scalar KnDelta = rEffDelta*rEffDelta*aNwCritDelta / (8.0*throatLength); + const auto epsilonPc = fluxVarsCache.epsilonPc(); + auto result = (KnDelta - Kn)/epsilonPc; + return result; + } }; } // end namespace NonWettingPhaseTransmissibility } // end namespace Dumux::PoreNetwork diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index f51264a1a8..10f2edf701 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -45,6 +45,11 @@ class TwoPFluxVariablesCache public: + TwoPFluxVariablesCache() + { + regPercent_ = getParam("Regularization.RegPercentage", 0.01); + } + template void update(const Problem& problem, @@ -57,18 +62,46 @@ public: const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); throatCrossSectionShape_ = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); throatShapeFactor_ = fvGeometry.gridGeometry().throatShapeFactor(eIdx); - pc_ = std::max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure()); pcEntry_ = problem.spatialParams().pcEntry(element, elemVolVars); + pcSnapoff_ = problem.spatialParams().pcSnapoff(element, elemVolVars); + + // take the average surface tension of both adjacent pores TODO: is this correct? + surfaceTension_ = 0.5*(elemVolVars[0].surfaceTension() + elemVolVars[1].surfaceTension()); - const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); + pc_ = std::max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure()); + + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); const auto& elemSol = elementSolution(element, elemVolVars, fvGeometry); const auto& spatialParams = problem.spatialParams(); - auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, outsideScv, elemSol); + auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, insideScv, elemSol); // open a file to write the analytical pcEntry and Sw for comparison purposes std::ofstream file("analytical_SwEntry"); file << std::setprecision(15) << "pcEntry: " << pcEntry_ << " Sw: " << fluidMatrixInteraction.sw(pcEntry_) << std::endl; - pcSnapoff_ = problem.spatialParams().pcSnapoff(element, elemVolVars); + regInvasionInterval_[0] = pcEntry_; + regSnapoffInterval_[1] = pcSnapoff_; + +#if REGULARIZATIONWITHPRESSURE + regInvasionInterval_[1] = pcEntry_ * (1 + regPercent_); + regSnapoffInterval_[0] = pcSnapoff_* (1 - regPercent_); +#else + const auto swEntry = fluidMatrixInteraction.sw(pcEntry_); + const auto swRegEntry = swEntry - regPercent_; + regInvasionInterval_[1] = fluidMatrixInteraction.pc(swRegEntry); + + const auto swSnapoff = fluidMatrixInteraction.sw(pcSnapoff_); + const auto swRegSnapoff = swSnapoff + regPercent_; + regSnapoffInterval_[0] = fluidMatrixInteraction.pc(swRegSnapoff); +#endif + regInvasionInterval_[2] = regInvasionInterval_[1] + epsilonPc_; + regSnapoffInterval_[2] = pcSnapoff_ + epsilonPc_; + + curvatureRadiusInvasion_[0] = surfaceTension_/regInvasionInterval_[1]; + curvatureRadiusInvasion_[1] = surfaceTension_/regInvasionInterval_[2]; + + curvatureRadiusSnapoff_[0] = surfaceTension_/regSnapoffInterval_[1]; + curvatureRadiusSnapoff_[1] = surfaceTension_/regSnapoffInterval_[2]; + throatInscribedRadius_ = problem.spatialParams().throatInscribedRadius(element, elemVolVars); throatLength_ = problem.spatialParams().throatLength(element, elemVolVars); invaded_ = invaded; @@ -78,34 +111,37 @@ public: using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem; nPhaseIdx_ = 1 - spatialParams.template wettingPhase(element, elemVolVars); - // take the average surface tension of both adjacent pores TODO: is this correct? - surfaceTension_ = 0.5*(elemVolVars[0].surfaceTension() + elemVolVars[1].surfaceTension()); const auto& cornerHalfAngles = spatialParams.cornerHalfAngles(element); wettingLayerArea_.clear(); wettingLayerArea_.resize(cornerHalfAngles.size()); const Scalar totalThroatCrossSectionalArea = spatialParams.throatCrossSectionalArea(element, elemVolVars); + entryWettingLayerArea_.clear(); entryWettingLayerArea_.resize(cornerHalfAngles.size()); + snapoffWettingLayerArea_.clear(); snapoffWettingLayerArea_.resize(cornerHalfAngles.size()); + epsilonEntryWettingLayerArea_.clear(); epsilonEntryWettingLayerArea_.resize(cornerHalfAngles.size()); + epsilonSnapoffWettingLayerArea_.clear(); epsilonSnapoffWettingLayerArea_.resize(cornerHalfAngles.size()); - if (invaded) // two-phase flow + const Scalar theta = spatialParams.contactAngle(element, elemVolVars); + for (int i = 0; i< cornerHalfAngles.size(); ++i) { - const Scalar theta = spatialParams.contactAngle(element, elemVolVars); - for (int i = 0; i< cornerHalfAngles.size(); ++i) - wettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadius(), theta, cornerHalfAngles[i]); - - // make sure the wetting phase area does not exceed the total cross-section area - throatCrossSectionalArea_[wPhaseIdx()] = std::min( - std::accumulate(wettingLayerArea_.begin(), wettingLayerArea_.end(), 0.0), - totalThroatCrossSectionalArea - ); - throatCrossSectionalArea_[nPhaseIdx()] = totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]; + wettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadius(), theta, cornerHalfAngles[i]); + entryWettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadiusInvasion(0), theta, cornerHalfAngles[i]); + snapoffWettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadiusSnapoff(0), theta, cornerHalfAngles[i]); + epsilonEntryWettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadiusInvasion(1), theta, cornerHalfAngles[i]); + epsilonSnapoffWettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadiusSnapoff(1), theta, cornerHalfAngles[i]); } - else // single-phase flow - { - for (int i = 0; i< cornerHalfAngles.size(); ++i) - wettingLayerArea_[i] = 0.0; - throatCrossSectionalArea_[wPhaseIdx()] = totalThroatCrossSectionalArea; - throatCrossSectionalArea_[nPhaseIdx()] = 0.0; - } + // make sure the wetting phase area does not exceed the total cross-section area + throatCrossSectionalArea_[wPhaseIdx()] = std::min( + std::accumulate(wettingLayerArea_.begin(), wettingLayerArea_.end(), 0.0), + totalThroatCrossSectionalArea + ); + throatCrossSectionalArea_[nPhaseIdx()] = totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]; + + regBoundaryWettingThroatAreaInvasion_[0] = std::min(std::accumulate(entryWettingLayerArea_.begin(), entryWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); + regBoundaryWettingThroatAreaInvasion_[1] = std::min(std::accumulate(epsilonEntryWettingLayerArea_.begin(), epsilonEntryWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); + regBoundaryNonwettingThroatAreaInvasion_[0] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[0]; + regBoundaryNonwettingThroatAreaInvasion_[1] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[1]; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { @@ -243,6 +279,80 @@ public: Scalar poreToPoreDistance() const { return poreToPoreDistance_; } + /*! + * \brief Returns the regularization Interval for invasion event + */ + Scalar regInvasionInterval(const int Idx) const + { return regInvasionInterval_[Idx]; } + + /*! + * \brief Returns the regularization Interval for snap-off event + */ + Scalar regSnapoffInterval(const int Idx) const + { return regSnapoffInterval_[Idx]; } + + /*! + * \brief Returns the curvature radius within the throat corresponds to entry pressure. + */ + Scalar curvatureRadiusInvasion(const int Idx) const + { return curvatureRadiusInvasion_[Idx];} + + /*! + * \brief Returns the curvature radius within the throat corresponds to snap off pressure. + */ + Scalar curvatureRadiusSnapoff(const int Idx) const + { return curvatureRadiusSnapoff_[Idx];} + + /*! + * \brief Returns delta pc used to calculate numerical derivative of K + */ + Scalar epsilonPc() const + { return epsilonPc_; } + + /*! + * \brief Returns wetting layer area for each corner + * at regularization pc right interval boundary for invasion + */ + Scalar entryWettingLayerArea(const int cornerIdx) const + { return entryWettingLayerArea_[cornerIdx]; } + + Scalar epsilonEntryWettingLayerArea(const int cornerIdx) const + { return epsilonEntryWettingLayerArea_[cornerIdx]; } + + /*! + * \brief Returns wetting layer area for each corner + * at regularization pc right interval boundary for snap-off + */ + Scalar snapoffWettingLayerArea(const int cornerIdx) const + { return snapoffWettingLayerArea_[cornerIdx]; } + + Scalar epsilonSnapoffWettingLayerArea(const int cornerIdx) const + { return epsilonSnapoffWettingLayerArea_[cornerIdx]; } + + /*! + * \brief Returns total wetting layer area for invasion + */ + Scalar regBoundaryWettingThroatAreaInvasion(const int Idx) const + { return regBoundaryWettingThroatAreaInvasion_[Idx]; } + + /*! + * \brief Returns total wetting layer area for snap-off + */ + Scalar regBoundaryWettingThroatAreaSnapoff(const int Idx) const + { return regBoundaryWettingThroatAreaSnapoff_[Idx]; } + + /*! + * \brief Returns total non-wetting throat area for invasion + */ + Scalar regBoundaryNonwettingThroatAreaInvasion(const int Idx) const + { return regBoundaryNonwettingThroatAreaInvasion_[Idx]; } + + /*! + * \brief Returns total non-wetting throat area for snap-off + */ + Scalar regBoundaryNonWettingThroatAreaSnapoff(const int Idx) const + { return regBoundaryNonwettingThroatAreaSnapoff_[Idx]; } + private: Throat::Shape throatCrossSectionShape_; Scalar throatShapeFactor_; @@ -262,6 +372,19 @@ private: typename AdvectionType::Transmissibility::SinglePhaseCache singlePhaseCache_; typename AdvectionType::Transmissibility::NonWettingPhaseCache nonWettingPhaseCache_; typename AdvectionType::Transmissibility::WettingLayerCache wettingLayerCache_; + + static constexpr double epsilonPc_ = 1e-8; + Scalar regPercent_; + std::array regInvasionInterval_; + std::array regSnapoffInterval_; + std::array curvatureRadiusInvasion_; + std::array curvatureRadiusSnapoff_; + NumCornerVector entryWettingLayerArea_, epsilonEntryWettingLayerArea_; + NumCornerVector snapoffWettingLayerArea_, epsilonSnapoffWettingLayerArea_; + std::array regBoundaryWettingThroatAreaInvasion_; + std::array regBoundaryWettingThroatAreaSnapoff_; + std::array regBoundaryNonwettingThroatAreaInvasion_; + std::array regBoundaryNonwettingThroatAreaSnapoff_; }; } // end Dumux::PoreNetwork diff --git a/test/porenetwork/2p/params.input b/test/porenetwork/2p/params.input index a798ce23cc..9448ddb867 100644 --- a/test/porenetwork/2p/params.input +++ b/test/porenetwork/2p/params.input @@ -28,3 +28,6 @@ Verbosity = true [SpatialParams] HighSwRegularizationMethod = Spline + +[Regularization] +RegPercentage = 0.1 \ No newline at end of file -- GitLab From ddd39b0664c1dc4bb3bf4492f875c6976f90c156 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Wed, 24 Aug 2022 10:21:24 +0200 Subject: [PATCH 2/6] [pnm] Do not block nonwetting flow out --- dumux/porenetwork/2p/invasionstate.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/porenetwork/2p/invasionstate.hh b/dumux/porenetwork/2p/invasionstate.hh index 7f18f06bbd..6322574eeb 100644 --- a/dumux/porenetwork/2p/invasionstate.hh +++ b/dumux/porenetwork/2p/invasionstate.hh @@ -207,7 +207,7 @@ private: }; // Block non-wetting phase flux out of the outlet - static const auto blockNonwettingPhase = getParamFromGroup>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector{Labels::outlet}); + static const auto blockNonwettingPhase = getParamFromGroup>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector{}); if (!blockNonwettingPhase.empty() && std::find(blockNonwettingPhase.begin(), blockNonwettingPhase.end(), gridGeometry.throatLabel(eIdx)) != blockNonwettingPhase.end()) { invadedCurrentTimeStep_[eIdx] = false; -- GitLab From b7f3030e679ac5093ac1d4c51e879218505d2306 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Wed, 24 Aug 2022 10:24:18 +0200 Subject: [PATCH 3/6] [pnm-reg] Fix CMakeLists.txt --- test/porenetwork/2p/CMakeLists.txt | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/test/porenetwork/2p/CMakeLists.txt b/test/porenetwork/2p/CMakeLists.txt index ad938d1a9b..5c4125bc13 100644 --- a/test/porenetwork/2p/CMakeLists.txt +++ b/test/porenetwork/2p/CMakeLists.txt @@ -5,21 +5,13 @@ add_input_file_links() dumux_add_test(NAME test_pnm_2p SOURCES main.cc LABELS porenetwork - COMPILE_DEFINITIONS ISOTHERMAL=1 + COMPILE_DEFINITIONS ISOTHERMAL=1 REGULARIZATIONWITHPRESSURE=1 COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )" - CMD_ARGS --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_2p-reference.vtp - ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2p-00107.vtp - --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2p") + CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )") dumux_add_test(NAME test_pnm_2pni SOURCES main.cc LABELS porenetwork - COMPILE_DEFINITIONS ISOTHERMAL=0 + COMPILE_DEFINITIONS ISOTHERMAL=0 REGULARIZATIONWITHPRESSURE=1 COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )" - CMD_ARGS --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_2pni-reference.vtp - ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pni-00133.vtp - --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pni params_ni.input -Problem.Name test_pnm_2pni") + CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )") -- GitLab From a3d43aefff18714386a4e480c0f97795f6871fc4 Mon Sep 17 00:00:00 2001 From: Maziar Veyskarami Date: Wed, 24 Aug 2022 12:25:32 +0200 Subject: [PATCH 4/6] [pnm-regularization][2p-test] increase the time limit of the test in the CMakeLists --- test/porenetwork/2p/CMakeLists.txt | 8 ++++++-- test/porenetwork/2p/params.input | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/test/porenetwork/2p/CMakeLists.txt b/test/porenetwork/2p/CMakeLists.txt index 5c4125bc13..c02a388b1c 100644 --- a/test/porenetwork/2p/CMakeLists.txt +++ b/test/porenetwork/2p/CMakeLists.txt @@ -6,12 +6,16 @@ dumux_add_test(NAME test_pnm_2p SOURCES main.cc LABELS porenetwork COMPILE_DEFINITIONS ISOTHERMAL=1 REGULARIZATIONWITHPRESSURE=1 + CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )" + TIMEOUT "1000" COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )") + CMD_ARGS --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2p params.input") dumux_add_test(NAME test_pnm_2pni SOURCES main.cc LABELS porenetwork COMPILE_DEFINITIONS ISOTHERMAL=0 REGULARIZATIONWITHPRESSURE=1 + CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )" + TIMEOUT "1000" COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )") + CMD_ARGS --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pni params_ni.input") diff --git a/test/porenetwork/2p/params.input b/test/porenetwork/2p/params.input index 9448ddb867..c406b2fac9 100644 --- a/test/porenetwork/2p/params.input +++ b/test/porenetwork/2p/params.input @@ -30,4 +30,4 @@ Verbosity = true HighSwRegularizationMethod = Spline [Regularization] -RegPercentage = 0.1 \ No newline at end of file +RegPercentage = 0.1 -- GitLab From 93b8b50e1de63a6bf6a8755bd37b82740433e563 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Thu, 25 Aug 2022 10:41:53 +0200 Subject: [PATCH 5/6] [pnm] Prevent negative non-wetting area --- dumux/porenetwork/2p/fluxvariablescache.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index 10f2edf701..8e35b4aa09 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -135,12 +135,12 @@ public: std::accumulate(wettingLayerArea_.begin(), wettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea ); - throatCrossSectionalArea_[nPhaseIdx()] = totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]; + throatCrossSectionalArea_[nPhaseIdx()] = std::max( (totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]), 0.0); regBoundaryWettingThroatAreaInvasion_[0] = std::min(std::accumulate(entryWettingLayerArea_.begin(), entryWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); regBoundaryWettingThroatAreaInvasion_[1] = std::min(std::accumulate(epsilonEntryWettingLayerArea_.begin(), epsilonEntryWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); - regBoundaryNonwettingThroatAreaInvasion_[0] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[0]; - regBoundaryNonwettingThroatAreaInvasion_[1] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[1]; + regBoundaryNonwettingThroatAreaInvasion_[0] = std::max( (totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[0]), 0.0); + regBoundaryNonwettingThroatAreaInvasion_[1] = std::max( (totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[1]), 0.0); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) -- GitLab From 1edd62623842ecc46f833f37aa7aa88f99e8dd93 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Thu, 25 Aug 2022 11:23:20 +0200 Subject: [PATCH 6/6] [pnm][bug-fix] Initialize snapoff-area --- dumux/porenetwork/2p/fluxvariablescache.hh | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index 8e35b4aa09..3b1f4aef19 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -135,12 +135,17 @@ public: std::accumulate(wettingLayerArea_.begin(), wettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea ); - throatCrossSectionalArea_[nPhaseIdx()] = std::max( (totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]), 0.0); + throatCrossSectionalArea_[nPhaseIdx()] = totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]; regBoundaryWettingThroatAreaInvasion_[0] = std::min(std::accumulate(entryWettingLayerArea_.begin(), entryWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); regBoundaryWettingThroatAreaInvasion_[1] = std::min(std::accumulate(epsilonEntryWettingLayerArea_.begin(), epsilonEntryWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); - regBoundaryNonwettingThroatAreaInvasion_[0] = std::max( (totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[0]), 0.0); - regBoundaryNonwettingThroatAreaInvasion_[1] = std::max( (totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[1]), 0.0); + regBoundaryNonwettingThroatAreaInvasion_[0] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[0]; + regBoundaryNonwettingThroatAreaInvasion_[1] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaInvasion_[1]; + + regBoundaryWettingThroatAreaSnapoff_[0] = std::min(std::accumulate(snapoffWettingLayerArea_.begin(), snapoffWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); + regBoundaryWettingThroatAreaSnapoff_[1] = std::min(std::accumulate(epsilonSnapoffWettingLayerArea_.begin(), epsilonSnapoffWettingLayerArea_.end(), 0.0), totalThroatCrossSectionalArea); + regBoundaryNonwettingThroatAreaSnapoff_[0] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaSnapoff_[0]; + regBoundaryNonwettingThroatAreaSnapoff_[1] = totalThroatCrossSectionalArea - regBoundaryWettingThroatAreaSnapoff_[1]; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) -- GitLab