diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/CMakeLists.txt b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b385e900529a7af019a12c75f330ee13256940b1
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/CMakeLists.txt
@@ -0,0 +1,25 @@
+add_input_file_links()
+
+dune_add_test(NAME test_mpnc_comparison_box
+ SOURCES test_mpnc_comparison_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=MPNCComparisonBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/mpnc_2p2c_box-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_mpnc_comparison_box-00009.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_mpnc_comparison_box test_mpnc_comparison_fv.input -Problem.Name test_mpnc_comparison_box")
+
+dune_add_test(NAME test_mpnc_comparison_tpfa
+ SOURCES test_mpnc_comparison_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=MPNCComparisonCCProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/mpnc_2p2c_tpfa-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_mpnc_comparison_tpfa-00009.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_mpnc_comparison_tpfa test_mpnc_comparison_fv.input -Problem.Name test_mpnc_comparison_tpfa")
+
+#install sources
+install(FILES
+mpnc_comparison_problem.hh
+test_mpnc_comparison_fv.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/mpnc/2p2ccomparison)
diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/mpnc_comparison_problem.hh b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/mpnc_comparison_problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e8b55df56b3b88bcfb702fdb9f9fd644f542d5db
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/mpnc_comparison_problem.hh
@@ -0,0 +1,355 @@
+// -*- 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
+ * \ingroup MPNCTests
+ * \brief Problem where air is injected in a unsaturated porous medium. This test compares a mpnc problem with a 2p2c problem
+ */
+#ifndef DUMUX_MPNC_TWOPTWOC_COMPARISON_OBSTACLEPROBLEM_HH
+#define DUMUX_MPNC_TWOPTWOC_COMPARISON_OBSTACLEPROBLEM_HH
+
+#include <dune/common/parametertreeparser.hh>
+
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/mpnc/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <test/porousmediumflow/2p2c/implicit/comparison/vtkoutputfields.hh>
+
+#include <dumux/material/fluidsystems/h2on2.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+
+#include "mpnc_comparison_spatialparams.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup MPNCTests
+ * \brief Problem where air is injected in a unsaturated porous medium. This test compares a mpnc problem with a 2p2c problem
+ */
+template <class TypeTag>
+class MPNCComparisonProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(MPNCComparisonProblem, INHERITS_FROM(MPNC, MPNCComparisonSpatialParams));
+NEW_TYPE_TAG(MPNCComparisonBoxProblem, INHERITS_FROM(BoxModel, MPNCComparisonProblem));
+NEW_TYPE_TAG(MPNCComparisonCCProblem, INHERITS_FROM(CCTpfaModel, MPNCComparisonProblem));
+
+// Set the grid type
+SET_TYPE_PROP(MPNCComparisonProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(MPNCComparisonProblem,
+ Problem,
+ MPNCComparisonProblem<TypeTag>);
+
+// Set fluid configuration
+SET_TYPE_PROP(MPNCComparisonProblem,
+ FluidSystem,
+ FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/false>);
+
+// decide which type to use for floating values (double / quad)
+SET_TYPE_PROP(MPNCComparisonProblem, Scalar, double);
+
+SET_BOOL_PROP(MPNCComparisonProblem, EnableMolecularDiffusion, true);
+
+SET_BOOL_PROP(MPNCComparisonProblem, UseMoles, true);
+SET_TYPE_PROP(MPNCComparisonProblem, VtkOutputFields, TwoPTwoCMPNCVtkOutputFields<TypeTag>);
+
+}
+
+
+/*!
+ * \ingroup MPNCTests
+ * \briefProblem where air is injected in a unsaturated porous medium. This test compares a mpnc problem with a 2p2c problem
+ *
+ */
+template <class TypeTag>
+class MPNCComparisonProblem
+ : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using ParameterCache = typename FluidSystem::ParameterCache;
+
+ // world dimension
+ enum {dimWorld = GridView::dimensionworld};
+ enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
+ enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
+ enum {nPhaseIdx = FluidSystem::nPhaseIdx};
+ enum {wPhaseIdx = FluidSystem::wPhaseIdx};
+ enum {wCompIdx = FluidSystem::wCompIdx};
+ enum {nCompIdx = FluidSystem::nCompIdx};
+ enum {fug0Idx = Indices::fug0Idx};
+ enum {s0Idx = Indices::s0Idx};
+ enum {p0Idx = Indices::p0Idx};
+
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ */
+ MPNCComparisonProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ temperature_ = 273.15 + 25; // -> 25°C
+
+ // initialize the tables of the fluid system
+ Scalar Tmin = temperature_ - 1.0;
+ Scalar Tmax = temperature_ + 1.0;
+ int nT = 3;
+
+ Scalar pmin = 1.0e5 * 0.75;
+ Scalar pmax = 2.0e5 * 1.25;
+ int np = 1000;
+
+ FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
+ name_ = getParam<std::string>("Problem.Name");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the problem name
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string name() const
+ { return name_; }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ *
+ * \param globalPos The global position
+ */
+ Scalar temperature() const
+ { return temperature_; }
+
+ /*!
+ * \brief Write a restart file?
+ */
+ bool shouldWriteRestartFile() const
+ {
+ return ParentType::shouldWriteRestartFile();
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+ /*!
+ * \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 bcTypes;
+ if (onOutlet_(globalPos))
+ bcTypes.setAllDirichlet();
+ else
+ bcTypes.setAllNeumann();
+ return bcTypes;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet
+ * boundary segment
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ return initial_(globalPos);
+ }
+
+ /*!
+ * \brief Evaluates 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 element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param intersection The intersection between element and boundary
+ * \param scvIdx The local index of the sub-control volume
+ * \param boundaryFaceIdx The index of the boundary face
+ *
+ * Negative values mean influx.
+ */
+ PrimaryVariables neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ const auto& globalPos = scvf.ipGlobal();
+ Scalar injectedAirMass = -1e-3;
+ Scalar injectedAirMolarMass = injectedAirMass/FluidSystem::molarMass(nCompIdx);
+ if (onInlet_(globalPos))
+ {
+ values[Indices::conti0EqIdx+1] = injectedAirMolarMass;
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all balance equations within a given
+ * sub-control-volume.
+ *
+ * \param values Stores the solution for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ *
+ * Positive values mean that mass is created, negative ones mean that it vanishes.
+ */
+ //! \copydoc Dumux::ImplicitProblem::source()
+ PrimaryVariables source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param values The initial values for the primary variables
+ * \param globalPos The center of the finite volume which ought to be set.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ return initial_(globalPos);
+ }
+
+ // \}
+
+private:
+ // the internal method for the initial condition
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ FluidState fs;
+
+ // set the fluid temperatures
+ fs.setTemperature(this->temperatureAtPos(globalPos));
+
+ // set water saturation
+ fs.setSaturation(wPhaseIdx, 0.8);
+ fs.setSaturation(nPhaseIdx, 1.0 - fs.saturation(wPhaseIdx));
+ // set pressure of the gas phase
+ fs.setPressure(nPhaseIdx, 1e5);
+ // calulate the capillary pressure
+ const auto& matParams =
+ this->spatialParams().materialLawParamsAtPos(globalPos);
+ PhaseVector pc;
+ MaterialLaw::capillaryPressures(pc, matParams, fs);
+ fs.setPressure(wPhaseIdx,
+ fs.pressure(nPhaseIdx) + pc[wPhaseIdx] - pc[nPhaseIdx]);
+
+ // make the fluid state consistent with local thermodynamic
+ // equilibrium
+ using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem, /*useKelvinEquation*/false>;
+
+ ParameterCache paramCache;
+ MiscibleMultiPhaseComposition::solve(fs,
+ paramCache,
+ /*setViscosity=*/true,
+ /*setEnthalpy=*/false);
+ ///////////
+ // assign the primary variables
+ ///////////
+
+ // all N component fugacities
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ values[fug0Idx + compIdx] = fs.fugacity(nPhaseIdx, compIdx);
+
+ // first M - 1 saturations
+ for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
+ values[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
+
+ // first pressure
+ values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
+ return values;
+ }
+
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar x = globalPos[0];
+ Scalar y = globalPos[1];
+ return x >= 60 - eps_ && y <= 10 + eps_;
+ }
+
+ bool onOutlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar x = globalPos[0];
+ Scalar y = globalPos[1];
+ return x < eps_ && y <= 10 + eps_;
+ }
+
+ Scalar temperature_;
+ static constexpr Scalar eps_ = 1e-6;
+ std::string name_;
+};
+} //end namespace
+
+#endif
diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/mpnc_comparison_spatialparams.hh b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/mpnc_comparison_spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1b1fcb737e2902991f61c2748c2767da61d7b1ea
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/mpnc_comparison_spatialparams.hh
@@ -0,0 +1,181 @@
+// -*- 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
+ * \ingroup TwoPTwoCTests
+ * \brief The spatial parameters for the TwoPTwoC MPNC comparison problem
+ */
+
+#ifndef DUMUX_MPNC_COMPARISON_SPATIAL_PARAMS_HH
+#define DUMUX_MPNC_COMPARISON_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+#include <dumux/material/fluidmatrixinteractions/mp/2padapter.hh>
+
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+namespace Dumux
+{
+/*!
+ * \ingroup MPNCTests
+ * \brief The spatial parameters for the ObstacleProblem
+ */
+//forward declaration
+template<class TypeTag>
+class MPNCComparisonSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(MPNCComparisonSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(MPNCComparisonSpatialParams, SpatialParams, MPNCComparisonSpatialParams<TypeTag>);
+
+// Set the material Law
+SET_PROP(MPNCComparisonSpatialParams, MaterialLaw)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ enum {wPhaseIdx = FluidSystem::wPhaseIdx};
+ // define the material law
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using EffMaterialLaw = RegularizedBrooksCorey<Scalar>;
+ using TwoPMaterialLaw = EffToAbsLaw<EffMaterialLaw>;
+public:
+ using type = TwoPAdapter<wPhaseIdx, TwoPMaterialLaw>;
+};
+}
+
+/**
+ * \ingroup MPNCModel
+ * \ingroup ImplicitTestProblems
+ * \brief Definition of the spatial params properties for the obstacle problem
+ *
+ */
+template<class TypeTag>
+class MPNCComparisonSpatialParams : public FVSpatialParams<TypeTag>
+{
+ using ParentType = FVSpatialParams<TypeTag>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
+
+ enum {dimWorld=GridView::dimensionworld};
+
+public:
+ using PermeabilityType = Scalar;
+
+
+ MPNCComparisonSpatialParams(const Problem &problem)
+ : ParentType(problem)
+ {
+ // intrinsic permeabilities
+ coarseK_ = 1e-12;
+ fineK_ = 1e-15;
+
+ // the porosity
+ porosity_ = 0.3;
+
+ // residual saturations
+ fineMaterialParams_.setSwr(0.2);
+ fineMaterialParams_.setSnr(0.0);
+ coarseMaterialParams_.setSwr(0.2);
+ coarseMaterialParams_.setSnr(0.0);
+
+ // parameters for the Brooks-Corey law
+ fineMaterialParams_.setPe(1e4);
+ coarseMaterialParams_.setPe(1e4);
+ fineMaterialParams_.setLambda(2.0);
+ coarseMaterialParams_.setLambda(2.0);
+ }
+
+
+
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ {
+ if (isFineMaterial_(scv.dofPosition()))
+ return fineK_;
+ else
+ return coarseK_;
+ }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the soil
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the porosity needs to be defined
+ */
+ Scalar porosity(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
+ {
+ return porosity_;
+ }
+
+ /*!
+ * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
+ *
+ * \param pos The global position of the sub-control volume.
+ * \return the material parameters object
+ */
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ {
+ if (isFineMaterial_(globalPos))
+ return fineMaterialParams_;
+ else
+ return coarseMaterialParams_;
+ }
+
+private:
+ /*!
+ * \brief Returns whether a given global position is in the
+ * fine-permeability region or not.
+ */
+ static bool isFineMaterial_(const GlobalPosition &pos)
+ {
+ return
+ 30 - eps_ <= pos[0] && pos[0] <= 50 + eps_ &&
+ 20 - eps_ <= pos[1] && pos[1] <= 40 + eps_;
+ }
+
+ Scalar coarseK_;
+ Scalar fineK_;
+ Scalar porosity_;
+ MaterialLawParams fineMaterialParams_;
+ MaterialLawParams coarseMaterialParams_;
+ static constexpr Scalar eps_ = 1e-6;
+};
+
+}
+
+#endif
diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/test_mpnc_comparison_fv.cc b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/test_mpnc_comparison_fv.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6b64b18c152a171e708da08da24478feca43baf9
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/test_mpnc_comparison_fv.cc
@@ -0,0 +1,252 @@
+// -*- 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 test for the mpnc porousmedium box 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 "mpnc_comparison_problem.hh"
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.LowerLeft Lower left corner coordinates\n"
+ "\t-Grid.UpperRight Upper right corner coordinates\n"
+ "\t-Grid.Cells Number of cells in respective coordinate directions\n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(TYPETAG);
+
+ // 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, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, 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 NewtonController = Dumux::NewtonController<Scalar>;
+ using NewtonMethod = NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+
+ // 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 newton controller
+ timeLoop->setTimeStepSize(newtonController->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;
+}
diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/test_mpnc_comparison_fv.input b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/test_mpnc_comparison_fv.input
new file mode 100644
index 0000000000000000000000000000000000000000..0ca290b1678da0e719fc2776801d042cb33af522
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/test_mpnc_comparison_fv.input
@@ -0,0 +1,14 @@
+[TimeLoop]
+DtInitial = 250 # [s]
+TEnd = 1e4 # [s]
+
+[Grid]
+UpperRight = 60 40
+Cells = 24 16
+
+[LinearSolver]
+ResidualReduction = 1e-12
+
+[Problem]
+Name = obstacle
+
diff --git a/test/porousmediumflow/mpnc/implicit/CMakeLists.txt b/test/porousmediumflow/mpnc/implicit/CMakeLists.txt
index c07298e4446ae14932a961f31fa3164767051340..276f04953f999fffe9abba6aa79f53973d91847f 100644
--- a/test/porousmediumflow/mpnc/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/mpnc/implicit/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory(2p2ccomparison)
add_input_file_links()
dune_add_test(NAME test_mpnc_box
@@ -17,7 +18,6 @@ dune_add_test(NAME test_mpnc_tpfa
--files ${CMAKE_SOURCE_DIR}/test/references/obstaclecc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ccmpnc-00009.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_mpnc_tpfa test_mpnc.input -Problem.Name test_ccmpnc")
-
# build target for the full kinetic test problem
dune_add_test(NAME test_mpnckinetic_box
SOURCES test_boxmpnckinetic.cc