Commit 86e06f89 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Kilian Weishaupt
Browse files

[porousmediumflow][velocity] Improve calculation in case of Neumann BCs

* check if we have non-compositional model
* estimate the volume flow by neumannFlow/density
parent 346e62cf
......@@ -26,6 +26,7 @@
#define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
#include <vector>
#include <type_traits>
#include <dune/common/fvector.hh>
#include <dune/common/float_cmp.hh>
......@@ -38,6 +39,27 @@
namespace Dumux {
#ifndef DOXYGEN
namespace Detail {
// helper structs and functions detecting if the model is compositional
template <typename T, typename ...Ts>
using MoleFractionDetector = decltype(std::declval<T>().moleFraction(std::declval<Ts>()...));
template<class T, typename ...Args>
static constexpr bool hasMoleFraction()
{ return Dune::Std::is_detected<MoleFractionDetector, T, Args...>::value; }
template <typename T, typename ...Ts>
using MassFractionDetector = decltype(std::declval<T>().massFraction(std::declval<Ts>()...));
template<class T, typename ...Args>
static constexpr bool hasMassFraction()
{ return Dune::Std::is_detected<MassFractionDetector, T, Args...>::value; }
} // end namespace Detail
#endif // DOXYGEN
/*!
* \ingroup PorousmediumflowModels
* \brief Velocity computation for implicit (porous media) models.
......@@ -70,6 +92,11 @@ class PorousMediumFlowVelocity
using Problem = typename GridVolumeVariables::Problem;
static constexpr bool modelIsCompositional = Detail::hasMoleFraction<typename GridVolumeVariables::VolumeVariables, int, int>() ||
Detail::hasMassFraction<typename GridVolumeVariables::VolumeVariables, int, int>();
static_assert(VolumeVariables::numFluidPhases() >= 1, "Velocity output only makes sense for models with fluid phases.");
public:
static constexpr int numFluidPhases = VolumeVariables::numFluidPhases();
using VelocityVector = std::vector<Dune::FieldVector<Scalar, dimWorld>>;
......@@ -303,13 +330,26 @@ public:
if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
// otherwise, we try some reconstruction (TODO: Can this be improved?)
// otherwise, we try some reconstruction
// for cubes
else if (dim == 1 || geomType.isCube())
{
const auto fIdx = scvfIndexInInside[localScvfIdx];
const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
if constexpr (!modelIsCompositional)
{
// We assume that the density at the face equals the one at the cell center and reconstruct a volume flux from the Neumann mass flux.
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto eqIdx = VolumeVariables::Indices::conti0EqIdx + phaseIdx;
scvfFluxes[fIdx] += neumannFlux[eqIdx] / insideVolVars.density(phaseIdx) * scvf.area() * insideVolVars.extrusionFactor();
}
else
{
// For compositional models, we generally can't reconstruct the volume flow from the Neumann flux (which is a species flux rather
// than a phase flux here). Instead, we use the velocity of the opposing face.
const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
}
}
// for simplices
......
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