[examples] Set up problem for documentation (2p, adaptive, point source)

test/CMakeLists.txt

 add_subdirectory(common)
+add_subdirectory(examples)
 add_subdirectory(geomechanics)
 add_subdirectory(freeflow)
 add_subdirectory(io)

test/examples/2pinfiltration/CMakeLists.txt

+dune_symlink_to_source_files(FILES "params.input")
+dune_symlink_to_source_files(FILES "initialsolutioncc.txt")
+
+# using tpfa and point source
+dune_add_test(NAME test_2p_pointsource_adaptive
+              LABELS porousmediumflow 2p
+              SOURCES main.cc
+              CMAKE_GUARD dune-alugrid_FOUND
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS --script fuzzy
+                       --files ${CMAKE_SOURCE_DIR}/test/references/test_2p_pointsource_adaptive_tpfa-reference.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p_pointsource_adaptive_tpfa-00001.vtu
+                       --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_pointsource_adaptive_tpfa params.input -Problem.Name test_2p_pointsource_adaptive_tpfa")

test/examples/2pinfiltration/README.md

+This tutorial was copied from tests/porousmediumflow/2p/adaptive and restricted to the cell-centered finite volume TPFA discretization scheme.
+You need ALUGrid[0] in order to compile and run it.
+
+# Two-phase flow with infiltration and adaptive grid
+
+## Problem set-up
+Soil contamination problem where DNAPL infiltrates a fully water saturated medium.
+
+...
+
+
+## Infiltration (point source)
+
+## Adaptive grid
+
+
+
+[0]: https://gitlab.dune-project.org/extensions/dune-alugrid

test/examples/2pinfiltration/main.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPTests
+ * \brief Test for the two-phase porous medium flow model
+ */
+
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/method.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/adaptive/adapt.hh>
+#include <dumux/adaptive/markelements.hh>
+#include <dumux/adaptive/initializationindicator.hh>
+#include <dumux/porousmediumflow/2p/griddatatransfer.hh>
+#include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
+
+#include <test/porousmediumflow/2p/implicit/incompressible/problem.hh>
+#include "problem.hh"
+
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using TypeTag = Properties::TTag::TwoPAdaptivePointSource;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // try to create a grid (from the given grid file or the input file)
+    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+    gridManager.init();
+
+    ////////////////////////////////////////////////////////////
+    // run instationary non-linear problem on this grid
+    ////////////////////////////////////////////////////////////
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+    fvGridGeometry->update();
+
+    // the problem (initial and boundary conditions)
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(fvGridGeometry);
+    problem->computePointSourceMap(); // enable point sources
+
+    // the solution vector
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    SolutionVector x(fvGridGeometry->numDofs());
+    problem->applyInitialSolution(x);
+    auto xOld = x;
+
+    // the grid variables
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+    gridVariables->init(x);
+
+    // instantiate indicator & data transfer, read parameters for indicator
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
+    const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
+    TwoPGridAdaptIndicator<TypeTag> indicator(fvGridGeometry);
+    TwoPGridDataTransfer<TypeTag> dataTransfer(problem, fvGridGeometry, gridVariables, x);
+
+    // Do initial refinement around sources/BCs
+    GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, fvGridGeometry, gridVariables);
+
+    // refine up to the maximum level
+    const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
+    for (std::size_t i = 0; i < maxLevel; ++i)
+    {
+        initIndicator.calculate(x);
+
+        bool wasAdapted = false;
+        if (markElements(gridManager.grid(), initIndicator))
+            wasAdapted = adapt(gridManager.grid(), dataTransfer);
+
+        // update grid data after adaption
+        if (wasAdapted)
+        {
+            xOld = x; //!< Overwrite the old solution with the new (resized & interpolated) one
+            gridVariables->updateAfterGridAdaption(x); //!< Initialize the secondary variables to the new (and "new old") solution
+            problem->computePointSourceMap(); //!< Update the point source map
+        }
+    }
+
+    // Do refinement for the initial conditions using our indicator
+    indicator.calculate(x, refineTol, coarsenTol);
+
+    bool wasAdapted = false;
+    if (markElements(gridManager.grid(), indicator))
+        wasAdapted = adapt(gridManager.grid(), dataTransfer);
+
+    // update grid data after adaption
+    if (wasAdapted)
+    {
+        xOld = x; //!< Overwrite the old solution with the new (resized & interpolated) one
+        gridVariables->updateAfterGridAdaption(x); //!< Initialize the secondary variables to the new (and "new old") solution
+        problem->computePointSourceMap(); //!< Update the point source map
+    }
+
+    // get some time loop parameters
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+    // intialize the vtk output module
+    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+    using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
+    vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
+    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+    vtkWriter.write(0.0);
+
+    // instantiate time loop
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
+    timeLoop->setMaxTimeStepSize(maxDt);
+
+    // the assembler with time loop for instationary problem
+    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+    // the linear solver
+    using LinearSolver = AMGBackend<TypeTag>;
+    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    // time loop
+    timeLoop->start(); do
+    {
+        // refine/coarsen only after first time step
+        if (timeLoop->time() > 0)
+        {
+            // compute refinement indicator
+            indicator.calculate(x, refineTol, coarsenTol);
+
+            // mark elements and maybe adapt grid
+            wasAdapted = false;
+            if (markElements(gridManager.grid(), indicator))
+                wasAdapted = adapt(gridManager.grid(), dataTransfer);
+
+            if (wasAdapted)
+            {
+                // Note that if we were using point sources, we would have to update the map here as well
+                xOld = x; //!< Overwrite the old solution with the new (resized & interpolated) one
+                assembler->setJacobianPattern(); //!< Tell the assembler to resize the matrix and set pattern
+                assembler->setResidualSize(); //!< Tell the assembler to resize the residual
+                gridVariables->updateAfterGridAdaption(x); //!< Initialize the secondary variables to the new (and "new old") solution
+                problem->computePointSourceMap(); //!< Update the point source map
+            }
+        }
+
+        // set previous solution for storage evaluations
+        assembler->setPreviousSolution(xOld);
+
+        // solve the non-linear system with time step control
+        nonLinearSolver.solve(x, *timeLoop);
+
+        // make the new solution the old solution
+        xOld = x;
+        gridVariables->advanceTimeStep();
+
+        // advance to the time loop to the next step
+        timeLoop->advanceTimeStep();
+
+        // write vtk output
+        vtkWriter.write(timeLoop->time());
+
+        // report statistics of this time step
+        timeLoop->reportTimeStep();
+
+        // set new dt as suggested by the newton solver
+        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+    } while (!timeLoop->finished());
+
+    timeLoop->finalize(leafGridView.comm());
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (Dune::DGFException & e)
+{
+    std::cerr << "DGF exception thrown (" << e <<
+                 "). Most likely, the DGF file name is wrong "
+                 "or the DGF file is corrupted, "
+                 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                 << " ---> Abort!" << std::endl;
+    return 2;
+}
+catch (Dune::Exception &e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 3;
+}
+catch (...)
+{
+    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+    return 4;
+}

