Commit b94fa769 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'fix/richards-analytic-storage-derivative' into 'master'

Fix/richards analytic storage derivative

Closes #981

See merge request !2574
parents 8d375bab b56ad1ff
Pipeline #3789 failed with stages
in 0 seconds
......@@ -34,6 +34,18 @@
#include <dumux/discretization/extrusion.hh>
#include <dumux/flux/referencesystemformulation.hh>
namespace Dumux::Detail {
// helper structs and functions detecting if the user-defined problem class
// implements addRobinFluxDerivatives
template <typename T, typename ...Ts>
using RobinDerivDetector = decltype(std::declval<T>().addRobinFluxDerivatives(std::declval<Ts>()...));
template<class T, typename ...Args>
static constexpr bool hasAddRobinFluxDerivatives()
{ return Dune::Std::is_detected<RobinDerivDetector, T, Args...>::value; }
} // end namespace Dumux::Detail
namespace Dumux {
/*!
......@@ -196,7 +208,7 @@ public:
static_assert(!FluidSystem::isCompressible(0),
"richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
const auto poreVolume = Extrusion::volume(scv)*curVolVars.porosity();
const auto poreVolume = Extrusion::volume(scv)*curVolVars.porosity()*curVolVars.extrusionFactor();
static const auto rho = curVolVars.density(0);
// partial derivative of storage term w.r.t. p_w
......@@ -485,7 +497,13 @@ public:
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{ /* TODO maybe forward to problem for the user to implement the Robin derivatives?*/ }
{
if constexpr(Detail::hasAddRobinFluxDerivatives<Problem,
PartialDerivativeMatrices, Element, FVElementGeometry,
ElementVolumeVariables, ElementFluxVariablesCache, SubControlVolumeFace>()
)
problem.addRobinFluxDerivatives(derivativeMatrices, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
}
private:
Implementation *asImp_()
......
......@@ -98,8 +98,7 @@ int main(int argc, char** argv)
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
// using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>; // TODO: something is wrong with analytic derivatives!
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
const auto dt = getParam<double>("TimeLoop.DtInitial");
const auto checkPoints = getParam<std::vector<double>>("TimeLoop.TEnd");
......
......@@ -175,6 +175,7 @@ public:
// kg/m^3 * m^2 * Pa / m / Pa / s = kg/s/m^2
auto criticalRate = density*K/viscosity*((cellPressure - criticalSurfacePressure_)/dist - density*gravity);
if (!std::signbit(criticalRate))
criticalRate *= useKrwAverage_ ? avgRelPerm : relPerm;
......@@ -229,6 +230,81 @@ public:
return rate*86400*1000/1000;
}
/*!
* \brief Adds Robin flux derivatives for wetting phase
*
* \param derivativeMatrices The matrices containing the derivatives
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curElemVolVars The current element volume variables
* \param elemFluxVarsCache The element flux variables cache
* \param scvf The sub control volume face
*/
template<class PartialDerivativeMatrices, class ElementVolumeVariables, class ElementFluxVariablesCache>
void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
const auto& globalPos = scvf.ipGlobal();
const auto insideFluidMatrixInteraction = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
//material law derivatives
const auto insideSw = insideVolVars.saturation(0);
const auto insidePc = insideVolVars.capillaryPressure();
const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
if (onUpperBoundary(globalPos))
{
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
const auto dist = (fvGeometry.scv(scvf.insideScvIdx()).center() - globalPos).two_norm();
const auto cellPressure = volVars.pressure(0);
const auto density = volVars.density(0);
const auto viscosity = volVars.viscosity(0);
const auto relPerm = useKrwAverage_ ? volVars.relativePermeability(0)*0.5 : volVars.relativePermeability(0);
const auto K = volVars.permeability();
const auto gravity = enableGravity_ ? 9.81 : 0.0;
// kg/m^3 * m^2 * Pa / m / Pa / s = kg/s/m^2
auto criticalRate = density*K/viscosity*((cellPressure - criticalSurfacePressure_)/dist - density*gravity);
if (!std::signbit(criticalRate))
criticalRate *= relPerm;
if (scenario_ == BenchmarkScenario::evaporation)
{
if (criticalRate <= potentialRate_)
derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0]
+= (density/viscosity*K*relPerm/dist + density*K/viscosity*((cellPressure - criticalSurfacePressure_)/dist - density*gravity) *dkrw_dsw_inside*dsw_dpw_inside)*surfaceArea_;
}
else //in case of infiltration no relative permeability is added in this term
{
if (criticalRate >= potentialRate_)
derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0] += density*K/viscosity/dist*surfaceArea_;
}
}
//free drainage (gravitational flux) for infiltration scenario
else if (onLowerBoundary(globalPos) && (scenario_ == BenchmarkScenario::infiltration))
{
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
const auto gravity = enableGravity_ ? 9.81 : 0.0;
const auto density = volVars.density(0);
const auto relPerm = volVars.relativePermeability(0);
const auto viscosity = volVars.viscosity(0);
const auto K = volVars.permeability();
derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0] += density*K*relPerm/viscosity*(density*gravity)*dkrw_dsw_inside*dsw_dpw_inside*surfaceArea_;
}
}
private:
Scalar initialPressure_, criticalSurfacePressure_, potentialRate_;
Scalar criticalSurfaceKrw_;
......
Markdown is supported
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