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>;
};