GCC Code Coverage Report

Directory:	../../../builds/dumux-repositories/
File:	/builds/dumux-repositories/dumux/test/porousmediumflow/richards/benchmarks/problem.hh
Date:	2024-09-21 20:52:54

	Exec	Total	Coverage
Lines:	105	114	92.1%
Functions:	6	9	66.7%
Branches:	154	248	62.1%

  
      Line
      Branch
      Exec
      Source
    
      // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
    
      // vi: set et ts=4 sw=4 sts=4:
    
      //
    
      // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
    
      // SPDX-License-Identifier: GPL-3.0-or-later
    
      //
    
      /*!
    
       * \file
    
       * \ingroup RichardsTests
    
       * \brief Richards benchmarks base problem
    
       *
    
       * Infiltration benchmark:
    
       * Root-soil benchmark paper Schnepf et al. (case M2.1, Eq. 4) https://doi.org/10.3389/fpls.2020.00316
    
       * based on Vanderborght 2005 (see Fig. 4abc and Eq. 56-60) https://doi.org/10.2113/4.1.206
    
       *
    
       * Evaporation benchmark:
    
       * Root-soil benchmark paper Schnepf et al. (case M2.2) https://doi.org/10.3389/fpls.2020.00316
    
       * based on Vanderborght 2005 (see Fig. 5abcd and Eq. 39-47) https://doi.org/10.2113/4.1.206
    
       */
    
      #ifndef DUMUX_RICHARDS_BECHMARK_PROBLEM_HH
    
      #define DUMUX_RICHARDS_BECHMARK_PROBLEM_HH
    
      #include <dumux/common/properties.hh>
    
      #include <dumux/common/parameters.hh>
    
      #include <dumux/common/boundarytypes.hh>
    
      #include <dumux/common/numeqvector.hh>
    
      #include <dumux/io/format.hh>
    
      #include <dumux/porousmediumflow/problem.hh>
    
      #include <dumux/material/components/simpleh2o.hh>
    
      namespace Dumux {
    
      enum class BenchmarkScenario {
    
          evaporation, infiltration
    
      };
    
      /*!
    
       * \ingroup RichardsTests
    
       * \brief Richards benchmarks base problem
    
       */
    
      template <class TypeTag>
    
      class RichardsBenchmarkProblem : public PorousMediumFlowProblem<TypeTag>
    
      {
    
          using ParentType = PorousMediumFlowProblem<TypeTag>;
    
          using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    
          using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    
          using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
    
          using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
    
          using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    
          using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    
          using FVElementGeometry = typename GridGeometry::LocalView;
    
          using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
    
          using GridView = typename GridGeometry::GridView;
    
          using Element = typename GridView::template Codim<0>::Entity;
    
          using GlobalPosition = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
    
          static constexpr int dimWorld = GridView::dimensionworld;
    
      public:
    
      7
          RichardsBenchmarkProblem(std::shared_ptr<const GridGeometry> gridGeometry)
    
        7/18✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 7 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 7 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 7 times.
✓ Branch 15 taken 7 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 7 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.

      21
          : ParentType(gridGeometry)
    
          {
    
        2/4✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 7 times.

      7
              name_ = getParam<std::string>("Problem.Name");
    
      7
              const auto density = Components::SimpleH2O<double>::liquidDensity(0,0);
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
              const auto initialHead = getParam<Scalar>("Problem.InitialHeadInCm")*0.01;
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
              const auto criticalHead = getParam<Scalar>("Problem.CriticalSurfaceHeadInCm")*0.01;
    
      7
              initialPressure_ = 1.0e5 + initialHead*9.81*density;
    
      7
              criticalSurfacePressure_ = 1.0e5 + criticalHead*9.81*density;
    
        3/6✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.

      21
              const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(0.0);
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
              const auto criticalSaturation = fm.sw(1.0e5 - criticalSurfacePressure_);
    
      7
              criticalSurfaceKrw_ = fm.krw(criticalSaturation);
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
              enableGravity_ = getParam<bool>("Problem.EnableGravity", true);
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
              const auto potentialRate = getParam<Scalar>("Problem.SurfaceFluxMilliMeterPerDay");
    
      7
              potentialRate_ = density*potentialRate/(1000*86400.0); // mm/day -> kg/s/m^2
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
              useKrwAverage_ = getParam<bool>("Problem.UseKrwAverage", false);
    
        1/4✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.

      7
              bottomDirichlet_ = getParam<bool>("Problem.BottomDirichlet", false);
    
        2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.

      7
              scenario_ = (potentialRate > 0) ? BenchmarkScenario::evaporation : BenchmarkScenario::infiltration;
    
        2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.

      7
              surfaceArea_ = (scenario_ == BenchmarkScenario::evaporation) ? 0.1*0.1 : 0.05*0.05;
    
      7
          }
    
          // output name
    
          const std::string& name() const
    
        1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.

      7
          { return name_; }
    
          // reference pressure
    
      ✗
          Scalar nonwettingReferencePressure() const
    
      ✗
          { return 1.0e5; };
    
          /*!
    
           * \brief Specifies which kind of boundary condition should be
    
           *        used for which equation on a given boundary segment.
    
           */
    
      114282
          BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    
          {
    
        2/2✓ Branch 0 taken 57141 times.
✓ Branch 1 taken 57141 times.

      114282
              BoundaryTypes bcTypes;
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57141 times.

      57141
              if (onLowerBoundary(globalPos) && !bottomDirichlet_)
    
                  bcTypes.setAllNeumann();
    
      ✗
              else if (onLowerBoundary(globalPos) && bottomDirichlet_)
    
                  bcTypes.setAllDirichlet();
    
      114282
              else if (onUpperBoundary(globalPos))
    
                  bcTypes.setAllNeumann();
    
              else
    
      ✗
                  DUNE_THROW(Dune::InvalidStateException, "Wrong boundary?");
    
      114282
              return bcTypes;
    
          }
    
          /*!
    
           * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
    
           */
    
      ✗
          PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
    
      ✗
          { return initialAtPos(globalPos); }
    
          /*!
    
           * \brief Evaluates the initial values for a control volume
    
           */
    
      ✗
          PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
    
          {
    
      4900
              PrimaryVariables values(0.0);
    
      4900
              values[Indices::pressureIdx] = initialPressure_;
    
      ✗
              return values;
    
          }
    
          /*!
    
           * \brief Evaluates the boundary conditions for a Neumann boundary segment.
    
           * Negative values mean influx.
    
           */
    
          template<class ElementVolumeVariables, class ElementFluxVariablesCache>
    
      79912
          NumEqVector neumann(const Element& element,
    
                              const FVElementGeometry& fvGeometry,
    
                              const ElementVolumeVariables& elemVolVars,
    
                              const ElementFluxVariablesCache& elemFluxVarsCache,
    
                              const SubControlVolumeFace& scvf) const
    
          {
    
        2/2✓ Branch 0 taken 19978 times.
✓ Branch 1 taken 19978 times.

      79912
              NumEqVector values(0.0);
    
        2/2✓ Branch 0 taken 19978 times.
✓ Branch 1 taken 19978 times.

      79912
              const auto& globalPos = scvf.ipGlobal();
    
      159824
              if (onUpperBoundary(globalPos))
    
              {
    
      79912
                  const auto& volVars = elemVolVars[scvf.insideScvIdx()];
    
        8/10✓ Branch 0 taken 19978 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 19978 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4427 times.
✓ Branch 5 taken 15551 times.
✓ Branch 6 taken 4427 times.
✓ Branch 7 taken 15551 times.
✓ Branch 8 taken 4427 times.
✓ Branch 9 taken 15551 times.

      199780
                  const auto dist = (fvGeometry.scv(scvf.insideScvIdx()).center() - globalPos).two_norm();
    
        2/2✓ Branch 0 taken 4427 times.
✓ Branch 1 taken 15551 times.

      39956
                  const auto cellPressure = volVars.pressure(0);
    
        2/2✓ Branch 0 taken 4427 times.
✓ Branch 1 taken 15551 times.

      39956
                  const auto density = volVars.density(0);
    
        2/2✓ Branch 0 taken 4427 times.
✓ Branch 1 taken 15551 times.

      39956
                  const auto viscosity = volVars.viscosity(0);
    
        2/2✓ Branch 0 taken 4427 times.
✓ Branch 1 taken 15551 times.

      39956
                  const auto relPerm = volVars.relativePermeability(0);
    
      39956
                  const auto K = volVars.permeability();
    
        2/2✓ Branch 0 taken 4427 times.
✓ Branch 1 taken 15551 times.

      39956
                  const auto gravity = enableGravity_ ? 9.81 : 0.0;
    
      39956
                  const auto avgRelPerm = 0.5*(relPerm + criticalSurfaceKrw_);
    
                  // kg/m^3 * m^2 * Pa / m / Pa / s = kg/s/m^2
    
      39956
                  auto criticalRate = density*K/viscosity*((cellPressure - criticalSurfacePressure_)/dist - density*gravity);
    
        4/4✓ Branch 0 taken 4354 times.
✓ Branch 1 taken 15624 times.
✓ Branch 2 taken 4354 times.
✓ Branch 3 taken 15624 times.

      79912
                  if (!std::signbit(criticalRate))
    
        1/2✓ Branch 0 taken 4354 times.
✗ Branch 1 not taken.

      17416
                      criticalRate *= useKrwAverage_ ? avgRelPerm : relPerm;
    
        2/2✓ Branch 0 taken 4427 times.
✓ Branch 1 taken 15551 times.

      39956
                  if (scenario_ == BenchmarkScenario::evaporation)
    
        2/2✓ Branch 0 taken 3318 times.
✓ Branch 1 taken 1109 times.

      15490
                      values[Indices::conti0EqIdx] = std::min(potentialRate_, criticalRate);
    
                  else
    
        2/2✓ Branch 0 taken 9121 times.
✓ Branch 1 taken 6430 times.

      49344
                      values[Indices::conti0EqIdx] = std::max(potentialRate_, criticalRate);
    
              }
    
              // free drainage (gravitational flux) for infiltration scenario
    
        2/2✓ Branch 0 taken 15551 times.
✓ Branch 1 taken 4427 times.

      39956
              else if (onLowerBoundary(globalPos) && (scenario_ == BenchmarkScenario::infiltration))
    
              {
    
      62204
                  const auto& volVars = elemVolVars[scvf.insideScvIdx()];
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15551 times.

      31102
                  const auto gravity = enableGravity_ ? 9.81 : 0.0;
    
      31102
                  const auto density = volVars.density(0);
    
      31102
                  const auto viscosity = volVars.viscosity(0);
    
      31102
                  const auto relPerm = volVars.relativePermeability(0);
    
      31102
                  const auto K = volVars.permeability();
    
      31102
                  values[Indices::conti0EqIdx] = density*K*relPerm/viscosity*(density*gravity);
    
              }
    
      79912
              return values;
    
          }
    
          bool onLowerBoundary(const GlobalPosition &globalPos) const
    
        16/26✓ Branch 0 taken 19047 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 19047 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 19047 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 19047 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 57141 times.
✓ Branch 9 taken 57141 times.
✓ Branch 10 taken 57141 times.
✓ Branch 11 taken 57141 times.
✓ Branch 12 taken 57141 times.
✓ Branch 13 taken 57141 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 57141 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 57141 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 57141 times.
✓ Branch 20 taken 931 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 931 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 931 times.
✗ Branch 25 not taken.

      421827
          { return globalPos[dimWorld-1] < this->gridGeometry().bBoxMin()[dimWorld-1] + eps_; }
    
          bool onUpperBoundary(const GlobalPosition &globalPos) const
    
        17/20✓ Branch 0 taken 19047 times.
✓ Branch 1 taken 19047 times.
✓ Branch 2 taken 19047 times.
✓ Branch 3 taken 19047 times.
✓ Branch 4 taken 19047 times.
✓ Branch 5 taken 19047 times.
✓ Branch 6 taken 19047 times.
✓ Branch 7 taken 19047 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 57141 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 57141 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 57141 times.
✓ Branch 14 taken 931 times.
✓ Branch 15 taken 931 times.
✓ Branch 16 taken 931 times.
✓ Branch 17 taken 931 times.
✓ Branch 18 taken 931 times.
✓ Branch 19 taken 931 times.

      329385
          { return globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_; }
    
          //! compute the actual evaporation/infiltration rate
    
          template<class SolutionVector, class GridVariables>
    
      931
          Scalar computeActualRate(const SolutionVector& sol, const GridVariables& gridVars, bool verbose = true) const
    
          {
    
      931
              Scalar rate = 0.0;
    
        2/4✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 931 times.
✗ Branch 5 not taken.

      2793
              auto fvGeometry = localView(this->gridGeometry());
    
        2/4✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 931 times.
✗ Branch 5 not taken.

      1862
              auto elemVolVars = localView(gridVars.curGridVolVars());
    
        5/10✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 931 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 931 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 931931 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 931000 times.
✗ Branch 13 not taken.

      1864793
              for (const auto& element : elements(this->gridGeometry().gridView()))
    
              {
    
        1/2✓ Branch 1 taken 931000 times.
✗ Branch 2 not taken.

      931000
                  fvGeometry.bindElement(element);
    
        1/2✓ Branch 1 taken 931000 times.
✗ Branch 2 not taken.

      931000
                  elemVolVars.bindElement(element, fvGeometry, sol);
    
        4/4✓ Branch 0 taken 1862000 times.
✓ Branch 1 taken 931000 times.
✓ Branch 2 taken 1862000 times.
✓ Branch 3 taken 931000 times.

      3724000
                  for (const auto& scvf : scvfs(fvGeometry))
    
        2/2✓ Branch 0 taken 1862 times.
✓ Branch 1 taken 1860138 times.

      1862000
                      if (scvf.boundary())
    
      1862
                          rate += this->neumann(element, fvGeometry, elemVolVars, 0.0, scvf)[0];
    
              }
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 931 times.

      931
              if (verbose)
    
      ✗
                  std::cout << Fmt::format("Actual rate: {:.5g} (mm/day)\n", rate*86400*1000/1000);
    
      1862
              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>
    
      38094
          void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
    
                                       const Element& element,
    
                                       const FVElementGeometry& fvGeometry,
    
                                       const ElementVolumeVariables& elemVolVars,
    
                                       const ElementFluxVariablesCache& elemFluxVarsCache,
    
                                       const SubControlVolumeFace& scvf) const
    
          {
    
      38094
              const auto insideScvIdx = scvf.insideScvIdx();
    
      38094
              const auto& insideVolVars = elemVolVars[insideScvIdx];
    
        2/2✓ Branch 0 taken 31323 times.
✓ Branch 1 taken 6771 times.

      38094
              const auto& globalPos = scvf.ipGlobal();
    
        4/4✓ Branch 0 taken 31323 times.
✓ Branch 1 taken 6771 times.
✓ Branch 2 taken 31323 times.
✓ Branch 3 taken 6771 times.

      76188
              const auto insideFluidMatrixInteraction = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
    
              //material law derivatives
    
        2/2✓ Branch 0 taken 31323 times.
✓ Branch 1 taken 6771 times.

      38094
              const auto insideSw = insideVolVars.saturation(0);
    
        2/2✓ Branch 0 taken 31323 times.
✓ Branch 1 taken 6771 times.

      38094
              const auto insidePc = insideVolVars.capillaryPressure();
    
      38094
              const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
    
      38094
              const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
    
        4/4✓ Branch 0 taken 19047 times.
✓ Branch 1 taken 19047 times.
✓ Branch 2 taken 19047 times.
✓ Branch 3 taken 19047 times.

      76188
              if (onUpperBoundary(globalPos))
    
              {
    
      38094
                  const auto& volVars = elemVolVars[scvf.insideScvIdx()];
    
        5/10✓ Branch 0 taken 19047 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 19047 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 19047 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 19047 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 19047 times.

      95235
                  const auto dist = (fvGeometry.scv(scvf.insideScvIdx()).center() - globalPos).two_norm();
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 19047 times.

      19047
                  const auto cellPressure = volVars.pressure(0);
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 19047 times.

      19047
                  const auto density = volVars.density(0);
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 19047 times.

      19047
                  const auto viscosity = volVars.viscosity(0);
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 19047 times.

      19047
                  const auto relPerm = useKrwAverage_ ?  volVars.relativePermeability(0)*0.5 : volVars.relativePermeability(0);
    
      19047
                  const auto K = volVars.permeability();
    
        2/2✓ Branch 0 taken 3496 times.
✓ Branch 1 taken 15551 times.

      19047
                  const auto gravity = enableGravity_ ? 9.81 : 0.0;
    
                  // kg/m^3 * m^2 * Pa / m / Pa / s = kg/s/m^2
    
      19047
                  auto criticalRate = density*K/viscosity*((cellPressure - criticalSurfacePressure_)/dist - density*gravity);
    
        4/4✓ Branch 0 taken 3423 times.
✓ Branch 1 taken 15624 times.
✓ Branch 2 taken 3423 times.
✓ Branch 3 taken 15624 times.

      38094
                  if (!std::signbit(criticalRate))
    
      3423
                      criticalRate *= relPerm;
    
        2/2✓ Branch 0 taken 3496 times.
✓ Branch 1 taken 15551 times.

      19047
                  if (scenario_ == BenchmarkScenario::evaporation)
    
                  {
    
        2/2✓ Branch 0 taken 2522 times.
✓ Branch 1 taken 974 times.

      3496
                      if (criticalRate <= potentialRate_)
    
      2522
                          derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0]
    
      2522
                              += (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
    
                  {
    
        2/2✓ Branch 0 taken 9121 times.
✓ Branch 1 taken 6430 times.

      15551
                      if (criticalRate >= potentialRate_)
    
      9121
                          derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0] += density*K/viscosity/dist*surfaceArea_;
    
                  }
    
              }
    
              //free drainage (gravitational flux) for infiltration scenario
    
        4/6✓ Branch 0 taken 19047 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 19047 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 15551 times.
✓ Branch 5 taken 3496 times.

      38094
              else if (onLowerBoundary(globalPos) && (scenario_ == BenchmarkScenario::infiltration))
    
              {
    
      31102
                  const auto& volVars = elemVolVars[scvf.insideScvIdx()];
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15551 times.

      15551
                  const auto gravity = enableGravity_ ? 9.81 : 0.0;
    
      15551
                  const auto density = volVars.density(0);
    
      15551
                  const auto relPerm = volVars.relativePermeability(0);
    
      15551
                  const auto viscosity = volVars.viscosity(0);
    
      15551
                  const auto K = volVars.permeability();
    
      15551
                  derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0] += density*K*relPerm/viscosity*(density*gravity)*dkrw_dsw_inside*dsw_dpw_inside*surfaceArea_;
    
              }
    
      38094
          }
    
      private:
    
          Scalar initialPressure_, criticalSurfacePressure_, potentialRate_;
    
          Scalar criticalSurfaceKrw_;
    
          static constexpr Scalar eps_ = 1.5e-7;
    
          std::string name_;
    
          bool enableGravity_;
    
          bool bottomDirichlet_;
    
          bool useKrwAverage_;
    
          BenchmarkScenario scenario_;
    
          Scalar surfaceArea_;
    
      };
    
      } // end namespace Dumux
    
      #endif

