[pmflow] add incompressible 2p test with analytic jacobian

dumux/porousmediumflow/2p/implicit/CMakeLists.txt

@@ -3,9 +3,11 @@
 install(FILES
 gridadaptindicator.hh
 indices.hh
+incompressiblelocalresidual.hh
 localresidual.hh
 model.hh
 properties.hh
 propertydefaults.hh
 volumevariables.hh
+vtkoutputfields.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/2p/implicit)

dumux/porousmediumflow/2p/implicit/incompressiblelocalresidual.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 2 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
+ *
+ * \brief Element-wise calculation of the residual and its derivatives
+ *        for a two-phase, incompressible test problem.
+ */
+#ifndef DUMUX_2P_INCOMPRESSIBLE_TEST_LOCAL_RESIDUAL_HH
+#define DUMUX_2P_INCOMPRESSIBLE_TEST_LOCAL_RESIDUAL_HH
+
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class TwoPIncompressibleLocalResidual : public ImmiscibleLocalResidual<TypeTag>
+{
+    using ParentType = ImmiscibleLocalResidual<TypeTag>;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    // first index for the mass balance
+    enum
+    {
+        contiWEqIdx = Indices::contiWEqIdx,
+        contiNEqIdx = Indices::contiNEqIdx,
+
+        pressureIdx = Indices::pressureIdx,
+        saturationIdx = Indices::saturationIdx
+    };
+
+    static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+
+public:
+    using ParentType::ParentType;
+
+    template<class PartialDerivativeMatrix>
+    void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
+                               const Problem& problem,
+                               const Element& element,
+                               const FVElementGeometry& fvGeometry,
+                               const VolumeVariables& curVolVars) const
+    {
+        // we know that these values are constant throughout the simulation
+        static const auto phi = curVolVars.porosity();
+        static const auto phi_rho_w = phi*curVolVars.density(FluidSystem::wPhaseIdx);
+        static const auto phi_rho_n = phi*curVolVars.density(FluidSystem::nPhaseIdx);
+
+        const auto volume = element.geometry().volume();
+
+        // partial derivative of wetting phase storage term w.r.t. p_w
+        partialDerivatives[contiWEqIdx][pressureIdx] += 0.0;
+        // partial derivative of wetting phase storage term w.r.t. S_n
+        partialDerivatives[contiWEqIdx][saturationIdx] -= volume*phi_rho_w/this->timeLoop().timeStepSize();
+        // partial derivative of non-wetting phase storage term w.r.t. p_w
+        partialDerivatives[contiNEqIdx][pressureIdx] += 0.0;
+        // partial derivative of non-wetting phase storage term w.r.t. S_n
+        partialDerivatives[contiNEqIdx][saturationIdx] += volume*phi_rho_n/this->timeLoop().timeStepSize();
+    }
+
+    template<class PartialDerivativeMatrix>
+    void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
+                              const Problem& problem,
+                              const Element& element,
+                              const FVElementGeometry& fvGeometry,
+                              const VolumeVariables& curVolVars) const {}
+
+    template<class PartialDerivativeMatrices>
+    void addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
+                            const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& curElemVolVars,
+                            const ElementFluxVariablesCache& elemFluxVarsCache,
+                            const SubControlVolumeFace& scvf) const
+    {
+        using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+        using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+        static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight);
+
+        // evaluate the current wetting phase Darcy flux and resulting upwind weights
+        const auto wettingflux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
+                                                     scvf, FluidSystem::wPhaseIdx, elemFluxVarsCache);
+        const auto insideUpwindWeightWetting = std::signbit(wettingflux) ? (1.0 - upwindWeight) : upwindWeight;
+        const auto outsideUpwindWeightWetting = 1.0 - insideUpwindWeightWetting;
+
+        const auto nonwettingflux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
+                                                        scvf, FluidSystem::nPhaseIdx, elemFluxVarsCache);
+        const auto insideUpwindWeightNonWetting = std::signbit(nonwettingflux) ? (1.0 - upwindWeight) : upwindWeight;
+        const auto outsideUpwindWeightNonWetting = 1.0 - insideUpwindWeightNonWetting;
+
+        // get references to the two participating vol vars & partial derivative matrices
+        const auto& insideVolVars = curElemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
+        auto& dI_dJ_inside = derivativeMatrices[scvf.insideScvIdx()];
+        auto& dI_dJ_outside = derivativeMatrices[scvf.outsideScvIdx()];
+
+        // some quantities to be reused (rho & mu are constant and thus equal for all cells)
+        static const auto rho_w = insideVolVars.density(FluidSystem::wPhaseIdx);
+        static const auto rho_n = insideVolVars.density(FluidSystem::nPhaseIdx);
+        static const auto rhow_muw = rho_w/insideVolVars.viscosity(FluidSystem::wPhaseIdx);
+        static const auto rhon_mun = rho_n/insideVolVars.viscosity(FluidSystem::nPhaseIdx);
+        const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(FluidSystem::wPhaseIdx);
+        const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(FluidSystem::nPhaseIdx);
+        const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(FluidSystem::wPhaseIdx);
+        const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(FluidSystem::nPhaseIdx);
+
+        // we know the material parameters are only position-dependent
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+        const auto& insideMaterialParams = problem.spatialParams().materialLawParamsAtPos(insideScv.center());
+        const auto& outsideMaterialParams = problem.spatialParams().materialLawParamsAtPos(outsideScv.center());
+
+        // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
+        const auto insideSw = insideVolVars.saturation(FluidSystem::wPhaseIdx);
+        const auto outsideSw = outsideVolVars.saturation(FluidSystem::wPhaseIdx);
+
+        const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
+        const auto dKrw_dSn_outside = -1.0*MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
+        const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
+        const auto dKrn_dSn_outside = -1.0*MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
+
+        const auto dpc_dSn_inside = -1.0*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
+        const auto dpc_dSn_outside = -1.0*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
+
+        const auto tij = elemFluxVarsCache[scvf].advectionTij();
+        // partial derivative of the wetting phase flux w.r.t. p_w
+        dI_dJ_inside[contiWEqIdx][pressureIdx] += tij*(rhowKrw_muw_inside*insideUpwindWeightWetting
+                                                       + rhowKrw_muw_outside*outsideUpwindWeightWetting);
+        dI_dJ_outside[contiWEqIdx][pressureIdx] -= tij*(rhowKrw_muw_inside*insideUpwindWeightWetting
+                                                        + rhowKrw_muw_outside*outsideUpwindWeightWetting);
+
+        // partial derivative of the wetting phase flux w.r.t. S_n
+        dI_dJ_inside[contiWEqIdx][saturationIdx] += rhow_muw*wettingflux*dKrw_dSn_inside*insideUpwindWeightWetting;
+        dI_dJ_outside[contiWEqIdx][saturationIdx] += rhow_muw*wettingflux*dKrw_dSn_outside*outsideUpwindWeightWetting;
+
+        // partial derivative of the non-wetting phase flux w.r.t. p_w
+        dI_dJ_inside[contiNEqIdx][pressureIdx] += tij*(rhonKrn_mun_inside*insideUpwindWeightNonWetting
+                                                       + rhonKrn_mun_outside*outsideUpwindWeightNonWetting);
+        dI_dJ_outside[contiNEqIdx][pressureIdx] -= tij*(rhonKrn_mun_inside*insideUpwindWeightNonWetting
+                                                        + rhonKrn_mun_outside*outsideUpwindWeightNonWetting);
+
+        // partial derivative of the non-wetting phase flux w.r.t. S_n (relative permeability derivative contribution)
+        dI_dJ_inside[contiNEqIdx][saturationIdx] += rhon_mun*nonwettingflux*dKrn_dSn_inside*insideUpwindWeightNonWetting;
+        dI_dJ_outside[contiNEqIdx][saturationIdx] += rhon_mun*nonwettingflux*dKrn_dSn_outside*outsideUpwindWeightNonWetting;
+
+        // partial derivative of the non-wetting phase flux w.r.t. S_n (capillary pressure derivative contribution)
+        dI_dJ_inside[contiNEqIdx][saturationIdx] += tij*dpc_dSn_inside*(rhonKrn_mun_inside*insideUpwindWeightNonWetting
+                                                                        + rhonKrn_mun_outside*outsideUpwindWeightNonWetting);
+        dI_dJ_outside[contiNEqIdx][saturationIdx] -= tij*dpc_dSn_outside*(rhonKrn_mun_inside*insideUpwindWeightNonWetting
+                                                                          + rhonKrn_mun_outside*outsideUpwindWeightNonWetting);
+    }
+
+    template<class PartialDerivativeMatrices>
+    void addDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
+                                     const Problem& problem,
+                                     const Element& element,
+                                     const FVElementGeometry& fvGeometry,
+                                     const ElementVolumeVariables& curElemVolVars,
+                                     const ElementFluxVariablesCache& elemFluxVarsCache,
+                                     const SubControlVolumeFace& scvf) const
+    {
+        using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+        using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+        static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight);
+
+        // evaluate the current wetting phase Darcy flux and resulting upwind weights
+        const auto wettingflux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
+                                                     scvf, FluidSystem::wPhaseIdx, elemFluxVarsCache);
+        const auto insideUpwindWeightWetting = std::signbit(wettingflux) ? (1.0 - upwindWeight) : upwindWeight;
+        const auto outsideUpwindWeightWetting = 1.0 - insideUpwindWeightWetting;
+
+        const auto nonwettingflux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
+                                                        scvf, FluidSystem::nPhaseIdx, elemFluxVarsCache);
+        const auto insideUpwindWeightNonWetting = std::signbit(nonwettingflux) ? (1.0 - upwindWeight) : upwindWeight;
+        const auto outsideUpwindWeightNonWetting = 1.0 - insideUpwindWeightNonWetting;
+
+        // get references to the two participating vol vars & partial derivative matrices
+        const auto& insideVolVars = curElemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
+        auto& dI_dJ_inside = derivativeMatrices[scvf.insideScvIdx()];
+
+        // some quantities to be reused (rho & mu are constant and thus equal for all cells)
+        static const auto rho_w = insideVolVars.density(FluidSystem::wPhaseIdx);
+        static const auto rho_n = insideVolVars.density(FluidSystem::nPhaseIdx);
+        static const auto rhow_muw = rho_w/insideVolVars.viscosity(FluidSystem::wPhaseIdx);
+        static const auto rhon_mun = rho_n/insideVolVars.viscosity(FluidSystem::nPhaseIdx);
+        const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(FluidSystem::wPhaseIdx);
+        const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(FluidSystem::nPhaseIdx);
+        const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(FluidSystem::wPhaseIdx);
+        const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(FluidSystem::nPhaseIdx);
+
+        // we know the material parameters are only position-dependent
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& insideMaterialParams = problem.spatialParams().materialLawParamsAtPos(insideScv.center());
+
+        // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
+        const auto insideSw = insideVolVars.saturation(FluidSystem::wPhaseIdx);
+        const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
+        const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
+        const auto dpc_dSn_inside = -1.0*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
+
+        const auto tij = elemFluxVarsCache[scvf].advectionTij();
+        // partial derivative of the wetting phase flux w.r.t. p_w
+        dI_dJ_inside[contiWEqIdx][pressureIdx] += tij*(rhowKrw_muw_inside*insideUpwindWeightWetting
+                                                       + rhowKrw_muw_outside*outsideUpwindWeightWetting);
+
+        // partial derivative of the wetting phase flux w.r.t. S_n
+        dI_dJ_inside[contiWEqIdx][saturationIdx] += rhow_muw*wettingflux*dKrw_dSn_inside*insideUpwindWeightWetting;
+
+        // partial derivative of the non-wetting phase flux w.r.t. p_w
+        dI_dJ_inside[contiNEqIdx][pressureIdx] += tij*(rhonKrn_mun_inside*insideUpwindWeightNonWetting
+                                                       + rhonKrn_mun_outside*outsideUpwindWeightNonWetting);
+
+        // partial derivative of the non-wetting phase flux w.r.t. S_n (relative permeability derivative contribution)
+        dI_dJ_inside[contiNEqIdx][saturationIdx] += rhon_mun*nonwettingflux*dKrn_dSn_inside*insideUpwindWeightNonWetting;
+
+        // partial derivative of the non-wetting phase flux w.r.t. S_n (capillary pressure derivative contribution)
+        dI_dJ_inside[contiNEqIdx][saturationIdx] += tij*dpc_dSn_inside*(rhonKrn_mun_inside*insideUpwindWeightNonWetting
+                                                                        + rhonKrn_mun_outside*outsideUpwindWeightNonWetting);
+    }
+
+    template<class PartialDerivativeMatrices>
+    void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
+                                 const Problem& problem,
+                                 const Element& element,
+                                 const FVElementGeometry& fvGeometry,
+                                 const ElementVolumeVariables& curElemVolVars,
+                                 const ElementFluxVariablesCache& elemFluxVarsCache,
+                                 const SubControlVolumeFace& scvf) const
+    {
+        // we have no Robin-type boundary conditions here. Do nothing.
+    }
+};
+
+} // end namespace Dumux
+
+#endif