test/examples/2pinfiltration/params.input

+[TimeLoop]
+DtInitial = 100 # [s]
+TEnd = 100 # [s]
+
+[Grid]
+UpperRight = 6 4
+Cells = 48 32
+Refinement = 1
+CellType = Cube
+ClosureType = Green
+
+[SpatialParams]
+LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
+LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
+
+[Problem]
+Name = 2padaptive # name passed to the output routines
+
+[Adaptive]
+RefineAtDirichletBC = 0
+RefineAtFluxBC = 1
+MinLevel = 0
+MaxLevel = 2
+CoarsenTolerance = 1e-4
+RefineTolerance = 1e-4

test/examples/2pinfiltration/problem.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPTests
+ * \brief Soil contamination problem where DNAPL infiltrates a fully
+ *        water saturated medium.
+ */
+
+#ifndef DUMUX_LENSPROBLEM_POINTSOURCE_ADAPTIVE_HH
+#define DUMUX_LENSPROBLEM_POINTSOURCE_ADAPTIVE_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/material/components/trichloroethene.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/fluidsystems/2pimmiscible.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/porousmediumflow/2p/model.hh>
+#include <dumux/porousmediumflow/2p/incompressiblelocalresidual.hh>
+
+#include "spatialparams.hh"
+
+#include <dumux/io/container.hh>
+
+namespace Dumux {
+  // forward declarations
+  template<class TypeTag> class PointSourceTestProblem;
+
+  namespace Properties {
+  // Create new type tags
+  namespace TTag {
+  struct TwoPAdaptivePointSource { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };
+  } // end namespace TTag
+
+  //! Use non-conforming refinement
+  // #if HAVE_DUNE_ALUGRID TODO
+  // template<class TypeTag>
+  // struct Grid<TypeTag, TTag::TwoPAdaptivePointSource> { using type = Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>; };
+  // #else
+  template<class TypeTag>
+  struct Grid<TypeTag, TTag::TwoPAdaptivePointSource> { using type = Dune::YaspGrid<2>; };
+  // #endif
+
+  template<class TypeTag>
+  struct Problem<TypeTag, TTag::TwoPAdaptivePointSource> { using type = PointSourceTestProblem<TypeTag>; };
+
+  // the local residual containing the analytic derivative methods
+  template<class TypeTag>
+  struct LocalResidual<TypeTag, TTag::TwoPAdaptivePointSource> { using type = TwoPIncompressibleLocalResidual<TypeTag>; };
+
+  // Set the fluid system
+  template<class TypeTag>
+  struct FluidSystem<TypeTag, TTag::TwoPAdaptivePointSource>
+  {
+      using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+      using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
+      using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
+      using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
+  };
+
+  // Set the spatial parameters
+  template<class TypeTag>
+  struct SpatialParams<TypeTag, TTag::TwoPAdaptivePointSource>
+  {
+  private:
+      using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+      using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+  public:
+      using type = TwoPTestSpatialParams<FVGridGeometry, Scalar>;
+  };
+
+  // Enable caching
+  template<class TypeTag>
+  struct EnableGridVolumeVariablesCache<TypeTag, TTag::TwoPAdaptivePointSource> { static constexpr bool value = false; };
+  template<class TypeTag>
+  struct EnableGridFluxVariablesCache<TypeTag, TTag::TwoPAdaptivePointSource> { static constexpr bool value = false; };
+  template<class TypeTag>
+  struct EnableFVGridGeometryCache<TypeTag, TTag::TwoPAdaptivePointSource> { static constexpr bool value = false; };
+
+  } // end namespace Properties
+
+/*!
+ * \ingroup TwoPTests
+ * \brief Soil contamination problem where DNAPL infiltrates a fully
+ *        water saturated medium.
+ */
+template <class TypeTag >
+class PointSourceTestProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Vertex = typename GridView::template Codim<GridView::dimensionworld>::Entity;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using PointSource =  GetPropType<TypeTag, Properties::PointSource>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    enum {
+        pressureH2OIdx = Indices::pressureIdx,
+        saturationDNAPLIdx = Indices::saturationIdx,
+        contiDNAPLEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx,
+        waterPhaseIdx = FluidSystem::phase0Idx,
+        dnaplPhaseIdx = FluidSystem::phase1Idx
+    };
+
+public:
+  PointSourceTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+  : ParentType(fvGridGeometry)
+  {
+    initialValues_ = readFileToContainer<std::vector<PrimaryVariables>>("initialsolutioncc.txt");
+  }
+
+  /*!
+   * \brief Specifies which kind of boundary condition should be
+   *        used for which equation on a given boundary segment
+   *
+   * \param globalPos The global position
+   */
+  BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+  {
+    BoundaryTypes values;
+    if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+        values.setAllDirichlet();
+    else
+        values.setAllNeumann();
+    return values;
+  }
+
+      /*!
+       * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
+       *
+       * \param globalPos The global position
+       */
+      PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+      {
+          PrimaryVariables values;
+          GetPropType<TypeTag, Properties::FluidState> fluidState;
+          fluidState.setTemperature(temperature());
+          fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
+          fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
+
+          Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
+
+          Scalar height = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
+          Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
+          Scalar alpha = 1 + 1.5/height;
+          Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
+          Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
+
+          // hydrostatic pressure scaled by alpha
+          values[pressureH2OIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
+          values[saturationDNAPLIdx] = 0.0;
+
+          return values;
+      }
+
+      /*!
+       * \brief Evaluates the boundary conditions for a Neumann boundary segment.
+       *
+       * \param globalPos The position of the integration point of the boundary segment.
+       *
+       * For this method, the \a values parameter stores the mass flux
+       * in normal direction of each phase. Negative values mean influx.
+       */
+      NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
+      {
+          NumEqVector values(0.0);
+          if (onInlet_(globalPos))
+              values[contiDNAPLEqIdx] = -0.04; // kg / (m * s)
+
+          // in the test with the oil wet lens, use higher injection rate
+          if (this->spatialParams().lensIsOilWet())
+              values[contiDNAPLEqIdx] *= 10;
+
+          return values;
+      }
+
+  /*!
+   * \brief Evaluates the initial value for an element for cell-centered models.
+   */
+  PrimaryVariables initial(const Element& element) const
+  {
+      const auto delta = 0.0625;
+      unsigned int cellsX = this->fvGridGeometry().bBoxMax()[0]/delta;
+      const auto globalPos = element.geometry().center();
+
+      // the input data corresponds to a uniform grid with discretization length deltaX_
+      unsigned int dataIdx = std::trunc(globalPos[1]/delta) * cellsX + std::trunc(globalPos[0]/delta);
+      return initialValues_[dataIdx];
+  }
+
+  /*!
+   * \brief Evaluates the initial values for a control volume.
+   *
+   * \param globalPos The global position
+   */
+  PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+  {
+      PrimaryVariables values;
+      GetPropType<TypeTag, Properties::FluidState> fluidState;
+      fluidState.setTemperature(temperature());
+      fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
+      fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
+
+      Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
+
+      Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
+
+      // hydrostatic pressure
+      values[pressureH2OIdx] = 1e5 - densityW*this->gravity()[1]*depth;
+      values[saturationDNAPLIdx] = 0;
+      return values;
+  }
+
+  /*!
+   * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
+   *
+   * This is not specific to the discretization. By default it just
+   * throws an exception so it must be overloaded by the problem if
+   * no energy equation is used.
+   */
+  Scalar temperature() const
+  {
+      return 293.15; // 10°C
+  }
+
+
+    /*!
+     * \brief Applies a vector of point sources. The point sources
+     *        are possibly solution dependent.
+     *
+     * \param pointSources A vector of PointSources that contain
+              source values for all phases and space positions.
+     *
+     * For this method, the \a values method of the point source
+     * has to return the absolute mass rate in untis
+     * \f$ [ \textnormal{unit of conserved quantity} / s ] \f$.
+     * Positive values mean that mass is created, negative ones mean that it vanishes.
+     */
+    void addPointSources(std::vector<PointSource>& pointSources) const
+    {
+        // inject 2 kg/s of non-wetting phase at position (1, 1);
+        pointSources.push_back(PointSource({0.502, 3.02}, {0, 0.1}));
+    }
+
+  private:
+    bool onLeftBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;
+    }
+
+    bool onRightBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
+    }
+
+    bool onLowerBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_;
+    }
+
+    bool onUpperBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
+    }
+
+    bool onInlet_(const GlobalPosition &globalPos) const
+    {
+        Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
+        Scalar lambda = (this->fvGridGeometry().bBoxMax()[0] - globalPos[0])/width;
+        return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
+    }
+
+    static constexpr Scalar eps_ = 1e-6;
+    std::vector<PrimaryVariables> initialValues_;
+};
+
+} // end namespace Dumux
+
+#endif

