Commit 85a9d3aa authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[pnm2p] Clean up invasionstate

parent 5fb834c0
......@@ -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
......
......@@ -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;
......
......@@ -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>;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment