diff --git a/dumux/porenetwork/CMakeLists.txt b/dumux/porenetwork/CMakeLists.txt
index 6618be881349769efbb0e06c1075d1457dcb79b2..1b9434334cc1eaf70d856f0b82c1298ffd92dee7 100644
--- a/dumux/porenetwork/CMakeLists.txt
+++ b/dumux/porenetwork/CMakeLists.txt
@@ -4,6 +4,7 @@
add_subdirectory(1p)
add_subdirectory(1pnc)
add_subdirectory(2p)
+add_subdirectory(2pnc)
add_subdirectory(common)
add_subdirectory(solidenergy)
add_subdirectory(util)
diff --git a/test/porenetwork/2pnc/CMakeLists.txt b/test/porenetwork/2pnc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..04791d2e1a4f6f0d88bc61073c238e6c1f4a385e
--- /dev/null
+++ b/test/porenetwork/2pnc/CMakeLists.txt
@@ -0,0 +1,27 @@
+# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+add_input_file_links()
+dune_symlink_to_source_files(FILES grids)
+
+dumux_add_test(NAME test_pnm_2pnc
+ SOURCES main.cc
+ LABELS porenetwork
+ COMPILE_DEFINITIONS ISOTHERMAL=1
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_2pnc-reference.vtp
+ ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc-00010.vtp
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc")
+
+dumux_add_test(NAME test_pnm_2pnc_ni
+ SOURCES main.cc
+ LABELS porenetwork
+ COMPILE_DEFINITIONS ISOTHERMAL=0
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_2pnc_ni-reference.vtp
+ ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc_ni-00004.vtp
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc_ni params_ni.input -Problem.Name test_pnm_2pnc_ni")
diff --git a/test/porenetwork/2pnc/grids/1dGrid.dgf b/test/porenetwork/2pnc/grids/1dGrid.dgf
new file mode 100644
index 0000000000000000000000000000000000000000..e3aa06359ff65177b9b5b1a2f6e737d0d1a25ff0
--- /dev/null
+++ b/test/porenetwork/2pnc/grids/1dGrid.dgf
@@ -0,0 +1,20 @@
+DGF
+% SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+% SPDX-License-Identifier: CC-BY-4.0
+% Vertex parameters: PoreInscribedRadius PoreLabel
+% Element parameters: ThroatInscribedRadius ThroatLength ThroatLabel
+Vertex % Coordinates and volumes of the pore bodies
+parameters 2
+0 0 0 0.0002263 2.0
+1 0 0 0.0002263 -1.0
+2 0 0 0.0002263 -1.0
+3 0 0 0.0002263 -1.0
+4 0 0 0.0002263 3.0
+#
+SIMPLEX % Connections of the pore bodies (pore throats)
+parameters 3
+0 1 3.3304e-05 6.6609e-05 2
+1 2 3.3304e-05 6.6609e-05 -1
+2 3 3.3304e-05 6.6609e-05 -1
+3 4 3.3304e-05 6.6609e-05 3
+#
diff --git a/test/porenetwork/2pnc/main.cc b/test/porenetwork/2pnc/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b12fb47681e4c5c235706ff633fb9858ea8636f2
--- /dev/null
+++ b/test/porenetwork/2pnc/main.cc
@@ -0,0 +1,162 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ *
+ * \brief Test for the two-phase n-components pore-network 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 <dumux/assembly/fvassembler.hh>
+#include <dumux/common/initialize.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/io/grid/porenetwork/gridmanager.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/porenetwork/common/pnmvtkoutputmodule.hh>
+#include <dumux/porenetwork/2p/newtonsolver.hh>
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ using TypeTag = Properties::TTag::DrainageProblem;
+
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments and input file
+ ////////////////////////////////////////////////////////////
+
+ // parse command line arguments
+ Parameters::init(argc, argv);
+
+ //////////////////////////////////////////////////////////////////////
+ // try to create a grid (from the given grid file or the input file)
+ /////////////////////////////////////////////////////////////////////
+
+ using GridManager = Dumux::PoreNetwork::GridManager<3>;
+ GridManager gridManager;
+ gridManager.init();
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+ auto gridData = gridManager.getGridData();
+
+ // create the finite volume grid geometry
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ auto gridGeometry = std::make_shared<GridGeometry>(leafGridView, *gridData);
+
+ // the spatial parameters
+ using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
+ auto spatialParams = std::make_shared<SpatialParams>(gridGeometry);
+
+ // the problem (boundary conditions)
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry, spatialParams);
+
+ // the solution vector
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ SolutionVector x(gridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+ gridVariables->init(x);
+
+ // 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");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = getParam<Scalar>("TimeLoop.Restart", 0);
+
+ // initialize the vtk output module
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+ IOFields::initOutputModule(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, gridGeometry, gridVariables, timeLoop, xOld);
+
+ // the linear solver
+ using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = PoreNetwork::TwoPNewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // try solving the non-linear system
+ 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
+ if(problem->shouldWriteOutput(timeLoop->timeStepIndex(), *gridVariables))
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+ timeLoop->finalize(leafGridView.comm());
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
diff --git a/test/porenetwork/2pnc/params.input b/test/porenetwork/2pnc/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..42f13ccf19889cd6d57701686e9257347898b0fc
--- /dev/null
+++ b/test/porenetwork/2pnc/params.input
@@ -0,0 +1,31 @@
+[TimeLoop]
+DtInitial = 1e-5 # [s]
+TEnd = 0.01 # [s]
+
+[Grid]
+File = ./grids/1dGrid.dgf
+PoreGeometry = Cube
+ThroatCrossSectionShape = Square
+
+[Problem]
+Name = test_pnm_2pnc
+VtpOutputFrequency = 10 # Write every n-th time step. 0 only writes a file if an invasion / snap-off occurred. -1 writes every step
+EnableGravity = false
+CapillaryPressure = 5000
+UseFixedPressureAndSaturationBoundary = false
+Source = 1e-5
+
+[Vtk]
+AddVelocity = 1
+
+[Newton]
+MaxSteps = 10
+TargetSteps = 4
+MaxRelativeShift = 1e-5
+
+[InvasionState]
+Verbosity = true
+#AccuracyCriterion = 0.99
+
+[SpatialParams]
+HighSwRegularizationMethod = Spline
diff --git a/test/porenetwork/2pnc/params_ni.input b/test/porenetwork/2pnc/params_ni.input
new file mode 100644
index 0000000000000000000000000000000000000000..90c30695d6c7874b0b685f3078ca8932611cc0f8
--- /dev/null
+++ b/test/porenetwork/2pnc/params_ni.input
@@ -0,0 +1,31 @@
+[TimeLoop]
+DtInitial = 1e-5 # [s]
+TEnd = 0.01 # [s]
+
+[Grid]
+File = ./grids/1dGrid.dgf
+PoreGeometry = Cube
+ThroatCrossSectionShape = Square
+
+[Problem]
+Name = test_pnm2pnc_ni
+VtpOutputFrequency = 10 # Write every n-th time step. 0 only writes a file if an invasion / snap-off occurred
+CapillaryPressure = 5000
+EnableGravity = false
+UseFixedPressureAndSaturationBoundary = true
+InletTemperature = 288.15
+OutletTemperature = 283.15
+Source = 1e-5
+
+[Vtk]
+AddVelocity = 1
+
+[Newton]
+MaxSteps = 10
+TargetSteps = 4
+MaxRelativeShift = 1e-5
+
+[Component]
+SolidHeatCapacity = 1.0 # for compatibility
+SolidDensity = 1.0 # for compatibility
+SolidThermalConductivity = 1.0 # for compatibility
diff --git a/test/porenetwork/2pnc/problem.hh b/test/porenetwork/2pnc/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bf71e5a8ab85c43992e75e23db594c16a576c2f6
--- /dev/null
+++ b/test/porenetwork/2pnc/problem.hh
@@ -0,0 +1,224 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ *
+ * \brief A test problem for the two-phase n-components pore network model.
+ */
+#ifndef DUMUX_PNM_2P_NC_PROBLEM_HH
+#define DUMUX_PNM_2P_NC_PROBLEM_HH
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/components/air.hh>
+
+namespace Dumux {
+
+template <class TypeTag>
+class DrainageProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+
+ // copy some indices for convenience
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using Labels = GetPropType<TypeTag, Properties::Labels>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim<GridView::dimension>::Entity;
+
+public:
+ template<class SpatialParams>
+ DrainageProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<SpatialParams> spatialParams)
+ : ParentType(gridGeometry, spatialParams)
+ {
+ vtpOutputFrequency_ = getParam<int>("Problem.VtpOutputFrequency");
+ useFixedPressureAndSaturationBoundary_ = getParam<bool>("Problem.UseFixedPressureAndSaturationBoundary", false);
+ pc_ = getParam<Scalar>("Problem.CapillaryPressure");
+ source_ = getParam<Scalar>("Problem.Source");
+ inletPressure_ = getParam<Scalar>("Problem.InletPressure", 1e5);
+ outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1e5);
+#if !ISOTHERMAL
+ inletTemperature_ = getParam<Scalar>("Problem.InletTemperature", 288.15);
+ outletTemperature_ = getParam<Scalar>("Problem.OutletTemperature", 283.15);
+#endif
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ bool shouldWriteOutput(const int timeStepIndex, const GridVariables& gridVariables) const
+ {
+ if (vtpOutputFrequency_ < 0)
+ return true;
+
+ if (vtpOutputFrequency_ == 0)
+ return (timeStepIndex == 0 || gridVariables.gridFluxVarsCache().invasionState().hasChanged());
+ else
+ return (timeStepIndex % vtpOutputFrequency_ == 0 || gridVariables.gridFluxVarsCache().invasionState().hasChanged());
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+ //! Specifies which kind of boundary condition should be used for
+ //! which equation for a sub control volume on the boundary.
+ BoundaryTypes boundaryTypes(const Element& element, const SubControlVolume& scv) const
+ {
+ BoundaryTypes bcTypes;
+
+ // If a global phase pressure difference (pn,inlet - pw,outlet) with fixed saturations is specified, use a Dirichlet BC here
+ if (useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
+ bcTypes.setAllDirichlet();
+ else if (!useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
+ bcTypes.setAllNeumann();
+ else if (isOutletPore_(scv))
+ bcTypes.setAllDirichlet();
+
+ return bcTypes;
+ }
+
+
+ //! Evaluate the boundary conditions for a Dirichlet control volume.
+ PrimaryVariables dirichlet(const Element& element,
+ const SubControlVolume& scv) const
+ {
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = 1e5;
+ values[Indices::switchIdx] = 0.0;
+
+ // If a global phase pressure difference (pn,inlet - pw,outlet) is specified and the saturation shall also be fixed, apply:
+ // pw,inlet = pw,outlet = 1e5; pn,outlet = pw,outlet + pc(S=0) = pw,outlet; pn,inlet = pw,inlet + pc_
+ if (useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
+ {
+ values.setState(Indices::bothPhases);
+ values[Indices::pressureIdx] = inletPressure_;
+ values[Indices::switchIdx] = 1.0 - this->spatialParams().fluidMatrixInteraction(element, scv, int()/*dummyElemsol*/).sw(pc_);
+#if !ISOTHERMAL
+ values[Indices::temperatureIdx] = inletTemperature_;
+#endif
+ }
+ else if (isOutletPore_(scv))
+ {
+ values.setState(Indices::firstPhaseOnly);
+ values[Indices::pressureIdx] = outletPressure_;
+ values[Indices::switchIdx] = 0.0;
+#if !ISOTHERMAL
+ values[Indices::temperatureIdx] = outletTemperature_;
+#endif
+ }
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ //! Evaluate the source term for all phases within a given sub-control-volume.
+ PrimaryVariables source(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume& scv) const
+ {
+ PrimaryVariables values(0.0);
+
+ // for isothermal case, we fix injection rate of non-wetting phase at inlet
+ // for non-isothermal case, we fix injection of air enthalpy at inlet
+ if (!useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
+ {
+ values[Indices::conti0EqIdx + 1] = source_/scv.volume();
+#if !ISOTHERMAL
+ const auto pressure = elemVolVars[scv].pressure(1);
+ const auto airEnthalpy = Components::Air<Scalar>::gasEnthalpy(inletTemperature_, pressure);
+ values[Indices::temperatureIdx] = airEnthalpy * source_ * Components::Air<Scalar>::molarMass()/scv.volume();
+#endif
+ }
+
+ return values;
+ }
+ // \}
+
+ //! Evaluate the initial value for a control volume.
+ PrimaryVariables initial(const Vertex& vertex) const
+ {
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = outletPressure_;
+
+ // get global index of pore
+ const auto dofIdxGlobal = this->gridGeometry().vertexMapper().index(vertex);
+ if (isInletPore_(dofIdxGlobal))
+ {
+ values.setState(Indices::firstPhaseOnly);
+ values[Indices::switchIdx] = 0.0;
+ }
+ else
+ {
+ values.setState(Indices::firstPhaseOnly);
+ values[Indices::switchIdx] = 0.0;
+ }
+
+#if !ISOTHERMAL
+ values[Indices::temperatureIdx] = outletTemperature_;
+#endif
+
+ return values;
+ }
+
+ //! Evaluate the initial invasion state of a pore throat
+ bool initialInvasionState(const Element& element) const
+ { return false; }
+
+ // \}
+
+private:
+
+ bool isInletPore_(const SubControlVolume& scv) const
+ {
+ return isInletPore_(scv.dofIndex());
+ }
+
+ bool isInletPore_(const std::size_t dofIdxGlobal) const
+ {
+ return this->gridGeometry().poreLabel(dofIdxGlobal) == Labels::inlet;
+ }
+
+ bool isOutletPore_(const SubControlVolume& scv) const
+ {
+ return this->gridGeometry().poreLabel(scv.dofIndex()) == Labels::outlet;
+ }
+
+ int vtpOutputFrequency_;
+ bool useFixedPressureAndSaturationBoundary_;
+ Scalar pc_;
+ Scalar source_;
+ Scalar inletPressure_;
+ Scalar outletPressure_;
+#if !ISOTHERMAL
+ Scalar inletTemperature_;
+ Scalar outletTemperature_;
+#endif
+};
+} //end namespace Dumux
+
+#endif
diff --git a/test/porenetwork/2pnc/properties.hh b/test/porenetwork/2pnc/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..19c02ecff02a5af0e6be0634600e43224a9b2946
--- /dev/null
+++ b/test/porenetwork/2pnc/properties.hh
@@ -0,0 +1,77 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ *
+ * \brief The properties for the two-phase n-components pore network model.
+ */
+#ifndef DUMUX_PNM_2P_NC_PROPERTIES_HH
+#define DUMUX_PNM_2P_NC_PROPERTIES_HH
+
+#include <dune/foamgrid/foamgrid.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/porenetwork/2pnc/model.hh>
+#include <dumux/porenetwork/2p/spatialparams.hh>
+#include <dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh>
+
+#include <dumux/common/properties.hh>
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/porenetwork/common/utilities.hh>
+
+#include "problem.hh"
+#include "spatialparams.hh"
+
+//////////
+// Specify the properties
+//////////
+namespace Dumux::Properties {
+
+// Create new type tags
+namespace TTag {
+#if ISOTHERMAL
+struct DrainageProblem { using InheritsFrom = std::tuple<PNMTwoPNC>; };
+#else
+struct DrainageProblem { using InheritsFrom = std::tuple<PNMTwoPNCNI>; };
+#endif
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DrainageProblem> { using type = DrainageProblem<TypeTag>; };
+
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DrainageProblem>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using H2OType = Components::TabulatedComponent<Components::H2O<Scalar> >;
+ using Policy = FluidSystems::H2OAirDefaultPolicy<false/*fastButSimplifiedRelations*/>;
+ using type = FluidSystems::H2OAir<Scalar, H2OType, Policy, true/*useKelvinVaporPressure*/>;
+
+};
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DrainageProblem>
+{
+private:
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using LocalRules = PoreNetwork::FluidMatrix::MultiShapeTwoPLocalRules<Scalar>;
+public:
+ using type = PoreNetwork::TwoPNCDrainageSpatialParams<GridGeometry, Scalar, LocalRules>;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DrainageProblem> { using type = Dune::FoamGrid<1, 3>; };
+
+} //end namespace Dumux::Properties
+
+#endif
diff --git a/test/porenetwork/2pnc/spatialparams.hh b/test/porenetwork/2pnc/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..86dc93f5b4c9a1b18ebe1d41e6d035e7de5a053d
--- /dev/null
+++ b/test/porenetwork/2pnc/spatialparams.hh
@@ -0,0 +1,69 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup PoreNetworkModels
+ * \ingroup SpatialParameters
+ * \brief Spatial parameters for a 2pnc pore-network model
+ */
+
+#ifndef DUMUX_PNM_2P_NC_DRAINAGE_SPATIAL_PARAMS_HH
+#define DUMUX_PNM_2P_NC_DRAINAGE_SPATIAL_PARAMS_HH
+
+#include <dumux/porenetwork/2p/spatialparams.hh>
+
+namespace Dumux::PoreNetwork {
+
+template<class GridGeometry, class Scalar, class MaterialLawT>
+class TwoPNCDrainageSpatialParams
+: public TwoPSpatialParams<GridGeometry, Scalar, MaterialLawT, TwoPNCDrainageSpatialParams<GridGeometry, Scalar, MaterialLawT>>
+{
+ using ParentType = TwoPSpatialParams<GridGeometry, Scalar, MaterialLawT, TwoPNCDrainageSpatialParams<GridGeometry, Scalar, MaterialLawT>>;
+
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+public:
+ using PermeabilityType = Scalar;
+ using ParentType::ParentType;
+
+ TwoPNCDrainageSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ temperature_ = getParam<Scalar>("SpatialParams.Temperature", 283.15);
+ gamma_ = getParam<Scalar>("SpatialParams.SurfaceTension", 0.0725); // default to surface tension of water/air
+ theta_ = getParam<Scalar>("SpatialParams.ContactAngle", 0.0);
+ }
+
+ //! \brief Function for defining the temperature
+ Scalar temperatureAtPos(const GlobalPosition& globalPos) const
+ { return temperature_; }
+
+ //! \brief Function for defining the wetting phase
+ template<class FluidSystem>
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ { return FluidSystem::phase0Idx; }
+
+
+ //! \brief Function for defining the contact angle
+ int contactAngleAtPos(const GlobalPosition& globalPos) const
+ { return theta_; }
+
+ //! \brief Function for defining the surface tension
+ Scalar surfaceTensionAtPos(const GlobalPosition& globalPos) const
+ { return gamma_; }
+
+private:
+ Scalar temperature_;
+ Scalar gamma_;
+ Scalar theta_;
+};
+} // end namespace Dumux::PoreNetwork
+
+#endif
diff --git a/test/porenetwork/CMakeLists.txt b/test/porenetwork/CMakeLists.txt
index 4cb5469ebdb3bb3d4f3ba1ca00c9383c560f2ca5..73edd9468deca4edf5a6fe4d4f66e18ffff3c4fb 100644
--- a/test/porenetwork/CMakeLists.txt
+++ b/test/porenetwork/CMakeLists.txt
@@ -5,3 +5,4 @@ add_subdirectory(1p)
add_subdirectory(1pnc)
add_subdirectory(2p)
add_subdirectory(solidenergy)
+add_subdirectory(2pnc)
diff --git a/test/references/test_pnm_2pnc-reference.vtp b/test/references/test_pnm_2pnc-reference.vtp
new file mode 100644
index 0000000000000000000000000000000000000000..5ad168b997c1c8647e8139f31177ac71791ac7b2
--- /dev/null
+++ b/test/references/test_pnm_2pnc-reference.vtp
@@ -0,0 +1,130 @@
+<?xml version="1.0"?>
+<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">
+ <PolyData>
+ <Piece NumberOfLines="4" NumberOfPoints="5">
+ <PointData Scalars="S_liq">
+ <DataArray type="Float32" Name="S_liq" NumberOfComponents="1" format="ascii">
+ 0.0374111 0.0474591 0.0343665 0.0391619 1
+ </DataArray>
+ <DataArray type="Float32" Name="p_liq" NumberOfComponents="1" format="ascii">
+ 99081.6 99126.6 97908.1 97768.2 100000
+ </DataArray>
+ <DataArray type="Float32" Name="rho_liq" NumberOfComponents="1" format="ascii">
+ 999.712 999.712 999.711 999.711 999.701
+ </DataArray>
+ <DataArray type="Float32" Name="mob_liq" NumberOfComponents="1" format="ascii">
+ 765.754 765.754 765.753 765.753 765.754
+ </DataArray>
+ <DataArray type="Float32" Name="S_gas" NumberOfComponents="1" format="ascii">
+ 0.962589 0.952541 0.965634 0.960838 0
+ </DataArray>
+ <DataArray type="Float32" Name="p_gas" NumberOfComponents="1" format="ascii">
+ 101923 101441 100971 100498 100000
+ </DataArray>
+ <DataArray type="Float32" Name="rho_gas" NumberOfComponents="1" format="ascii">
+ 1.24808 1.24214 1.23636 1.23055 0.00940657
+ </DataArray>
+ <DataArray type="Float32" Name="mob_gas" NumberOfComponents="1" format="ascii">
+ 56885.7 56887.5 56889.3 56891.1 103741
+ </DataArray>
+ <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+ 2841.63 2314.36 3062.68 2730.15 0
+ </DataArray>
+ <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="x^H2O_liq" NumberOfComponents="1" format="ascii">
+ 0.999981 0.999982 0.999982 0.999982 1
+ </DataArray>
+ <DataArray type="Float32" Name="X^H2O_liq" NumberOfComponents="1" format="ascii">
+ 0.99997 0.99997 0.99997 0.999971 1
+ </DataArray>
+ <DataArray type="Float32" Name="x^Air_liq" NumberOfComponents="1" format="ascii">
+ 1.8543e-05 1.84542e-05 1.83676e-05 1.82806e-05 0
+ </DataArray>
+ <DataArray type="Float32" Name="X^Air_liq" NumberOfComponents="1" format="ascii">
+ 2.98082e-05 2.96654e-05 2.95262e-05 2.93864e-05 0
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar_liq" NumberOfComponents="1" format="ascii">
+ 55492.1 55492.1 55492.1 55492.1 55492.2
+ </DataArray>
+ <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
+ 0.0120499 0.0121072 0.0121635 0.0122207 0.0122818
+ </DataArray>
+ <DataArray type="Float32" Name="X^H2O_gas" NumberOfComponents="1" format="ascii">
+ 0.00753016 0.00756613 0.00760153 0.00763743 0.0122818
+ </DataArray>
+ <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
+ 0.98795 0.987893 0.987836 0.987779 0
+ </DataArray>
+ <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
+ 0.99247 0.992434 0.992399 0.992363 0
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar_gas" NumberOfComponents="1" format="ascii">
+ 43.2939 43.089 42.8893 42.6886 0.522147
+ </DataArray>
+ <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
+ 3 3 3 3 1
+ </DataArray>
+ <DataArray type="Float32" Name="poreInscribedRadius" NumberOfComponents="1" format="ascii">
+ 0.0002263 0.0002263 0.0002263 0.0002263 0.0002263
+ </DataArray>
+ <DataArray type="Float32" Name="coordinationNumber" NumberOfComponents="1" format="ascii">
+ 1 2 2 2 1
+ </DataArray>
+ <DataArray type="Float32" Name="poreLabel" NumberOfComponents="1" format="ascii">
+ 2 -1 -1 -1 3
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank" Vectors="velocity_liq (m/s)">
+ <DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
+ -0.00296711 -0 -0 0.0692308 0 0 0.00794862 0 0 -0.159574 -0 -0
+ </DataArray>
+ <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
+ 60.2882 0 0 59.377 0 0 59.6628 0 0 61.9152 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="throatLabel" NumberOfComponents="1" format="ascii">
+ 2 -1 -1 3
+ </DataArray>
+ <DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
+ 3.3304e-05 3.3304e-05 3.3304e-05 3.3304e-05
+ </DataArray>
+ <DataArray type="Float32" Name="throatLength" NumberOfComponents="1" format="ascii">
+ 6.6609e-05 6.6609e-05 6.6609e-05 6.6609e-05
+ </DataArray>
+ <DataArray type="Float32" Name="pcEntry" NumberOfComponents="1" format="ascii">
+ 4106.16 4106.16 4106.16 4106.16
+ </DataArray>
+ <DataArray type="Float32" Name="transmissibilityW" NumberOfComponents="1" format="ascii">
+ 4.81607e-17 3.56908e-17 3.56908e-17 5.65214e-17
+ </DataArray>
+ <DataArray type="Float32" Name="transmissibilityN" NumberOfComponents="1" format="ascii">
+ 8.5211e-15 8.78119e-15 8.78119e-15 8.36682e-15
+ </DataArray>
+ <DataArray type="Float32" Name="volumeFluxW" NumberOfComponents="1" format="ascii">
+ 1.65794e-12 3.33017e-11 3.82348e-12 9.65956e-11
+ </DataArray>
+ <DataArray type="Float32" Name="volumeFluxN" NumberOfComponents="1" format="ascii">
+ 2.33789e-07 2.34872e-07 2.36002e-07 2.37215e-07
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 1 0 0 2 0 0 3 0 0
+ 4 0 0
+ </DataArray>
+ </Points>
+ <Lines>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 1 2 2 3 3 4
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 2 4 6 8
+ </DataArray>
+ </Lines>
+ </Piece>
+ </PolyData>
+</VTKFile>
diff --git a/test/references/test_pnm_2pnc_ni-reference.vtp b/test/references/test_pnm_2pnc_ni-reference.vtp
new file mode 100644
index 0000000000000000000000000000000000000000..be60712cea5a4e52236baa9770af90938e5835ed
--- /dev/null
+++ b/test/references/test_pnm_2pnc_ni-reference.vtp
@@ -0,0 +1,133 @@
+<?xml version="1.0"?>
+<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">
+ <PolyData>
+ <Piece NumberOfLines="4" NumberOfPoints="5">
+ <PointData Scalars="S_liq">
+ <DataArray type="Float32" Name="S_liq" NumberOfComponents="1" format="ascii">
+ 0.0200785 0.0261951 0.948722 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="p_liq" NumberOfComponents="1" format="ascii">
+ 100000 101052 104283 102142 100000
+ </DataArray>
+ <DataArray type="Float32" Name="rho_liq" NumberOfComponents="1" format="ascii">
+ 999.111 999.711 999.715 999.705 999.701
+ </DataArray>
+ <DataArray type="Float32" Name="mob_liq" NumberOfComponents="1" format="ascii">
+ 879.067 766.229 765.823 765.781 765.754
+ </DataArray>
+ <DataArray type="Float32" Name="S_gas" NumberOfComponents="1" format="ascii">
+ 0.979922 0.973805 0.0512775 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="p_gas" NumberOfComponents="1" format="ascii">
+ 105000 104963 104925 102142 100000
+ </DataArray>
+ <DataArray type="Float32" Name="rho_gas" NumberOfComponents="1" format="ascii">
+ 1.26143 1.28536 1.28499 0.379949 0.00940657
+ </DataArray>
+ <DataArray type="Float32" Name="mob_gas" NumberOfComponents="1" format="ascii">
+ 56207.5 56871.4 56874.1 57633.3 103741
+ </DataArray>
+ <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+ 5000 3911.23 641.727 0 -0
+ </DataArray>
+ <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="x^H2O_liq" NumberOfComponents="1" format="ascii">
+ 0.999983 0.999981 0.999981 0.999994 1
+ </DataArray>
+ <DataArray type="Float32" Name="X^H2O_liq" NumberOfComponents="1" format="ascii">
+ 0.999973 0.999969 0.999969 0.999991 1
+ </DataArray>
+ <DataArray type="Float32" Name="x^Air_liq" NumberOfComponents="1" format="ascii">
+ 1.70916e-05 1.90931e-05 1.90944e-05 5.54692e-06 0
+ </DataArray>
+ <DataArray type="Float32" Name="X^Air_liq" NumberOfComponents="1" format="ascii">
+ 2.74751e-05 3.06925e-05 3.06946e-05 8.91683e-06 0
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar_liq" NumberOfComponents="1" format="ascii">
+ 55458.8 55492.1 55492.3 55492.2 55492.2
+ </DataArray>
+ <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
+ 0.0162449 0.0117178 0.0117075 0.0120252 0.0122818
+ </DataArray>
+ <DataArray type="Float32" Name="X^H2O_gas" NumberOfComponents="1" format="ascii">
+ 0.0101679 0.00732173 0.00731526 0.00759294 0.0122818
+ </DataArray>
+ <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
+ 0.983755 0.988282 0.988293 0.294909 0
+ </DataArray>
+ <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
+ 0.989832 0.992678 0.992685 0.299341 0
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar_gas" NumberOfComponents="1" format="ascii">
+ 43.8272 44.5816 44.5685 13.3171 0.522147
+ </DataArray>
+ <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
+ 3 3 3 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="poreInscribedRadius" NumberOfComponents="1" format="ascii">
+ 0.0002263 0.0002263 0.0002263 0.0002263 0.0002263
+ </DataArray>
+ <DataArray type="Float32" Name="T" NumberOfComponents="1" format="ascii">
+ 288.15 283.172 283.153 283.151 283.15
+ </DataArray>
+ <DataArray type="Float32" Name="coordinationNumber" NumberOfComponents="1" format="ascii">
+ 1 2 2 2 1
+ </DataArray>
+ <DataArray type="Float32" Name="poreLabel" NumberOfComponents="1" format="ascii">
+ 2 -1 -1 -1 3
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank" Vectors="velocity_liq (m/s)">
+ <DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
+ -0.0224304 -0 -0 -0.112597 -0 -0 3.83914 0 0 3.83908 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
+ 4.82473 0 0 4.88488 0 0 0 0 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="throatLabel" NumberOfComponents="1" format="ascii">
+ 2 -1 -1 3
+ </DataArray>
+ <DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
+ 3.3304e-05 3.3304e-05 3.3304e-05 3.3304e-05
+ </DataArray>
+ <DataArray type="Float32" Name="throatLength" NumberOfComponents="1" format="ascii">
+ 6.6609e-05 6.6609e-05 6.6609e-05 6.6609e-05
+ </DataArray>
+ <DataArray type="Float32" Name="pcEntry" NumberOfComponents="1" format="ascii">
+ 4106.16 4106.16 4106.16 4106.16
+ </DataArray>
+ <DataArray type="Float32" Name="transmissibilityW" NumberOfComponents="1" format="ascii">
+ 5.02435e-18 1.34185e-17 1.03853e-14 1.03853e-14
+ </DataArray>
+ <DataArray type="Float32" Name="transmissibilityN" NumberOfComponents="1" format="ascii">
+ 9.81542e-15 9.41615e-15 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="volumeFluxW" NumberOfComponents="1" format="ascii">
+ 4.04824e-12 3.32099e-11 1.70328e-08 1.70325e-08
+ </DataArray>
+ <DataArray type="Float32" Name="volumeFluxN" NumberOfComponents="1" format="ascii">
+ 2.05347e-08 2.02316e-08 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 1 0 0 2 0 0 3 0 0
+ 4 0 0
+ </DataArray>
+ </Points>
+ <Lines>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 1 2 2 3 3 4
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 2 4 6 8
+ </DataArray>
+ </Lines>
+ </Piece>
+ </PolyData>
+</VTKFile>