test/examples/2pinfiltration/spatialparams.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTests
+ * \brief The spatial params for the incompressible 2p test.
+ */
+
+#ifndef DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+#include <dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TwoPTests
+ * \brief The spatial params for the incompressible 2p test.
+ */
+template<class FVGridGeometry, class Scalar>
+class TwoPTestSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar, TwoPTestSpatialParams<FVGridGeometry, Scalar>>
+{
+    using GridView = typename FVGridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using ThisType = TwoPTestSpatialParams<FVGridGeometry, Scalar>;
+    using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
+
+    static constexpr int dimWorld = GridView::dimensionworld;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+    using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
+
+public:
+    using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+    using MaterialLawParams = typename MaterialLaw::Params;
+    using PermeabilityType = Scalar;
+
+    TwoPTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        lensIsOilWet_ = getParam<bool>("SpatialParams.LensIsOilWet", false);
+
+        lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
+        lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
+
+        // residual saturations
+        lensMaterialParams_.setSwr(0.18);
+        lensMaterialParams_.setSnr(0.0);
+        outerMaterialParams_.setSwr(0.05);
+        outerMaterialParams_.setSnr(0.0);
+
+        // parameters for the Van Genuchten law
+        // alpha and n
+        lensMaterialParams_.setVgAlpha(0.00045);
+        lensMaterialParams_.setVgn(7.3);
+        outerMaterialParams_.setVgAlpha(0.0037);
+        outerMaterialParams_.setVgn(4.7);
+
+        lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
+        outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
+    }
+
+    /*!
+     * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+     *        In this test, we use element-wise distributed permeabilities.
+     *
+     * \param element The current element
+     * \param scv The sub-control volume inside the element.
+     * \param elemSol The solution at the dofs connected to the element.
+     * \return The permeability
+     */
+    template<class ElementSolution>
+    PermeabilityType permeability(const Element& element,
+                                  const SubControlVolume& scv,
+                                  const ElementSolution& elemSol) const
+    {
+
+        // do not use a less permeable lens in the test with inverted wettability
+        if (isInLens_(element.geometry().center()) && !lensIsOilWet_)
+            return lensK_;
+        return outerK_;
+    }
+
+    /*!
+     * \brief Returns the porosity \f$[-]\f$
+     *
+     * \param globalPos The global position
+     */
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
+    { return 0.4; }
+
+    /*!
+     * \brief Returns the parameter object for the Brooks-Corey material law.
+     *
+     * In this test, we use element-wise distributed material parameters.
+     *
+     * \param element The current element
+     * \param scv The sub-control volume inside the element.
+     * \param elemSol The solution at the dofs connected to the element.
+     * \return The material parameters object
+     */
+    template<class ElementSolution>
+    const MaterialLawParams& materialLawParams(const Element& element,
+                                               const SubControlVolume& scv,
+                                               const ElementSolution& elemSol) const
+    {
+        // do not use different parameters in the test with inverted wettability
+        if (isInLens_(element.geometry().center()) && !lensIsOilWet_)
+            return lensMaterialParams_;
+        return outerMaterialParams_;
+    }
+
+    /*!
+     * \brief Function for defining which phase is to be considered as the wetting phase.
+     *
+     * \param globalPos The global position
+     * \return The wetting phase index
+     */
+    template<class FluidSystem>
+    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+    {
+        if (isInLens_(globalPos) && lensIsOilWet_)
+            return FluidSystem::phase1Idx;
+        return FluidSystem::phase0Idx;
+    }
+
+    //! Updates the map of which material parameters are associated with a nodal dof.
+    template<class SolutionVector>
+    void updateMaterialInterfaceParams(const SolutionVector& x)
+    {
+        if (FVGridGeometry::discMethod == DiscretizationMethod::box)
+            materialInterfaceParams_.update(this->fvGridGeometry(), *this, x);
+    }
+
+    //! Returns the material parameters associated with a nodal dof
+    const BoxMaterialInterfaceParams<ThisType>& materialInterfaceParams() const
+    { return materialInterfaceParams_; }
+
+    //! Returns whether or not the lens is oil wet
+    bool lensIsOilWet() const { return lensIsOilWet_; }
+
+private:
+    bool isInLens_(const GlobalPosition &globalPos) const
+    {
+        for (int i = 0; i < dimWorld; ++i) {
+            if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
+                return false;
+        }
+        return true;
+    }
+
+    bool lensIsOilWet_;
+    GlobalPosition lensLowerLeft_;
+    GlobalPosition lensUpperRight_;
+
+    Scalar lensK_;
+    Scalar outerK_;
+    MaterialLawParams lensMaterialParams_;
+    MaterialLawParams outerMaterialParams_;
+
+    // Determines the parameters associated with the dofs at material interfaces
+    BoxMaterialInterfaceParams<ThisType> materialInterfaceParams_;
+
+    static constexpr Scalar eps_ = 1.5e-7;
+};
+
+} // end namespace Dumux
+
+#endif

test/examples/CMakeLists.txt

+add_subdirectory(2pinfiltration)