diff --git a/dumux/porenetwork/2pnc/CMakeLists.txt b/dumux/porenetwork/2pnc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3de09eeabc5386d3152c4e3b9a581734571246ad
--- /dev/null
+++ b/dumux/porenetwork/2pnc/CMakeLists.txt
@@ -0,0 +1,6 @@
+# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+file(GLOB DUMUX_PORENETWORK_2P_NC_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PORENETWORK_2P_NC_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porenetwork/2pnc)
diff --git a/dumux/porenetwork/2pnc/iofields.hh b/dumux/porenetwork/2pnc/iofields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f14cb1f9a15f02fb38d168ef290623aca983c480
--- /dev/null
+++ b/dumux/porenetwork/2pnc/iofields.hh
@@ -0,0 +1,67 @@
+// -*- 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 PNMTwoPNCModel
+ * \copydoc Dumux::PoreNetwork::TwoPNCIOFields
+ */
+#ifndef DUMUX_PNM_2P_NC_IO_FIELDS_HH
+#define DUMUX_PNM_2P_NC_IO_FIELDS_HH
+
+#include <dumux/porousmediumflow/2pnc/iofields.hh>
+#include <dumux/porenetwork/common/iofields.hh>
+
+namespace Dumux::PoreNetwork {
+
+/*!
+ * \ingroup PNMTwoPNCModel
+ * \brief Adds output fields specific to the PNM 2pnc model
+ */
+class TwoPNCIOFields
+{
+public:
+ template <class OutputModule>
+ static void initOutputModule(OutputModule& out)
+ {
+ // use default fields from the 2pnc model
+ Dumux::TwoPNCIOFields::initOutputModule(out);
+
+ using VolumeVariables = typename OutputModule::VolumeVariables;
+ using FS = typename VolumeVariables::FluidSystem;
+
+ CommonIOFields::initOutputModule(out);
+
+ out.addFluxVariable([](const auto& fluxVars, const auto& fluxVarsCache)
+ { return fluxVarsCache.pcEntry(); }, "pcEntry");
+
+ out.addFluxVariable([](const auto& fluxVars, const auto& fluxVarsCache)
+ { return fluxVarsCache.transmissibility(FS::phase0Idx); }, "transmissibilityW");
+
+ out.addFluxVariable([](const auto& fluxVars, const auto& fluxVarsCache)
+ { return fluxVarsCache.transmissibility(FS::phase1Idx); }, "transmissibilityN");
+
+ auto volumeFluxW = [](const auto& fluxVars, const auto& fluxVarsCache)
+ {
+ auto upwindTerm = [](const auto& volVars) { return volVars.mobility(FS::phase0Idx); };
+ using std::abs;
+ return abs(fluxVars.advectiveFlux(FS::phase0Idx, upwindTerm));
+ };
+ out.addFluxVariable(volumeFluxW, "volumeFluxW");
+
+ auto volumeFluxN = [](const auto& fluxVars, const auto& fluxVarsCache)
+ {
+ auto upwindTerm = [](const auto& volVars) { return volVars.mobility(FS::phase1Idx); };
+ using std::abs;
+ return abs(fluxVars.advectiveFlux(FS::phase1Idx, upwindTerm));
+ };
+ out.addFluxVariable(volumeFluxN, "volumeFluxN");
+ }
+};
+
+} // end namespace Dumux::PoreNetwork
+
+#endif
diff --git a/dumux/porenetwork/2pnc/model.hh b/dumux/porenetwork/2pnc/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5a73becb375e6b865995eccbcc7665c1a4fbe6c1
--- /dev/null
+++ b/dumux/porenetwork/2pnc/model.hh
@@ -0,0 +1,259 @@
+// -*- 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 PNMTwoPNCModel
+ * \brief A two-phase multi-compositional pore-network model using fully implicit scheme.
+ *
+ * A mass balance equation is formulated for each pore body \f$i\f$ and each component \f$\kappa\f$:
+ *
+ * \f[
+ * V_i \frac{\partial (\sum_{\alpha} (x^\kappa_\alpha \varrho_\alpha S_\alpha)_i) }{\partial t} + \sum_\alpha \sum_j (x^\kappa_\alpha \varrho_\alpha Q_\alpha)_{ij}
+ * + \sum_\alpha \sum_j (A_\alpha j^\kappa_\alpha)_{ij}= (V q^\kappa)_i ~.
+ * \f]
+ *
+ * \f$V_i\f$ is the pore body volume, \f$x^\kappa_\alpha\f$ resepresents the mole fraction of component \f$\kappa\f$ in phase \f$\alpha\f$ , the advective mass flow \f$(\varrho_\alpha Q_\alpha)_{ij}\f$
+ * through throat \f$ij\f$ can be based on the fluid phase density \f$\varrho\f$ either of the upstream pore body \f$i\f$ or \f$j\f$ (upwinding) or on the respective averaged value.
+ * \f$j^\kappa_\alpha\f$ represents the diffuive flux vector and \f$A_\alpha\f$ the area of throat cross section occupied by phase \f$\alpha\f$. \f$q_\alpha\f$ is a mass sink or source term defined on pore body \f$i\f$.
+ *
+ * Per default, the volume flow rate \f$Q_{\alpha,ij}\f$ follows a linear Hagen-Poiseuille-type law (PoreNetworkModel::CreepingFlow) which is only valid for \f$Re < 1\f$:
+ *
+ * \f[
+ * Q_{\alpha,ij} = g_{\alpha, ij} (p_{\alpha, i} - p_{\alpha, j} + \Psi_\alpha) ~.
+ * \f]
+ *
+ * \f$g_{\alpha,ij}\f$ is a suitable throat conductance value that takes into account the presence/saturation of the individual phases while \f$p_{\alpha,i}\f$ and \f$p_{\alpha,j}\f$ are averaged pore body phase pressures.
+ *
+ * The (optional) influence of gravity is given by
+ *
+ * \f[
+ * \Psi_\alpha = \varrho_\alpha \mathbf{g} \cdot (\mathbf{x}_i - \mathbf{x}_j),
+ * \f]
+ *
+ * where \f$\mathbf{x}_i - \mathbf{x}_j\f$ is the distance vector between the centers of pore bodies \f$i\f$ and \f$j\f$ and \f$\mathbf{g}\f$ is the gravity vector.
+ *
+ * The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
+ * problem file. Make sure that the according units are used in the problem setup. useMoles is set to true by default.
+ *
+ * The primary variables are the pressure \f$p_\alpha\f$ for phase \f$\alpha\f$ and the mole or mass fraction of dissolved components in phase \f$\alpha\f$ \f$x^\kappa_\alpha\f$ or \f$X^\kappa_\alpha\f$.
+ */
+
+#ifndef DUMUX_PNM_2P_NC_MODEL_HH
+#define DUMUX_PNM_2P_NC_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/porenetwork/properties.hh>
+#include <dumux/flux/porenetwork/fickslaw.hh>
+
+#include <dumux/porousmediumflow/nonisothermal/model.hh>
+#include <dumux/porousmediumflow/nonisothermal/indices.hh>
+#include <dumux/porousmediumflow/nonisothermal/iofields.hh>
+
+#include <dumux/flux/porenetwork/advection.hh>
+#include <dumux/porousmediumflow/2pnc/model.hh>
+#include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
+#include <dumux/material/fluidmatrixinteractions/porenetwork/throat/transmissibility1p.hh>
+#include <dumux/material/fluidmatrixinteractions/porenetwork/throat/transmissibility2p.hh>
+#include <dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh>
+
+#include <dumux/porenetwork/2p/fluxvariablescache.hh>
+#include <dumux/porenetwork/2p/gridfluxvariablescache.hh>
+#include <dumux/porenetwork/2p/spatialparams.hh>
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+
+
+#include "iofields.hh"
+#include "volumevariables.hh"
+
+namespace Dumux::Properties{
+// \{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit two-phase n-components problems
+// Create new type tags
+namespace TTag {
+struct PNMTwoPNC { using InheritsFrom = std::tuple<PoreNetworkModel, TwoPNC>; };
+
+//! The type tags for the corresponding non-isothermal problems
+struct PNMTwoPNCNI { using InheritsFrom = std::tuple<PNMTwoPNC>; };
+} // end namespace TTag
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+//! Set the volume variables property
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::PNMTwoPNC>
+{
+private:
+ using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
+ using FST = GetPropType<TypeTag, Properties::FluidState>;
+ using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
+ using SST = GetPropType<TypeTag, Properties::SolidState>;
+ using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
+ using MT = GetPropType<TypeTag, Properties::ModelTraits>;
+ using DM = typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod;
+ static constexpr bool enableIS = getPropValue<TypeTag, Properties::EnableBoxInterfaceSolver>();
+ // class used for scv-wise reconstruction of nonwetting phase saturations
+ using SR = TwoPScvSaturationReconstruction<DM, enableIS>;
+ using BaseTraits = TwoPVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, SR>;
+
+ using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+ using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+
+ template<class BaseTraits, class DT, class EDM>
+ struct NCTraits : public BaseTraits
+ {
+ using DiffusionType = DT;
+ using EffectiveDiffusivityModel = EDM;
+ };
+public:
+ using type = PoreNetwork::TwoPNCVolumeVariables<NCTraits<BaseTraits, DT, EDM>>;
+};
+
+//! The primary variables vector for the 2pnc model
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::PNMTwoPNC>
+{
+private:
+ using PrimaryVariablesVector = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
+ GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+public:
+ using type = SwitchablePrimaryVariables<PrimaryVariablesVector, int>;
+};
+
+//! The spatial parameters to be employed.
+//! Use PNMTwoPSpatialParams by default.
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::PNMTwoPNC>
+{
+private:
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using LocalRules = PoreNetwork::FluidMatrix::MultiShapeTwoPLocalRules<Scalar>;
+public:
+ using type = PoreNetwork::TwoPDefaultSpatialParams<GridGeometry, Scalar, LocalRules>;
+};
+
+//! The advection type
+template<class TypeTag>
+struct AdvectionType<TypeTag, TTag::PNMTwoPNC>
+{
+private:
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using S = PoreNetwork::TransmissibilityPatzekSilin<Scalar>;
+ using W = PoreNetwork::WettingLayerTransmissibility::RansohoffRadke<Scalar>;
+ using N = PoreNetwork::NonWettingPhaseTransmissibility::BakkeOren<Scalar>;
+
+public:
+ using type = PoreNetwork::CreepingFlow<Scalar, S, W, N>;
+};
+
+template<class TypeTag>
+struct IOFields<TypeTag, TTag::PNMTwoPNC> { using type = PoreNetwork::TwoPNCIOFields; };
+
+//! The grid flux variables cache vector class
+template<class TypeTag>
+struct GridFluxVariablesCache<TypeTag, TTag::PNMTwoPNC>
+{
+private:
+ static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluxVariablesCache = GetPropTypeOr<TypeTag,
+ Properties::FluxVariablesCache, FluxVariablesCaching::EmptyCache<Scalar>
+ >;
+ using Traits = PoreNetwork::PNMTwoPDefaultGridFVCTraits<Problem, FluxVariablesCache>;
+public:
+ using type = PoreNetwork::PNMTwoPGridFluxVariablesCache<Problem, FluxVariablesCache, enableCache, Traits>;
+};
+
+//! The flux variables cache
+template<class TypeTag>
+struct FluxVariablesCache<TypeTag, TTag::PNMTwoPNC>
+{ using type = PoreNetwork::TwoPFluxVariablesCache<GetPropType<TypeTag, Properties::AdvectionType>>; };
+
+//! We use fick's law as the default for the diffusive fluxes
+template<class TypeTag>
+struct MolecularDiffusionType<TypeTag, TTag::PNMTwoPNC>
+{
+private:
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+public:
+ using type = PoreNetwork::PNMFicksLaw<Scalar, ModelTraits::numFluidPhases(), ModelTraits::numFluidComponents()>;
+};
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::PNMTwoPNCNI>
+{
+private:
+ //! we use the number of components specified by the fluid system here
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ static_assert(FluidSystem::numPhases == 2, "Only fluid systems with 2 fluid phases are supported by the 2p-nc model!");
+ using IsothermalTraits = TwoPNCModelTraits<FluidSystem::numComponents,
+ getPropValue<TypeTag, Properties::UseMoles>(),
+ getPropValue<TypeTag, Properties::SetMoleFractionsForFirstPhase>(),
+ getPropValue<TypeTag, Properties::Formulation>(), getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
+public:
+ using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
+}; //!< The model traits of the non-isothermal model
+
+//! Set the volume variables property
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::PNMTwoPNCNI>
+{
+private:
+ using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
+ using FST = GetPropType<TypeTag, Properties::FluidState>;
+ using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
+ using SST = GetPropType<TypeTag, Properties::SolidState>;
+ using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
+ using MT = GetPropType<TypeTag, Properties::ModelTraits>;
+ using DM = typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod;
+ static constexpr bool enableIS = getPropValue<TypeTag, Properties::EnableBoxInterfaceSolver>();
+ // class used for scv-wise reconstruction of nonwetting phase saturations
+ using SR = TwoPScvSaturationReconstruction<DM, enableIS>;
+ using BaseTraits = TwoPVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, SR>;
+
+ using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+ using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+ using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
+
+ template<class BaseTraits, class DT, class EDM, class ETCM>
+ struct NCNITraits : public BaseTraits
+ {
+ using DiffusionType = DT;
+ using EffectiveDiffusivityModel = EDM;
+ using EffectiveThermalConductivityModel = ETCM;
+ };
+public:
+ using type = PoreNetwork::TwoPNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
+};
+
+//! Set non-isothermal output fields
+template<class TypeTag>
+struct IOFields<TypeTag, TTag::PNMTwoPNCNI> { using type = EnergyIOFields<PoreNetwork::TwoPNCIOFields>; };
+
+// Somerton is used as default model to compute the effective thermal heat conductivity
+template<class TypeTag>
+struct ThermalConductivityModel<TypeTag, TTag::PNMTwoPNCNI>
+{
+ using type = ThermalConductivitySomertonTwoP<GetPropType<TypeTag, Properties::Scalar>>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/dumux/porenetwork/2pnc/volumevariables.hh b/dumux/porenetwork/2pnc/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b6fd98bebd498e502a260cf3207a016b6b3da0c0
--- /dev/null
+++ b/dumux/porenetwork/2pnc/volumevariables.hh
@@ -0,0 +1,89 @@
+// -*- 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 PNMTwoPNCModel
+ * \copydoc Dumux::PoreNetwork::TwoPNCVolumeVariables
+ */
+
+#ifndef DUMUX_PNM_2P_NC_VOLUME_VARIABLES_HH
+#define DUMUX_PNM_2P_NC_VOLUME_VARIABLES_HH
+
+#include <dumux/porousmediumflow/2pnc/volumevariables.hh>
+
+namespace Dumux::PoreNetwork {
+
+/*!
+ * \ingroup PNMTwoPNCModel
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the two-phase n-components model.
+ */
+template <class Traits>
+class TwoPNCVolumeVariables
+: public Dumux::TwoPNCVolumeVariables<Traits>
+{
+ using ParentType = Dumux::TwoPNCVolumeVariables<Traits>;
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+
+public:
+ //! Export type of fluid system
+ using FluidSystem = typename Traits::FluidSystem;
+ //! Export type of fluid state
+ using FluidState = typename Traits::FluidState;
+ //! Export type of solid state
+ using SolidState = typename Traits::SolidState;
+ //! Export type of solid system
+ using SolidSystem = typename Traits::SolidSystem;
+
+ /*!
+ * \brief Updates all quantities for a given control volume.
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to
+ * be simulated
+ * \param element An element which contains part of the control volume
+ * \param scv The sub control volume
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ void update(const ElemSol &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const Scv& scv)
+ {
+ ParentType::update(elemSol, problem, element, scv);
+ poreInscribedRadius_ = problem.spatialParams().poreInscribedRadius(element, scv, elemSol);
+ poreVolume_ = problem.gridGeometry().poreVolume(scv.dofIndex()) * this->porosity();
+ surfaceTension_ = problem.spatialParams().surfaceTension(element, scv, elemSol);
+ }
+
+ /*!
+ * \brief Returns the pore's inscribed radius.
+ */
+ Scalar poreInscribedRadius() const
+ { return poreInscribedRadius_; }
+
+ /*!
+ * \brief Returns the pore volume. // TODO should this be a fraction only?
+ */
+ Scalar poreVolume() const
+ { return poreVolume_; }
+
+ /*!
+ * \brief Returns the surface tension.
+ */
+ Scalar surfaceTension() const
+ { return surfaceTension_; }
+
+protected:
+ Scalar poreInscribedRadius_;
+ Scalar poreVolume_;
+ Scalar surfaceTension_;
+};
+
+} // end namespace Dumux::PoreNetwork
+
+#endif
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>