Line	Branch	Exec	Source
1			// -- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 --
2			// vi: set et ts=4 sw=4 sts=4:
3			//
4			// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5			// SPDX-License-Identifier: GPL-3.0-or-later
6			//
7			/*!
8			* \file
9			* \ingroup RichardsTests
10			* \brief Richards benchmarks base problem
11			*
12			* Infiltration benchmark:
13			* Root-soil benchmark paper Schnepf et al. (case M2.1, Eq. 4) https://doi.org/10.3389/fpls.2020.00316
14			* based on Vanderborght 2005 (see Fig. 4abc and Eq. 56-60) https://doi.org/10.2113/4.1.206
15			*
16			* Evaporation benchmark:
17			* Root-soil benchmark paper Schnepf et al. (case M2.2) https://doi.org/10.3389/fpls.2020.00316
18			* based on Vanderborght 2005 (see Fig. 5abcd and Eq. 39-47) https://doi.org/10.2113/4.1.206
19			*/
20
21			#ifndef DUMUX_RICHARDS_BECHMARK_PROBLEM_HH
22			#define DUMUX_RICHARDS_BECHMARK_PROBLEM_HH
23
24			#include <dumux/common/properties.hh>
25			#include <dumux/common/parameters.hh>
26			#include <dumux/common/boundarytypes.hh>
27			#include <dumux/common/numeqvector.hh>
28			#include <dumux/io/format.hh>
29
30			#include <dumux/porousmediumflow/problem.hh>
31			#include <dumux/material/components/simpleh2o.hh>
32
33			namespace Dumux {
34
35			enum class BenchmarkScenario {
36			evaporation, infiltration
37			};
38
39			/*!
40			* \ingroup RichardsTests
41			* \brief Richards benchmarks base problem
42			*/
43			template <class TypeTag>
44			class RichardsBenchmarkProblem : public PorousMediumFlowProblem<TypeTag>
45			{
46			using ParentType = PorousMediumFlowProblem<TypeTag>;
47
48			using Scalar = GetPropType<TypeTag, Properties::Scalar>;
49			using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
50			using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
51			using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
52
53			using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
54
55			using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
56			using FVElementGeometry = typename GridGeometry::LocalView;
57			using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
58			using GridView = typename GridGeometry::GridView;
59			using Element = typename GridView::template Codim<0>::Entity;
60			using GlobalPosition = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
61
62			static constexpr int dimWorld = GridView::dimensionworld;
63
64			public:
65		7	RichardsBenchmarkProblem(std::shared_ptr<const GridGeometry> gridGeometry)
66	7/18 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 7 times. ✗ Branch 5 not taken. ✓ Branch 8 taken 7 times. ✗ Branch 9 not taken. ✓ Branch 10 taken 7 times. ✗ Branch 11 not taken. ✗ Branch 12 not taken. ✓ Branch 13 taken 7 times. ✓ Branch 15 taken 7 times. ✗ Branch 16 not taken. ✓ Branch 18 taken 7 times. ✗ Branch 19 not taken. ✗ Branch 20 not taken. ✗ Branch 21 not taken. ✗ Branch 22 not taken. ✗ Branch 23 not taken.	21	: ParentType(gridGeometry)
67			{
68	2/4 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken. ✗ Branch 4 not taken. ✓ Branch 5 taken 7 times.	7	name_ = getParam<std::string>("Problem.Name");
69		7	const auto density = Components::SimpleH2O<double>::liquidDensity(0,0);
70	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	const auto initialHead = getParam<Scalar>("Problem.InitialHeadInCm")*0.01;
71	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	const auto criticalHead = getParam<Scalar>("Problem.CriticalSurfaceHeadInCm")*0.01;
72		7	initialPressure_ = 1.0e5 + initialHead9.81density;
73		7	criticalSurfacePressure_ = 1.0e5 + criticalHead9.81density;
74	3/6 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 7 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 7 times. ✗ Branch 8 not taken.	21	const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(0.0);
75	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	const auto criticalSaturation = fm.sw(1.0e5 - criticalSurfacePressure_);
76		7	criticalSurfaceKrw_ = fm.krw(criticalSaturation);
77	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	enableGravity_ = getParam<bool>("Problem.EnableGravity", true);
78
79	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	const auto potentialRate = getParam<Scalar>("Problem.SurfaceFluxMilliMeterPerDay");
80		7	potentialRate_ = densitypotentialRate/(100086400.0); // mm/day -> kg/s/m^2
81	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	useKrwAverage_ = getParam<bool>("Problem.UseKrwAverage", false);
82	1/4 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken. ✗ Branch 3 not taken. ✗ Branch 4 not taken.	7	bottomDirichlet_ = getParam<bool>("Problem.BottomDirichlet", false);
83	2/2 ✓ Branch 0 taken 3 times. ✓ Branch 1 taken 4 times.	7	scenario_ = (potentialRate > 0) ? BenchmarkScenario::evaporation : BenchmarkScenario::infiltration;
84	2/2 ✓ Branch 0 taken 3 times. ✓ Branch 1 taken 4 times.	7	surfaceArea_ = (scenario_ == BenchmarkScenario::evaporation) ? 0.10.1 : 0.050.05;
85		7	}
86
87			// output name
88			const std::string& name() const
89	1/2 ✓ Branch 1 taken 7 times. ✗ Branch 2 not taken.	7	{ return name_; }
90
91			// reference pressure
92		✗	Scalar nonwettingReferencePressure() const
93		✗	{ return 1.0e5; };
94
95			/*!
96			* \brief Specifies which kind of boundary condition should be
97			* used for which equation on a given boundary segment.
98			*/
99		114282	BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
100			{
101	2/2 ✓ Branch 0 taken 57141 times. ✓ Branch 1 taken 57141 times.	114282	BoundaryTypes bcTypes;
102	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 57141 times.	57141	if (onLowerBoundary(globalPos) && !bottomDirichlet_)
103			bcTypes.setAllNeumann();
104		✗	else if (onLowerBoundary(globalPos) && bottomDirichlet_)
105			bcTypes.setAllDirichlet();
106		114282	else if (onUpperBoundary(globalPos))
107			bcTypes.setAllNeumann();
108			else
109		✗	DUNE_THROW(Dune::InvalidStateException, "Wrong boundary?");
110		114282	return bcTypes;
111			}
112
113			/*!
114			* \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
115			*/
116		✗	PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
117		✗	{ return initialAtPos(globalPos); }
118
119			/*!
120			* \brief Evaluates the initial values for a control volume
121			*/
122		✗	PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
123			{
124		4900	PrimaryVariables values(0.0);
125		4900	values[Indices::pressureIdx] = initialPressure_;
126		✗	return values;
127			}
128
129			/*!
130			* \brief Evaluates the boundary conditions for a Neumann boundary segment.
131			* Negative values mean influx.
132			*/
133			template<class ElementVolumeVariables, class ElementFluxVariablesCache>
134		79912	NumEqVector neumann(const Element& element,
135			const FVElementGeometry& fvGeometry,
136			const ElementVolumeVariables& elemVolVars,
137			const ElementFluxVariablesCache& elemFluxVarsCache,
138			const SubControlVolumeFace& scvf) const
139			{
140	2/2 ✓ Branch 0 taken 19978 times. ✓ Branch 1 taken 19978 times.	79912	NumEqVector values(0.0);
141	2/2 ✓ Branch 0 taken 19978 times. ✓ Branch 1 taken 19978 times.	79912	const auto& globalPos = scvf.ipGlobal();
142		159824	if (onUpperBoundary(globalPos))
143			{
144		79912	const auto& volVars = elemVolVars[scvf.insideScvIdx()];
145
146	8/10 ✓ Branch 0 taken 19978 times. ✗ Branch 1 not taken. ✓ Branch 2 taken 19978 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 4427 times. ✓ Branch 5 taken 15551 times. ✓ Branch 6 taken 4427 times. ✓ Branch 7 taken 15551 times. ✓ Branch 8 taken 4427 times. ✓ Branch 9 taken 15551 times.	199780	const auto dist = (fvGeometry.scv(scvf.insideScvIdx()).center() - globalPos).two_norm();
147	2/2 ✓ Branch 0 taken 4427 times. ✓ Branch 1 taken 15551 times.	39956	const auto cellPressure = volVars.pressure(0);
148	2/2 ✓ Branch 0 taken 4427 times. ✓ Branch 1 taken 15551 times.	39956	const auto density = volVars.density(0);
149	2/2 ✓ Branch 0 taken 4427 times. ✓ Branch 1 taken 15551 times.	39956	const auto viscosity = volVars.viscosity(0);
150	2/2 ✓ Branch 0 taken 4427 times. ✓ Branch 1 taken 15551 times.	39956	const auto relPerm = volVars.relativePermeability(0);
151		39956	const auto K = volVars.permeability();
152	2/2 ✓ Branch 0 taken 4427 times. ✓ Branch 1 taken 15551 times.	39956	const auto gravity = enableGravity_ ? 9.81 : 0.0;
153		39956	const auto avgRelPerm = 0.5*(relPerm + criticalSurfaceKrw_);
154
155			// kg/m^3 * m^2 * Pa / m / Pa / s = kg/s/m^2
156		39956	auto criticalRate = densityK/viscosity((cellPressure - criticalSurfacePressure_)/dist - density*gravity);
157
158	4/4 ✓ Branch 0 taken 4354 times. ✓ Branch 1 taken 15624 times. ✓ Branch 2 taken 4354 times. ✓ Branch 3 taken 15624 times.	79912	if (!std::signbit(criticalRate))
159	1/2 ✓ Branch 0 taken 4354 times. ✗ Branch 1 not taken.	17416	criticalRate *= useKrwAverage_ ? avgRelPerm : relPerm;
160
161	2/2 ✓ Branch 0 taken 4427 times. ✓ Branch 1 taken 15551 times.	39956	if (scenario_ == BenchmarkScenario::evaporation)
162	2/2 ✓ Branch 0 taken 3318 times. ✓ Branch 1 taken 1109 times.	15490	values[Indices::conti0EqIdx] = std::min(potentialRate_, criticalRate);
163			else
164	2/2 ✓ Branch 0 taken 9121 times. ✓ Branch 1 taken 6430 times.	49344	values[Indices::conti0EqIdx] = std::max(potentialRate_, criticalRate);
165			}
166
167			// free drainage (gravitational flux) for infiltration scenario
168	2/2 ✓ Branch 0 taken 15551 times. ✓ Branch 1 taken 4427 times.	39956	else if (onLowerBoundary(globalPos) && (scenario_ == BenchmarkScenario::infiltration))
169			{
170		62204	const auto& volVars = elemVolVars[scvf.insideScvIdx()];
171	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 15551 times.	31102	const auto gravity = enableGravity_ ? 9.81 : 0.0;
172		31102	const auto density = volVars.density(0);
173		31102	const auto viscosity = volVars.viscosity(0);
174		31102	const auto relPerm = volVars.relativePermeability(0);
175		31102	const auto K = volVars.permeability();
176
177		31102	values[Indices::conti0EqIdx] = densityKrelPerm/viscosity(densitygravity);
178			}
179
180		79912	return values;
181			}
182
183			bool onLowerBoundary(const GlobalPosition &globalPos) const
184	16/26 ✓ Branch 0 taken 19047 times. ✗ Branch 1 not taken. ✓ Branch 2 taken 19047 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 19047 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 19047 times. ✗ Branch 7 not taken. ✓ Branch 8 taken 57141 times. ✓ Branch 9 taken 57141 times. ✓ Branch 10 taken 57141 times. ✓ Branch 11 taken 57141 times. ✓ Branch 12 taken 57141 times. ✓ Branch 13 taken 57141 times. ✗ Branch 14 not taken. ✓ Branch 15 taken 57141 times. ✗ Branch 16 not taken. ✓ Branch 17 taken 57141 times. ✗ Branch 18 not taken. ✓ Branch 19 taken 57141 times. ✓ Branch 20 taken 931 times. ✗ Branch 21 not taken. ✓ Branch 22 taken 931 times. ✗ Branch 23 not taken. ✓ Branch 24 taken 931 times. ✗ Branch 25 not taken.	421827	{ return globalPos[dimWorld-1] < this->gridGeometry().bBoxMin()[dimWorld-1] + eps_; }
185
186			bool onUpperBoundary(const GlobalPosition &globalPos) const
187	17/20 ✓ Branch 0 taken 19047 times. ✓ Branch 1 taken 19047 times. ✓ Branch 2 taken 19047 times. ✓ Branch 3 taken 19047 times. ✓ Branch 4 taken 19047 times. ✓ Branch 5 taken 19047 times. ✓ Branch 6 taken 19047 times. ✓ Branch 7 taken 19047 times. ✗ Branch 8 not taken. ✓ Branch 9 taken 57141 times. ✗ Branch 10 not taken. ✓ Branch 11 taken 57141 times. ✗ Branch 12 not taken. ✓ Branch 13 taken 57141 times. ✓ Branch 14 taken 931 times. ✓ Branch 15 taken 931 times. ✓ Branch 16 taken 931 times. ✓ Branch 17 taken 931 times. ✓ Branch 18 taken 931 times. ✓ Branch 19 taken 931 times.	329385	{ return globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_; }
188
189			//! compute the actual evaporation/infiltration rate
190			template<class SolutionVector, class GridVariables>
191		931	Scalar computeActualRate(const SolutionVector& sol, const GridVariables& gridVars, bool verbose = true) const
192			{
193		931	Scalar rate = 0.0;
194
195	2/4 ✓ Branch 1 taken 931 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 931 times. ✗ Branch 5 not taken.	2793	auto fvGeometry = localView(this->gridGeometry());
196	2/4 ✓ Branch 1 taken 931 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 931 times. ✗ Branch 5 not taken.	1862	auto elemVolVars = localView(gridVars.curGridVolVars());
197
198	5/10 ✓ Branch 1 taken 931 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 931 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 931 times. ✗ Branch 8 not taken. ✓ Branch 9 taken 931931 times. ✗ Branch 10 not taken. ✓ Branch 12 taken 931000 times. ✗ Branch 13 not taken.	1864793	for (const auto& element : elements(this->gridGeometry().gridView()))
199			{
200	1/2 ✓ Branch 1 taken 931000 times. ✗ Branch 2 not taken.	931000	fvGeometry.bindElement(element);
201	1/2 ✓ Branch 1 taken 931000 times. ✗ Branch 2 not taken.	931000	elemVolVars.bindElement(element, fvGeometry, sol);
202	4/4 ✓ Branch 0 taken 1862000 times. ✓ Branch 1 taken 931000 times. ✓ Branch 2 taken 1862000 times. ✓ Branch 3 taken 931000 times.	3724000	for (const auto& scvf : scvfs(fvGeometry))
203	2/2 ✓ Branch 0 taken 1862 times. ✓ Branch 1 taken 1860138 times.	1862000	if (scvf.boundary())
204		1862	rate += this->neumann(element, fvGeometry, elemVolVars, 0.0, scvf)[0];
205			}
206
207	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 931 times.	931	if (verbose)
208		✗	std::cout << Fmt::format("Actual rate: {:.5g} (mm/day)\n", rate864001000/1000);
209
210		1862	return rate864001000/1000;
211			}
212
213			/*!
214			* \brief Adds Robin flux derivatives for wetting phase
215			*
216			* \param derivativeMatrices The matrices containing the derivatives
217			* \param element The element
218			* \param fvGeometry The finite volume element geometry
219			* \param curElemVolVars The current element volume variables
220			* \param elemFluxVarsCache The element flux variables cache
221			* \param scvf The sub control volume face
222			*/
223			template<class PartialDerivativeMatrices, class ElementVolumeVariables, class ElementFluxVariablesCache>
224		38094	void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
225			const Element& element,
226			const FVElementGeometry& fvGeometry,
227			const ElementVolumeVariables& elemVolVars,
228			const ElementFluxVariablesCache& elemFluxVarsCache,
229			const SubControlVolumeFace& scvf) const
230			{
231
232		38094	const auto insideScvIdx = scvf.insideScvIdx();
233		38094	const auto& insideVolVars = elemVolVars[insideScvIdx];
234	2/2 ✓ Branch 0 taken 31323 times. ✓ Branch 1 taken 6771 times.	38094	const auto& globalPos = scvf.ipGlobal();
235	4/4 ✓ Branch 0 taken 31323 times. ✓ Branch 1 taken 6771 times. ✓ Branch 2 taken 31323 times. ✓ Branch 3 taken 6771 times.	76188	const auto insideFluidMatrixInteraction = this->spatialParams().fluidMatrixInteractionAtPos(globalPos);
236
237			//material law derivatives
238	2/2 ✓ Branch 0 taken 31323 times. ✓ Branch 1 taken 6771 times.	38094	const auto insideSw = insideVolVars.saturation(0);
239	2/2 ✓ Branch 0 taken 31323 times. ✓ Branch 1 taken 6771 times.	38094	const auto insidePc = insideVolVars.capillaryPressure();
240		38094	const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
241		38094	const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
242
243	4/4 ✓ Branch 0 taken 19047 times. ✓ Branch 1 taken 19047 times. ✓ Branch 2 taken 19047 times. ✓ Branch 3 taken 19047 times.	76188	if (onUpperBoundary(globalPos))
244			{
245		38094	const auto& volVars = elemVolVars[scvf.insideScvIdx()];
246
247	5/10 ✓ Branch 0 taken 19047 times. ✗ Branch 1 not taken. ✓ Branch 2 taken 19047 times. ✗ Branch 3 not taken. ✗ Branch 4 not taken. ✓ Branch 5 taken 19047 times. ✗ Branch 6 not taken. ✓ Branch 7 taken 19047 times. ✗ Branch 8 not taken. ✓ Branch 9 taken 19047 times.	95235	const auto dist = (fvGeometry.scv(scvf.insideScvIdx()).center() - globalPos).two_norm();
248	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 19047 times.	19047	const auto cellPressure = volVars.pressure(0);
249	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 19047 times.	19047	const auto density = volVars.density(0);
250	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 19047 times.	19047	const auto viscosity = volVars.viscosity(0);
251	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 19047 times.	19047	const auto relPerm = useKrwAverage_ ? volVars.relativePermeability(0)*0.5 : volVars.relativePermeability(0);
252		19047	const auto K = volVars.permeability();
253	2/2 ✓ Branch 0 taken 3496 times. ✓ Branch 1 taken 15551 times.	19047	const auto gravity = enableGravity_ ? 9.81 : 0.0;
254
255			// kg/m^3 * m^2 * Pa / m / Pa / s = kg/s/m^2
256		19047	auto criticalRate = densityK/viscosity((cellPressure - criticalSurfacePressure_)/dist - density*gravity);
257
258	4/4 ✓ Branch 0 taken 3423 times. ✓ Branch 1 taken 15624 times. ✓ Branch 2 taken 3423 times. ✓ Branch 3 taken 15624 times.	38094	if (!std::signbit(criticalRate))
259		3423	criticalRate *= relPerm;
260
261	2/2 ✓ Branch 0 taken 3496 times. ✓ Branch 1 taken 15551 times.	19047	if (scenario_ == BenchmarkScenario::evaporation)
262			{
263	2/2 ✓ Branch 0 taken 2522 times. ✓ Branch 1 taken 974 times.	3496	if (criticalRate <= potentialRate_)
264		2522	derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0]
265		2522	+= (density/viscosityKrelPerm/dist + densityK/viscosity((cellPressure - criticalSurfacePressure_)/dist - densitygravity) dkrw_dsw_insidedsw_dpw_inside)surfaceArea_;
266			}
267			else //in case of infiltration no relative permeability is added in this term
268			{
269	2/2 ✓ Branch 0 taken 9121 times. ✓ Branch 1 taken 6430 times.	15551	if (criticalRate >= potentialRate_)
270		9121	derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0] += densityK/viscosity/distsurfaceArea_;
271			}
272			}
273
274			//free drainage (gravitational flux) for infiltration scenario
275	4/6 ✓ Branch 0 taken 19047 times. ✗ Branch 1 not taken. ✓ Branch 2 taken 19047 times. ✗ Branch 3 not taken. ✓ Branch 4 taken 15551 times. ✓ Branch 5 taken 3496 times.	38094	else if (onLowerBoundary(globalPos) && (scenario_ == BenchmarkScenario::infiltration))
276			{
277		31102	const auto& volVars = elemVolVars[scvf.insideScvIdx()];
278	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 15551 times.	15551	const auto gravity = enableGravity_ ? 9.81 : 0.0;
279		15551	const auto density = volVars.density(0);
280		15551	const auto relPerm = volVars.relativePermeability(0);
281		15551	const auto viscosity = volVars.viscosity(0);
282		15551	const auto K = volVars.permeability();
283
284		15551	derivativeMatrices[insideScvIdx][Indices::conti0EqIdx][0] += densityKrelPerm/viscosity(densitygravity)dkrw_dsw_insidedsw_dpw_inside*surfaceArea_;
285			}
286		38094	}
287
288			private:
289			Scalar initialPressure_, criticalSurfacePressure_, potentialRate_;
290			Scalar criticalSurfaceKrw_;
291			static constexpr Scalar eps_ = 1.5e-7;
292			std::string name_;
293			bool enableGravity_;
294			bool bottomDirichlet_;
295			bool useKrwAverage_;
296			BenchmarkScenario scenario_;
297			Scalar surfaceArea_;
298			};
299
300			} // end namespace Dumux
301
302			#endif
303