[tests][richards] change benchmark test to analytic derivatives

[richards][test] add robin flux derivative for analytic differentiation in lens and benchmark problem

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_()

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");

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_;

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
     {