diff --git a/dumux/flux/porenetwork/fourierslaw.hh b/dumux/flux/porenetwork/fourierslaw.hh
index 6db79c10cc4b2405b5f1880f1364ed5ae54ed418..2f091f6a20db667fce4a24f20ba5c3b608aed6cd 100644
--- a/dumux/flux/porenetwork/fourierslaw.hh
+++ b/dumux/flux/porenetwork/fourierslaw.hh
@@ -26,13 +26,22 @@
#define DUMUX_FLUX_PNM_FOURIERS_LAW_HH
#include <dumux/common/math.hh>
+#include <dumux/flux/referencesystemformulation.hh>
+#include <type_traits>
namespace Dumux::PoreNetwork {
+namespace Detail {
+
+struct NoDiffusionType {};
+
+} // end namespace Detail
+
/*!
* \ingroup PoreNetworkModels
* \brief Specialization of Fourier's Law for the pore-network model.
*/
+template<class MolecularDiffusionType = Detail::NoDiffusionType>
struct PNMFouriersLaw
{
@@ -68,7 +77,41 @@ struct PNMFouriersLaw
const Scalar gradT = deltaT/fluxVarsCache.throatLength();
heatflux += thermalConductivity*gradT*area;
+
+ if constexpr (!std::is_same_v<MolecularDiffusionType, Detail::NoDiffusionType>)
+ heatflux += componentEnthalpyFlux_(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache, phaseIdx);
+ }
+
+ return heatflux;
+ }
+
+private:
+ template<class Problem, class Element, class FVElementGeometry,
+ class ElementVolumeVariables, class ElementFluxVariablesCache>
+ static auto componentEnthalpyFlux_(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const int phaseIdx)
+ {
+ using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
+ Scalar heatflux = 0.0;
+ using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem;
+ const auto diffusiveFlux = MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
+ for (int compIdx = 0; compIdx < ElementVolumeVariables::VolumeVariables::numFluidComponents(); ++compIdx)
+ {
+ const bool insideIsUpstream = diffusiveFlux[compIdx] > 0.0;
+ const auto& upstreamVolVars = insideIsUpstream ? elemVolVars[scvf.insideScvIdx()] : elemVolVars[scvf.outsideScvIdx()];
+ const Scalar componentEnthalpy = FluidSystem::componentEnthalpy(upstreamVolVars.fluidState(), phaseIdx, compIdx);
+
+ if (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+ heatflux += diffusiveFlux[compIdx] * componentEnthalpy;
+ else
+ heatflux += diffusiveFlux[compIdx] * FluidSystem::molarMass(compIdx) * componentEnthalpy;
}
+
return heatflux;
}
};
diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh
index e13600e8b6969f5ef85c130d609107e4fb78c0c4..c97b28fcc132065e4704ec0b226b47a7aa272033 100644
--- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh
+++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh
@@ -57,10 +57,8 @@ struct PlatonicBodyParams
const ElemSol& elemSol)
: shape_(spatialParams.gridGeometry().poreGeometry()[scv.dofIndex()])
, radius_(spatialParams.poreInscribedRadius(element, scv, elemSol))
- {
- static const Scalar surfaceTension = getParam<Scalar>("SpatialParams.SurfaceTension", 0.0725); // TODO
- surfaceTension_ = surfaceTension;
- }
+ , surfaceTension_(spatialParams.surfaceTension(element, scv, elemSol))
+ {}
template<class SpatialParams, class Element, class SubControlVolume, class ElemSol>
void update(const SpatialParams& spatialParams,
@@ -71,9 +69,7 @@ struct PlatonicBodyParams
const auto& gridGeometry = spatialParams.gridGeometry();
shape_ = gridGeometry.poreGeometry()[scv.dofIndex()];
radius_ = spatialParams.poreInscribedRadius(element, scv, elemSol);
-
- static const Scalar surfaceTension = getParam<Scalar>("SpatialParams.SurfaceTension", 0.0725); // TODO
- surfaceTension_ = surfaceTension;
+ surfaceTension_ = spatialParams.surfaceTension(element, scv, elemSol);
}
Pore::Shape poreShape() const { return shape_; }
@@ -321,11 +317,6 @@ private:
if (sw < pcLowSw_)
return pcLowSwPcValue_() + pcDerivativeLowSw_() * (sw - pcLowSw_);
- auto linearCurveForHighSw = [&]()
- {
- return pcDerivativeHighSwEnd_()*(sw - 1.0);
- };
-
if (sw <= pcHighSw_)
return {}; // standard
else if (sw < 1.0) // regularized part below sw = 1.0
@@ -335,7 +326,10 @@ private:
return pcHighSwPcValue_() * pow(((1.0-sw)/(1.0-pcHighSw_)), 1.0/3.0);
else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear)
- return linearCurveForHighSw();
+ {
+ const Scalar slope = -pcHighSwPcValue_() / (1.0 - pcHighSw_);
+ return pcHighSwPcValue_() + (sw - pcHighSw_) * slope;
+ }
else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::spline)
return pcSpline_().eval(sw);
@@ -344,7 +338,7 @@ private:
DUNE_THROW(Dune::NotImplemented, "Regularization not method not implemented");
}
else // regularized part above sw = 1.0
- return linearCurveForHighSw();
+ return pcDerivativeHighSwEnd_()*(sw - 1.0);
}
/*!
@@ -516,6 +510,9 @@ private:
baseLawParamsPtr_ = &m->basicParams();
+ static const auto highSwFixedSlopeInput = getParamFromGroup<Scalar>(paramGroup, "RegularizationHighSwFixedSlope", std::numeric_limits<Scalar>::quiet_NaN());
+ highSwFixedSlope_ = highSwFixedSlopeInput;
+
// Note: we do not pre-calculate all end-point values (e.g., pcLowSwPcValue_ and pcDerivativeLowSw_)
// as done in, e.g., VanGenuchten. This is because this law will generally be instantiated only as
// a temporary object since all pore bodies generally have different parameters.
@@ -550,7 +547,11 @@ private:
// dpc_dsw at Sw = 1.0
Scalar pcDerivativeHighSwEnd_() const
{
- return (0.0 - pcHighSwPcValue_()) / (1.0 - pcHighSw_);
+ using std::isnan;
+ if (!isnan(highSwFixedSlope_))
+ return highSwFixedSlope_;
+ else
+ return (0.0 - pcHighSwPcValue_()) / (1.0 - pcHighSw_);
}
const Spline<Scalar>& pcSpline_() const
@@ -574,6 +575,7 @@ private:
const BaseLawParams* baseLawParamsPtr_;
mutable std::optional<Spline<Scalar>> optionalPcSpline_;
bool highSwSplineZeroSlope_;
+ Scalar highSwFixedSlope_;
};
/*!
diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh
index e248654efbd0c1dc622956f53e5c3b0118bf0240..5f3bab10a793a4fb85b8f5de947f634259a9d2f6 100644
--- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh
+++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh
@@ -34,11 +34,11 @@ namespace Dumux::PoreNetwork::FluidMatrix {
template<class ScalarT>
struct LocalRulesTraits
{
- using Tetrahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::tetrahedron>;
- using Cube = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::cube>;
- using Octahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::octahedron>;
- using Icosahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::icosahedron>;
- using Dodecahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::dodecahedron>;
+ using Tetrahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::tetrahedron, ScalarT>;
+ using Cube = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::cube, ScalarT>;
+ using Octahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::octahedron, ScalarT>;
+ using Icosahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::icosahedron, ScalarT>;
+ using Dodecahedron = TwoPLocalRulesPlatonicBodyDefault<Pore::Shape::dodecahedron, ScalarT>;
};
template<class ScalarT>
@@ -102,6 +102,7 @@ public:
{ return 2; }
template<class SpatialParams, class Element, class SubControlVolume, class ElemSol>
+ [[deprecated("Use constructor of BasicParams directly. Will be removed before release 3.4")]]
static BasicParams makeParams(const SpatialParams& spatialParams,
const Element& element,
const SubControlVolume& scv,
diff --git a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh
index 415bca66c6535c20789234566f3bb0a390720a7c..8ec48f9de5649117a959cabd70ec2f3382672d1f 100644
--- a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh
+++ b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh
@@ -157,9 +157,23 @@ public:
}
/*!
- * \brief Returns the parameter object for the Brooks-Corey material law.
+ * \brief Returns the surface tension \f$ N/m\f$
*
- * In this test, we use element-wise distributed material parameters.
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ */
+ template<class ElementSolution>
+ Scalar surfaceTension(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ static const Scalar gamma = getParam<Scalar>("SpatialParams.SurfaceTension", 0.0725); // default to surface tension of water/air
+ return gamma;
+ }
+
+ /*!
+ * \brief Returns the parameter object for the pore-local pc-Sw law
*
* \param element The current element
* \param scv The sub-control volume inside the element.
@@ -171,7 +185,7 @@ public:
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
- const auto params = LocalRules::makeParams(*this, element, scv, elemSol);
+ const auto params = typename LocalRules::BasicParams(*this, element, scv, elemSol);
return makeFluidMatrixInteraction(LocalRules(params, "SpatialParams"));
}
diff --git a/dumux/porenetwork/1pnc/model.hh b/dumux/porenetwork/1pnc/model.hh
index 807f6e20b61a958a5a3f10d298c7f202595051c2..7619f8c2dc1322d65332042a1cbf595d98841659 100644
--- a/dumux/porenetwork/1pnc/model.hh
+++ b/dumux/porenetwork/1pnc/model.hh
@@ -242,5 +242,12 @@ struct ThermalConductivityModel<TypeTag, TTag::PNMOnePNCNI>
using type = ThermalConductivityAverage<GetPropType<TypeTag, Properties::Scalar>>;
}; //!< Use the average for effective conductivities
+// template<class TypeTag>
+// struct HeatConductionType<TypeTag, TTag::PNMOnePNCNI>
+// {
+// TODO uncomment this as soon as there is a generalized approach for component enthalpies in all fluid systems
+// using type = Dumux::PoreNetwork::PNMFouriersLaw<GetPropType<TypeTag, MolecularDiffusionType>>;
+// }; //!< Use Fourier's law and also consider enthalpy transport by molecular diffusion
+
} // end namespace Dumux::Properties
#endif
diff --git a/dumux/porenetwork/2p/gridfluxvariablescache.hh b/dumux/porenetwork/2p/gridfluxvariablescache.hh
index af5d4afa8805c910009d0b34556fd7302026a36b..a6a1969aa39ec9a7224baa0851f62ac3807bec44 100644
--- a/dumux/porenetwork/2p/gridfluxvariablescache.hh
+++ b/dumux/porenetwork/2p/gridfluxvariablescache.hh
@@ -36,13 +36,12 @@ namespace Dumux::PoreNetwork {
* \ingroup PNMTwoPModel
* \brief Flux variable caches traits
*/
-template<class P, class FVC, class I, class L>
+template<class P, class FVC, class IS = TwoPInvasionState<P>>
struct PNMTwoPDefaultGridFVCTraits
{
using Problem = P;
using FluxVariablesCache = FVC;
- using Indices = I;
- using Labels = L;
+ using InvasionState = IS;
template<class GridFluxVariablesCache, bool cachingEnabled>
using LocalView = PNMTwoPElementFluxVariablesCache<GridFluxVariablesCache, cachingEnabled>;
@@ -69,9 +68,9 @@ class PNMTwoPGridFluxVariablesCache<P, FVC, true, Traits>
{
using Problem = typename Traits::Problem;
using ThisType = PNMTwoPGridFluxVariablesCache<P, FVC, true, Traits>;
- using InvasionState = TwoPInvasionState<Problem>;
+ using InvasionState = typename Traits::InvasionState;
- public:
+public:
//! export the flux variable cache type
using FluxVariablesCache = typename Traits::FluxVariablesCache;
@@ -140,7 +139,7 @@ class PNMTwoPGridFluxVariablesCache<P, FVC, false, Traits>
{
using Problem = typename Traits::Problem;
using ThisType = PNMTwoPGridFluxVariablesCache<P, FVC, false, Traits>;
- using InvasionState = TwoPInvasionState<Problem>;
+ using InvasionState = typename Traits::InvasionState;
public:
//! export the flux variable cache type
diff --git a/dumux/porenetwork/2p/invasionstate.hh b/dumux/porenetwork/2p/invasionstate.hh
index e6fce845b847e3f6482120d1682d872156fba694..8100a46bb500d937ab2e32423a4309173bbca493 100644
--- a/dumux/porenetwork/2p/invasionstate.hh
+++ b/dumux/porenetwork/2p/invasionstate.hh
@@ -19,42 +19,37 @@
/*!
* \file
* \ingroup PNMTwoPModel
- * \brief Global flux variable cache
+ * \brief Invasion state class for the two-phase PNM.
*/
#ifndef DUMUX_PNM_2P_INVASIONSTATE_HH
#define DUMUX_PNM_2P_INVASIONSTATE_HH
#include <vector>
#include <type_traits>
+#include <dune/common/std/type_traits.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/porenetwork/common/labels.hh>
namespace Dumux::PoreNetwork {
-#ifndef DOXYGEN
-namespace Detail {
-// helper struct detecting if the user-defined problem params class has a globalCapillaryPressure function
-struct hasGlobalCapillaryPressure
-{
- template<class Problem>
- auto operator()(const Problem& p)
- -> decltype(p.globalCapillaryPressure())
- {}
-};
-} // end namespace Detail
-#endif
-
/*!
* \ingroup PNMTwoPModel
- * \brief The grid flux variables cache for the two-phase PNM hodling the invasion state of the throats
- * \note The flux caches of the gridview are stored which is memory intensive but faster
+ * \brief This class updates the invasion state for the two-phase PNM.
*/
template<class P>
class TwoPInvasionState
{
using Problem = P;
+ template <class T>
+ using GlobalCapillaryPressureDetector = decltype(std::declval<T>().globalCapillaryPressure());
+
+ template<class T>
+ static constexpr bool hasGlobalCapillaryPressure()
+ { return Dune::Std::is_detected<GlobalCapillaryPressureDetector, T>::value; }
+
+ enum class EventType {invasion, snapOff, none};
+
public:
TwoPInvasionState(const Problem& problem) : problem_(problem)
@@ -71,10 +66,10 @@ public:
}
numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true);
- verbose_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.InvasionStateVerbosity", true);
+ verbose_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.Verbosity", true);
restrictToGlobalCapillaryPressure_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.RestrictInvasionToGlobalCapillaryPressure", false);
- if (decltype(isValid(Detail::hasGlobalCapillaryPressure())(problem)){})
+ if constexpr (hasGlobalCapillaryPressure<Problem>())
{
if (restrictToGlobalCapillaryPressure_)
std::cout << "\n *** Invasion behavior is restricted by a global capillary pressure defined in the problem! *** \n" << std::endl;
@@ -115,11 +110,14 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
// checks if invasion or snap-off occured after Newton iteration step
- if (invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf]))
+ if (const auto invasionResult = invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf]); invasionResult)
{
hasChangedInCurrentIteration_ = true;
- static constexpr auto cachingEnabled = std::integral_constant<bool, GridFluxVariablesCache::cachingEnabled>{};
- updateChangedFluxVarCache_(element, fvGeometry, gridFluxVarsCache, elemVolVars, scvf, cachingEnabled);
+ if constexpr (GridFluxVariablesCache::cachingEnabled)
+ {
+ const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
+ gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invadedCurrentIteration_[eIdx]);
+ }
}
}
}
@@ -194,36 +192,46 @@ private:
//! The switch for determining the invasion state of a pore throat. Called at the end of each Newton step.
template<class Element, class ElementVolumeVariables, class FluxVariablesCache>
- bool invasionSwitch_(const Element& element,
+ auto invasionSwitch_(const Element& element,
const ElementVolumeVariables& elemVolVars,
const FluxVariablesCache& fluxVarsCache)
{
+ using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
const auto& gridGeometry = problem_.gridGeometry();
const auto& spatialParams = problem_.spatialParams();
const auto eIdx = gridGeometry.elementMapper().index(element);
bool invadedBeforeSwitch = invadedCurrentIteration_[eIdx];
bool invadedAfterSwitch = invadedBeforeSwitch;
- const auto wPhaseIdx = spatialParams.template wettingPhase<typename ElementVolumeVariables::VolumeVariables::FluidSystem>(element, elemVolVars);
+ // Result type, containing the local scv index of the pore from which the invasion/snap-off occurred
+ // Evaluates to 'false' if no invasion/snap-off occurred
+ struct Result
+ {
+ std::uint8_t localScvIdxWithCriticalPc;
+ Scalar criticalPc;
+ EventType event = EventType::none;
+
+ operator bool() const
+ { return event != EventType::none; }
+ };
// Block non-wetting phase flux out of the outlet
- static const bool blockOutlet = getParamFromGroup<bool>(problem_.paramGroup(), "InvasionState.BlockOutletForNonwettingPhase", true);
- if (blockOutlet && gridGeometry.throatLabel(eIdx) == Labels::outlet)
+ static const auto blockNonwettingPhase = getParamFromGroup<std::vector<int>>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector<int>{Labels::outlet});
+ if (!blockNonwettingPhase.empty() && std::find(blockNonwettingPhase.begin(), blockNonwettingPhase.end(), gridGeometry.throatLabel(eIdx)) != blockNonwettingPhase.end())
{
invadedCurrentIteration_[eIdx] = false;
- return false;
+ return Result{}; // nothing happened
}
//Determine whether throat gets invaded or snap-off occurs
- using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
const std::array<Scalar, 2> pc = { elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure() };
const auto pcMax = std::max_element(pc.begin(), pc.end());
const Scalar pcEntry = fluxVarsCache.pcEntry();
const Scalar pcSnapoff = fluxVarsCache.pcSnapoff();
// check if there is a user-specified global capillary pressure which needs to be obeyed
- if (maybeRestrictToGlobalCapillaryPressure_(pcEntry, decltype(isValid(Detail::hasGlobalCapillaryPressure())(problem_)){}))
+ if (maybeRestrictToGlobalCapillaryPressure_(pcEntry))
{
if (*pcMax > pcEntry)
{
@@ -232,7 +240,7 @@ private:
}
invadedCurrentIteration_[eIdx] = false;
- return false;
+ return Result{}; //nothing happened
}
if (*pcMax > pcEntry)
@@ -242,68 +250,51 @@ private:
invadedCurrentIteration_[eIdx] = invadedAfterSwitch;
- if (invadedAfterSwitch != invadedBeforeSwitch && verbose_)
+ if (invadedBeforeSwitch == invadedAfterSwitch)
+ return Result{}; // nothing happened
+ else
{
- const std::array<Scalar, 2> sw = { elemVolVars[0].saturation(wPhaseIdx), elemVolVars[1].saturation(wPhaseIdx) };
- const auto scvIdx = pcMax - pc.begin();
- const auto vIdx = gridGeometry.gridView().indexSet().subIndex(element, scvIdx, 1);
- if (!invadedBeforeSwitch && invadedAfterSwitch)
- {
- std::cout << "Throat " << eIdx << " was invaded from pore " << vIdx << " :";
- std::cout << " pc: " << *pcMax;
- std::cout << ", pcEntry: " << spatialParams.pcEntry(element, elemVolVars);
- std::cout << ", sw: " << sw[scvIdx] << std::endl;
- }
- else if (invadedBeforeSwitch && !invadedAfterSwitch)
- {
- std::cout << "Snap-off occured at: " << eIdx << " from pore " << vIdx << " :";
- std::cout << " pc: " << *pcMax;
- std::cout << ", pcSnapoff: " << spatialParams.pcSnapoff(element, elemVolVars);
- std::cout << ", sw: " << sw[scvIdx] << std::endl;
- }
- else
- DUNE_THROW(Dune::InvalidStateException, "Invalid Process ");
- }
- return invadedBeforeSwitch != invadedAfterSwitch;
- }
+ Result result;
+ result.localScvIdxWithCriticalPc = std::distance(pc.begin(), pcMax);
+ result.criticalPc = *pcMax;
+ result.event = !invadedBeforeSwitch && invadedAfterSwitch ? EventType::invasion : EventType::snapOff;
- //! Update the fluxVarsCache after an invasion/snap-off occured (only for caching enabled)
- template<class Element, class FVElementGeometry, class GridFluxVariablesCache,
- class ElementVolumeVariables, class SubControlVolumeFace>
- void updateChangedFluxVarCache_(const Element& element,
- const FVElementGeometry& fvGeometry,
- GridFluxVariablesCache& gridFluxVarsCache,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- std::true_type)
- {
- const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
- gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invaded(element));
- }
+ if (verbose_)
+ {
+ const auto wPhaseIdx = spatialParams.template wettingPhase<typename ElementVolumeVariables::VolumeVariables::FluidSystem>(element, elemVolVars);
+ const std::array sw = { elemVolVars[0].saturation(wPhaseIdx), elemVolVars[1].saturation(wPhaseIdx) };
+ const auto vIdx = gridGeometry.gridView().indexSet().subIndex(element, result.localScvIdxWithCriticalPc, 1);
+ if (result.event == EventType::invasion)
+ {
+ std::cout << "Throat " << eIdx << " was invaded from pore " << vIdx << " :";
+ std::cout << " pc: " << *pcMax;
+ std::cout << ", pcEntry: " << spatialParams.pcEntry(element, elemVolVars);
+ std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl;
+ }
+ else
+ {
+ std::cout << "Snap-off occured at throat " << eIdx << " from pore " << vIdx << " :";
+ std::cout << " pc: " << *pcMax;
+ std::cout << ", pcSnapoff: " << spatialParams.pcSnapoff(element, elemVolVars);
+ std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl;
+ }
+ }
- //! Do nothing if the fluxVarCaches are not cached globally
- template<class Element, class FVElementGeometry, class GridFluxVariablesCache,
- class ElementVolumeVariables, class SubControlVolumeFace>
- void updateChangedFluxVarCache_(const Element& element,
- const FVElementGeometry& fvGeometry,
- GridFluxVariablesCache& gridFluxVarsCache,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- std::false_type) {}
+ return result;
+ }
+ }
//! If the user has specified a global capillary pressure, check if it is lower than the given entry capillary pressure.
//! This may be needed to exactly reproduce pc-S curves given by static network models.
template<class Scalar>
- bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pcEntry, std::true_type) const
+ bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pcEntry) const
{
- return restrictToGlobalCapillaryPressure_ && (pcEntry > problem_.globalCapillaryPressure());
+ if constexpr (hasGlobalCapillaryPressure<Problem>())
+ return restrictToGlobalCapillaryPressure_ && (pcEntry > problem_.globalCapillaryPressure());
+ else
+ return false;
}
- //! Do nothing here.
- template<class Scalar>
- bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pc, std::false_type) const
- { return false; }
-
std::vector<bool> invadedCurrentIteration_;
std::vector<bool> invadedPreviousTimeStep_;
bool hasChangedInCurrentIteration_ = false;
diff --git a/dumux/porenetwork/2p/model.hh b/dumux/porenetwork/2p/model.hh
index 8d676cf2f0e81b14cff8683c430b6961fd2bfcda..16591c003b4b0813a75e0586e7f5a3993b7dffeb 100644
--- a/dumux/porenetwork/2p/model.hh
+++ b/dumux/porenetwork/2p/model.hh
@@ -130,10 +130,8 @@ struct GridFluxVariablesCache<TypeTag, TTag::PNMTwoP>
private:
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
using Problem = GetPropType<TypeTag, Properties::Problem>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
- using Labels = GetPropType<TypeTag, Properties::Labels>;
using FluxVariablesCache = GetPropType<TypeTag, Properties::FluxVariablesCache>;
- using Traits = Dumux::PoreNetwork::PNMTwoPDefaultGridFVCTraits<Problem, FluxVariablesCache, Indices, Labels>;
+ using Traits = Dumux::PoreNetwork::PNMTwoPDefaultGridFVCTraits<Problem, FluxVariablesCache>;
public:
using type = Dumux::PoreNetwork::PNMTwoPGridFluxVariablesCache<Problem, FluxVariablesCache, enableCache, Traits>;
};
diff --git a/dumux/porenetwork/2p/volumevariables.hh b/dumux/porenetwork/2p/volumevariables.hh
index dcb859bc171e433df525b08a19c4986b4a98e5f8..f05245028c7d1731a79f41ac72b9380d0c3b77cf 100644
--- a/dumux/porenetwork/2p/volumevariables.hh
+++ b/dumux/porenetwork/2p/volumevariables.hh
@@ -76,10 +76,7 @@ public:
ParentType::update(elemSol, problem, element, scv);
poreInscribedRadius_ = problem.spatialParams().poreInscribedRadius(element, scv, elemSol);
poreVolume_ = problem.gridGeometry().poreVolume(scv.dofIndex()) * this->porosity();
-
- // the value of water/air TODO make general in fluid system
- static const Scalar gamma = getParamFromGroup<Scalar>(problem.paramGroup(), "Problem.SurfaceTension", 0.0725);
- surfaceTension_ = gamma;
+ surfaceTension_ = problem.spatialParams().surfaceTension(element, scv, elemSol);
}
/*!
diff --git a/dumux/porenetwork/properties.hh b/dumux/porenetwork/properties.hh
index 628182f78bb93359ae9ecf54b2919fde000102a3..dec06d8a0a3f017e099cf77919a7cbaa78d946cb 100644
--- a/dumux/porenetwork/properties.hh
+++ b/dumux/porenetwork/properties.hh
@@ -68,7 +68,7 @@ public:
};
template<class TypeTag>
-struct HeatConductionType<TypeTag, TTag::PoreNetworkModel> { using type = Dumux::PoreNetwork::PNMFouriersLaw; };
+struct HeatConductionType<TypeTag, TTag::PoreNetworkModel> { using type = Dumux::PoreNetwork::PNMFouriersLaw<>; };
//! The labels
template<class TypeTag>