diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh
index 6e58d391e27214bd8ace91eeb3419cb73d4c73ab..42d231f2faf8eb4cec367b5553ea1c1a733aa3d7 100644
--- a/dumux/porousmediumflow/richards/localresidual.hh
+++ b/dumux/porousmediumflow/richards/localresidual.hh
@@ -485,7 +485,11 @@ public:
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
- { /* TODO maybe forward to problem for the user to implement the Robin derivatives?*/ }
+ {
+ /* forward to problem for the user to implement the Robin derivatives*/
+ problem.addRobinFluxDerivatives(derivativeMatrices, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
+
+ }
private:
Implementation *asImp_()
diff --git a/test/porousmediumflow/richards/benchmarks/main.cc b/test/porousmediumflow/richards/benchmarks/main.cc
index bca77d5a3029f1c35a987e2932f5c5685d972703..aeb226cd0dbcd1e62a28850d5414807d7a5363c9 100644
--- a/test/porousmediumflow/richards/benchmarks/main.cc
+++ b/test/porousmediumflow/richards/benchmarks/main.cc
@@ -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");
diff --git a/test/porousmediumflow/richards/benchmarks/problem.hh b/test/porousmediumflow/richards/benchmarks/problem.hh
index 477f84114b11fe99834ffabe439e856d201b38b2..d3a9e89de33bc6c3525e8fb603cbb1738fcc9f2d 100644
--- a/test/porousmediumflow/richards/benchmarks/problem.hh
+++ b/test/porousmediumflow/richards/benchmarks/problem.hh
@@ -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_;
diff --git a/test/porousmediumflow/richards/lens/problem.hh b/test/porousmediumflow/richards/lens/problem.hh
index b7f8e093597909b910043bfd39b438b338e592d8..40599e10ffa1a505a088d9da00ed174bfea1011a 100644
--- a/test/porousmediumflow/richards/lens/problem.hh
+++ b/test/porousmediumflow/richards/lens/problem.hh
@@ -73,6 +73,8 @@ class RichardsLensProblem : public PorousMediumFlowProblem<TypeTag>
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
enum {
// copy some indices for convenience
pressureIdx = Indices::pressureIdx,
@@ -192,6 +194,25 @@ public:
// \}
+ /*!
+ * \brief Adds Robin flux derivatives for wetting phase. This is needed in case of solution dependent Neumann conditions.
+ *
+ * \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
+ {}
+
private:
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{