dumux/porousmediumflow/2p/implicit/propertydefaults.hh

@@ -32,6 +32,7 @@
 #include "model.hh"
 #include "indices.hh"
 #include "volumevariables.hh"
+#include "vtkoutputfields.hh"
 
 #include <dumux/porousmediumflow/immiscible/localresidual.hh>
 #include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
@@ -93,6 +94,9 @@ SET_TYPE_PROP(TwoP,
 //! Use ImplicitSpatialParams by default.
 SET_TYPE_PROP(TwoP, SpatialParams, ImplicitSpatialParams<TypeTag>);
 
+//! Set the vtk output fields specific to the twop model
+SET_TYPE_PROP(TwoP, VtkOutputFields, TwoPVtkOutputFields<TypeTag>);
+
 /*!
  * \brief Set the property for the material parameters by extracting
  *        it from the material law.

dumux/porousmediumflow/2p/implicit/vtkoutputfields.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 2 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
+ * \brief Adds vtk output fields specific to the twop model
+ */
+#ifndef DUMUX_TWOP_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_TWOP_VTK_OUTPUT_FIELDS_HH
+
+#include <dumux/implicit/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoP, InputOutput
+ * \brief Adds vtk output fields specific to the twop model
+ */
+template<class TypeTag>
+class TwoPVtkOutputFields
+{
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+public:
+    template <class VtkOutputModule>
+    static void init(VtkOutputModule& vtk)
+    {
+        vtk.addSecondaryVariable("Sw", [](const VolumeVariables& v){ return v.saturation(Indices::wPhaseIdx); });
+        vtk.addSecondaryVariable("Sn", [](const VolumeVariables& v){ return v.saturation(Indices::nPhaseIdx); });
+        vtk.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(Indices::wPhaseIdx); });
+        vtk.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(Indices::nPhaseIdx); });
+        vtk.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
+        vtk.addSecondaryVariable("rhoW", [](const VolumeVariables& v){ return v.density(Indices::wPhaseIdx); });
+        vtk.addSecondaryVariable("rhoN", [](const VolumeVariables& v){ return v.density(Indices::nPhaseIdx); });
+        vtk.addSecondaryVariable("mobW", [](const VolumeVariables& v){ return v.mobility(Indices::wPhaseIdx); });
+        vtk.addSecondaryVariable("mobN", [](const VolumeVariables& v){ return v.mobility(Indices::nPhaseIdx); });
+        vtk.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+        vtk.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
+    }
+};
+
+} // end namespace Dumux
+
+#endif

