Commit 0fef5ff3 authored by Hanchuan Wu's avatar Hanchuan Wu
Browse files

[pnm] Fix invasion state logic, update at the end of time step

parent 31346f62
......@@ -527,6 +527,7 @@ public:
// is capable of handling MultiType matrices or not
bool converged = solveLinearSystem_(deltaU);
// make sure all processes converged
int convergedRemote = converged;
if (comm_.size() > 1)
......@@ -670,13 +671,15 @@ public:
// When the Newton iterations are done: ask the model to check whether it makes sense
// TODO: how do we realize this? -> do this here in the Newton solver
// model_().checkPlausibility();
}
/*!
* \brief Called if the Newton method ended
* (not known yet if we failed or succeeded)
* (for pore-network model the invasion status will be updated here)
*/
virtual void newtonEnd() {}
virtual void newtonEnd(Variables &vars, const SolutionVector &uLastIter) {}
/*!
* \brief Returns true if the error of the solution is below the
......@@ -1028,7 +1031,7 @@ private:
}
// tell solver we are done
newtonEnd();
newtonEnd(vars, uLastIter);
// reset state if Newton failed
if (!newtonConverged())
......
......@@ -55,17 +55,17 @@ public:
TwoPInvasionState(const Problem& problem) : problem_(problem)
{
// initialize the invasion state
invadedCurrentIteration_.resize(problem.gridGeometry().gridView().size(0));
invadedCurrentTimeStep_.resize(problem.gridGeometry().gridView().size(0));
invadedPreviousTimeStep_.resize(problem.gridGeometry().gridView().size(0));
for (auto&& element : elements(problem.gridGeometry().gridView()))
{
const auto eIdx = problem.gridGeometry().elementMapper().index(element);
invadedCurrentIteration_[eIdx] = problem.initialInvasionState(element);
invadedPreviousTimeStep_[eIdx] = invadedCurrentIteration_[eIdx];
invadedCurrentTimeStep_[eIdx] = problem.initialInvasionState(element);
invadedPreviousTimeStep_[eIdx] = invadedCurrentTimeStep_[eIdx];
}
numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true);
numThroatsInvaded_ = std::count(invadedCurrentTimeStep_.begin(), invadedCurrentTimeStep_.end(), true);
verbose_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.Verbosity", true);
restrictToGlobalCapillaryPressure_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.RestrictInvasionToGlobalCapillaryPressure", false);
......@@ -84,7 +84,7 @@ public:
bool invaded(const Element& element) const
{
const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
return invadedCurrentIteration_[eIdx];
return invadedCurrentTimeStep_[eIdx];
}
//! Return the number of currently invaded throats
......@@ -93,9 +93,8 @@ public:
//! Update the invasion state of all throats. This is done after each Newton step by a call from the Newton solver.
template<class SolutionVector, class GridVolumeVariables, class GridFluxVariablesCache>
bool update(const SolutionVector& sol, const GridVolumeVariables& gridVolVars, GridFluxVariablesCache& gridFluxVarsCache)
void update(const SolutionVector& sol, const GridVolumeVariables& gridVolVars, GridFluxVariablesCache& gridFluxVarsCache)
{
hasChangedInCurrentIteration_ = false;
auto fvGeometry = localView(problem_.gridGeometry());
auto elemVolVars = localView(gridVolVars);
auto elemFluxVarsCache = localView(gridFluxVarsCache);
......@@ -110,39 +109,34 @@ public:
// checks if invasion or snap-off occured after Newton iteration step
if (const auto invasionResult = invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf]); invasionResult)
{
hasChangedInCurrentIteration_ = true;
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]);
gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invadedCurrentTimeStep_[eIdx]);
}
}
}
}
numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true);
return hasChangedInCurrentIteration_;
numThroatsInvaded_ = std::count(invadedCurrentTimeStep_.begin(), invadedCurrentTimeStep_.end(), true);
}
//! Restore the old invasion state after a Newton iteration has failed.
void reset()
{
hasChangedInCurrentIteration_ = false;
invadedCurrentIteration_ = invadedPreviousTimeStep_;
invadedCurrentTimeStep_ = invadedPreviousTimeStep_;
}
//! Return whether an invasion or snap-off occurred anywhere. Can be used, e.g., for output file writing control.
bool hasChanged() const
{ return hasChangedComparedToPreviousTimestep_; }
//! Return whether an invasion or snap-off occurred anywhere during the current Newton iteration.
bool hasChangedInCurrentIteration() const
{ return hasChangedInCurrentIteration_; }
//! This is called after the Newton method has successfully finished one time step.
void advance()
{
hasChangedComparedToPreviousTimestep_ = (invadedPreviousTimeStep_ != invadedCurrentIteration_);
invadedPreviousTimeStep_ = invadedCurrentIteration_;
hasChangedComparedToPreviousTimestep_ = (invadedPreviousTimeStep_ != invadedCurrentTimeStep_);
invadedPreviousTimeStep_ = invadedCurrentTimeStep_;
}
template<class SolutionVector, class GridVolumeVariables, class GridFluxVariablesCache>
......@@ -163,7 +157,7 @@ public:
{
// Only consider throats which have been invaded during the current time step
const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
if (!invadedCurrentIteration_[eIdx] || invadedPreviousTimeStep_[eIdx] == invadedCurrentIteration_[eIdx])
if (!invadedCurrentTimeStep_[eIdx] || invadedPreviousTimeStep_[eIdx] == invadedCurrentTimeStep_[eIdx])
continue;
fvGeometry.bindElement(element);
......@@ -186,7 +180,7 @@ public:
private:
//! The switch for determining the invasion state of a pore throat. Called at the end of each Newton step.
//! The switch for determining the invasion state of a pore throat. Called at the end of each time step
template<class Element, class ElementVolumeVariables, class FluxVariablesCache>
auto invasionSwitch_(const Element& element,
const ElementVolumeVariables& elemVolVars,
......@@ -197,7 +191,7 @@ private:
const auto& gridGeometry = problem_.gridGeometry();
const auto& spatialParams = problem_.spatialParams();
const auto eIdx = gridGeometry.elementMapper().index(element);
bool invadedBeforeSwitch = invadedCurrentIteration_[eIdx];
bool invadedBeforeSwitch = invadedCurrentTimeStep_[eIdx];
bool invadedAfterSwitch = invadedBeforeSwitch;
// Result type, containing the local scv index of the pore from which the invasion/snap-off occurred
......@@ -216,7 +210,7 @@ private:
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;
invadedCurrentTimeStep_[eIdx] = false;
return Result{}; // nothing happened
}
......@@ -235,7 +229,7 @@ private:
std::cout << ". pcEntry: " << spatialParams.pcEntry(element, elemVolVars) << std::endl;
}
invadedCurrentIteration_[eIdx] = false;
invadedCurrentTimeStep_[eIdx] = false;
return Result{}; //nothing happened
}
......@@ -244,7 +238,7 @@ private:
else if (*pcMax <= pcSnapoff)
invadedAfterSwitch = false;
invadedCurrentIteration_[eIdx] = invadedAfterSwitch;
invadedCurrentTimeStep_[eIdx] = invadedAfterSwitch;
if (invadedBeforeSwitch == invadedAfterSwitch)
return Result{}; // nothing happened
......@@ -291,9 +285,8 @@ private:
return false;
}
std::vector<bool> invadedCurrentIteration_;
std::vector<bool> invadedCurrentTimeStep_;
std::vector<bool> invadedPreviousTimeStep_;
bool hasChangedInCurrentIteration_ = false;
bool hasChangedComparedToPreviousTimestep_ = false;
std::size_t numThroatsInvaded_;
bool verbose_;
......
......@@ -36,37 +36,40 @@ namespace Dumux::PoreNetwork {
* which allows the newton method to abort quicker if the solution is
* way out of bounds.
*/
template<class Assembler, class LinearSolver,
template<class, class> class NewtonConsistencyChecks = TwoPNewtonConsistencyChecks>
class TwoPNewtonSolver : public Dumux::NewtonSolver<Assembler, LinearSolver>
template<class Assembler, class LinearSolver, bool useRegularization = true>
class TwoPNewtonSolver;
template<class Assembler, class LinearSolver>
class TwoPNewtonSolver<Assembler, LinearSolver, true>: public Dumux::NewtonSolver<Assembler, LinearSolver>
{
using ParentType = Dumux::NewtonSolver<Assembler, LinearSolver>;
using SolutionVector = typename Assembler::ResidualType;
public:
using ParentType::ParentType;
/*!
* \brief Called after each Newton update
*
* \param uCurrentIter The current global solution vector
* \param uLastIter The previous global solution vector
*/
void newtonEndStep(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter) final
void newtonEnd(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter) final
{
// call the method of the base class
ParentType::newtonEndStep(uCurrentIter, uLastIter);
ParentType::newtonEnd(uCurrentIter, uLastIter);
auto& gridVariables = this->assembler().gridVariables();
auto& invasionState = gridVariables.gridFluxVarsCache().invasionState();
switchedInLastIteration_ = invasionState.update(uCurrentIter, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
// here the invasion state is updated
invasionState.update(uCurrentIter, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
// If the solution is about to be accepted, check for accuracy and trigger a retry
// with a decreased time step size if necessary.
if (newtonConverged())
{
NewtonConsistencyChecks<typename Assembler::GridVariables, SolutionVector> checks;
TwoPNewtonConsistencyChecks<typename Assembler::GridVariables, SolutionVector> checks;
checks.performChecks(gridVariables, uCurrentIter, this->assembler().prevSol());
}
}
......@@ -77,12 +80,7 @@ public:
* switch was triggered in the last iteration.
*/
bool newtonConverged() const final
{
if (switchedInLastIteration_)
return false;
return ParentType::newtonConverged();
}
{ return ParentType::newtonConverged(); }
/*!
* \brief Called if the Newton method broke down.
......@@ -104,11 +102,7 @@ public:
auto& gridVariables = this->assembler().gridVariables();
gridVariables.gridFluxVarsCache().invasionState().advance();
}
private:
bool switchedInLastIteration_{false};
};
} // end namespace Dumux::PoreNetwork
#endif
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