test/porousmediumflow/2p/implicit/CMakeLists.txt

 add_subdirectory(pointsources)
+add_subdirectory(incompressible)
 add_input_file_links()
 
 # isothermal tests

test/porousmediumflow/2p/implicit/incompressible/CMakeLists.txt

+dune_symlink_to_source_files(FILES "test_2pincompressible.input")
+
+# incompressible
+add_dumux_test(test_2pincompressible test_2pincompressible test_2pincompressible.cc
+               ${CMAKE_CURRENT_BINARY_DIR}/test_2pincompressible)
+
+set(CMAKE_BUILD_TYPE Release)
\ No newline at end of file

test/porousmediumflow/2p/implicit/incompressible/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 2 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
+ * \brief The properties for the incompressible test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
+
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/components/dnapl.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/2p/implicit/propertydefaults.hh>
+#include <dumux/porousmediumflow/2p/implicit/incompressiblelocalresidual.hh>
+
+#include "spatialparams.hh"
+
+namespace Dumux
+{
+// forward declarations
+template<class TypeTag> class TwoPTestProblem;
+
+namespace Properties
+{
+
+NEW_PROP_TAG(EnableFVGridGeometryCache);
+NEW_PROP_TAG(FVGridGeometry);
+
+NEW_TYPE_TAG(IncompressibleTestProblem, INHERITS_FROM(CCTpfaModel, TwoP, SpatialParams));
+
+// Set the grid type
+SET_TYPE_PROP(IncompressibleTestProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem type
+SET_TYPE_PROP(IncompressibleTestProblem, Problem, TwoPTestProblem<TypeTag>);
+
+// the local residual containing the analytic derivative methods
+SET_TYPE_PROP(IncompressibleTestProblem, LocalResidual, TwoPIncompressibleLocalResidual<TypeTag>);
+
+// Set the wetting phase
+SET_PROP(IncompressibleTestProblem, WettingPhase)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+    typedef FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(IncompressibleTestProblem, NonwettingPhase)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+    typedef FluidSystems::LiquidPhase<Scalar, DNAPL<Scalar> > type;
+};
+
+// Enable caching
+SET_BOOL_PROP(IncompressibleTestProblem, EnableGlobalVolumeVariablesCache, false);
+SET_BOOL_PROP(IncompressibleTestProblem, EnableGlobalFluxVariablesCache, false);
+SET_BOOL_PROP(IncompressibleTestProblem, EnableFVGridGeometryCache, false);
+} // end namespace Properties
+
+template<class TypeTag>
+class TwoPTestProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
+    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    static constexpr int dimWorld = GridView::dimensionworld;
+
+public:
+    TwoPTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry) {}
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment
+     *
+     * \param values Stores the value of the boundary type
+     * \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 values Stores the Dirichlet values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} ] \f$
+     * \param globalPos The global position
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+    {
+        PrimaryVariables values;
+        typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
+        fluidState.setTemperature(temperature());
+        fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
+        fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
+
+        Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
+
+        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[Indices::pwIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
+        values[Indices::snIdx] = 0.0;
+
+        return values;
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a neumann
+     *        boundary segment.
+     *
+     * \param values Stores the Neumann values for the conservation equations in
+     *               \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
+     * \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.
+     */
+    NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
+    {
+        NeumannFluxes values(0.0);
+        if (onInlet_(globalPos))
+            values[Indices::contiNEqIdx] = -0.04; // kg / (m * s)
+        return values;
+    }
+
+    /*!
+     * \brief Evaluates the initial values for a control volume
+     *
+     * \param values Stores the initial values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variables} ] \f$
+     * \param globalPos The global position
+     */
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    {
+        PrimaryVariables values;
+        typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
+        fluidState.setTemperature(temperature());
+        fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
+        fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
+
+        Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
+
+        Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
+
+        // hydrostatic pressure
+        values[Indices::pwIdx] = 1e5 - densityW*this->gravity()[1]*depth;
+        values[Indices::snIdx] = 0.0;
+        return values;