diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index aee63d335e022d3bfee25089bbc6b30891132c7c..6fb5a53ff6429d3c00264ba2d283060a080720a2 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -124,6 +124,9 @@ NEW_PROP_TAG(UseKelvinEquation); //!< If we use Kelvin equatio
// or at the scvf center for analytical permeability fields (e.g. convergence studies)
NEW_PROP_TAG(EvaluatePermeabilityAtScvfIP);
+// When using the box method in a multi-phase context, an interface solver might be necessary
+NEW_PROP_TAG(EnableBoxInterfaceSolver);
+
//////////////////////////////////////////////////////////////
// Additional properties used by the 2pnc and 2pncmin models:
//////////////////////////////////////////////////////////////
diff --git a/dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh b/dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4c50018aedbfbbe71b542a746a11d66dd1b75065
--- /dev/null
+++ b/dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh
@@ -0,0 +1,116 @@
+// -*- 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 TwoPModel
+ * \brief copydoc Dumux::TwoPScvSaturationReconstruction
+ */
+#ifndef DUMUX_2P_BOX_MATERIAL_INTERFACE_PARAMS_HH
+#define DUMUX_2P_BOX_MATERIAL_INTERFACE_PARAMS_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/box/elementsolution.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TwoPModel
+ * \brief Class that determines the material with the lowest capillary
+ * pressure (under fully water-saturated conditions) around the nodes
+ * of a grid. These parameters are then associated with the global degree
+ * of freedom. On the basis of these parameters, the saturations in the
+ * remaining sub-control volumes connected to the vertex can be reconstructed.
+ */
+template<class SpatialParams>
+class BoxMaterialInterfaceParams
+{
+public:
+ using MaterialLawParams = typename SpatialParams::MaterialLaw::Params;
+
+ /*!
+ * \brief Update the scv -> dofparameter map
+ *
+ * \param The finite volume grid geometry
+ * \param SpatialParams Class encapsulating the spatial parameters
+ * \param x The current state of the solution vector
+ */
+ template<class FVGridGeometry, class SolutionVector>
+ void update(const FVGridGeometry& fvGridGeometry,
+ const SpatialParams& spatialParams,
+ const SolutionVector& x)
+ {
+ using MaterialLaw = typename SpatialParams::MaterialLaw;
+
+ // Make sure the spatial params return a const ref and no copy!
+ using Elem = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+ using ElemSol = decltype( elementSolution(Elem(), x, fvGridGeometry) );
+ using Scv = typename FVGridGeometry::SubControlVolume;
+ using ReturnType = decltype(spatialParams.materialLawParams(Elem(), Scv(), ElemSol()));
+ static_assert(std::is_same<ReturnType, const MaterialLawParams&>::value,
+ "In order to use the box-interface solver please provide access "
+ "to the material law parameters via returning (const) references");
+
+ // make sure this is only called for geometries of the box method!
+ if (FVGridGeometry::discMethod != DiscretizationMethod::box)
+ DUNE_THROW(Dune::InvalidStateException, "Determination of the interface material parameters with "
+ "this class only makes sense when using the box method!");
+
+ isOnMaterialInterface_.resize(fvGridGeometry.numDofs(), true);
+ dofParams_.resize(fvGridGeometry.numDofs(), nullptr);
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ const auto elemSol = elementSolution(element, x, fvGridGeometry);
+
+ auto fvGeometry = localView(fvGridGeometry);
+ fvGeometry.bind(element);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto& params = spatialParams.materialLawParams(element, scv, elemSol);
+
+ // is no parameters had been set, set them now
+ if (dofParams_[scv.dofIndex()] == nullptr)
+ dofParams_[scv.dofIndex()] = ¶ms;
+
+ // otherwise only use the current ones if endPointPc (e.g. Brooks-Corey entry pressure) is lower
+ else if (MaterialLaw::endPointPc( params ) < MaterialLaw::endPointPc( *(dofParams_[scv.dofIndex()]) ))
+ dofParams_[scv.dofIndex()] = ¶ms;
+ }
+ }
+ }
+
+ //! Return if this scv is connected to a material interface
+ template<class Scv>
+ bool isOnMaterialInterface(const Scv& scv) const
+ { return isOnMaterialInterface_[scv.dofIndex()]; }
+
+ //! Return the material parameters associated with the dof
+ template<class Scv>
+ const MaterialLawParams& getDofParams(const Scv& scv) const
+ { return *(dofParams_[scv.dofIndex()]); }
+
+private:
+ std::vector<bool> isOnMaterialInterface_;
+ std::vector<const MaterialLawParams*> dofParams_;
+};
+
+}
+
+#endif
diff --git a/dumux/porousmediumflow/2p/model.hh b/dumux/porousmediumflow/2p/model.hh
index 1aacdb680e02b1c7f0269bbb7d5ca5e21d2d9426..c6bfff71e0dd89cac8b6e911c2b1b72ea93c30e7 100644
--- a/dumux/porousmediumflow/2p/model.hh
+++ b/dumux/porousmediumflow/2p/model.hh
@@ -75,6 +75,7 @@
#include "indices.hh"
#include "volumevariables.hh"
#include "vtkoutputfields.hh"
+#include "saturationreconstruction.hh"
namespace Dumux
{
@@ -108,8 +109,10 @@ struct TwoPModelTraits
* \tparam FST The fluid state type
* \tparam PT The type used for permeabilities
* \tparam MT The model traits
+ * \tparam SR The class used for reconstruction of
+ * non-wetting phase saturations in scvs
*/
-template<class PV, class FSY, class FST, class PT, class MT>
+template<class PV, class FSY, class FST, class PT, class MT, class SR>
struct TwoPVolumeVariablesTraits
{
using PrimaryVariables = PV;
@@ -117,6 +120,7 @@ struct TwoPVolumeVariablesTraits
using FluidState = FST;
using PermeabilityType = PT;
using ModelTraits = MT;
+ using SaturationReconstruction = SR;
};
////////////////////////////////
@@ -160,7 +164,12 @@ private:
using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using PT = typename GET_PROP_TYPE(TypeTag, SpatialParams)::PermeabilityType;
- using Traits = TwoPVolumeVariablesTraits<PV, FSY, FST, PT, MT>;
+ static constexpr auto DM = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod;
+ static constexpr bool enableIS = GET_PROP_VALUE(TypeTag, EnableBoxInterfaceSolver);
+ // class used for scv-wise reconstruction of non-wetting phase saturations
+ using SR = TwoPScvSaturationReconstruction<DM, enableIS>;
+
+ using Traits = TwoPVolumeVariablesTraits<PV, FSY, FST, PT, MT, SR>;
public:
using type = TwoPVolumeVariables<Traits>;
};
diff --git a/dumux/porousmediumflow/2p/saturationreconstruction.hh b/dumux/porousmediumflow/2p/saturationreconstruction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..54ddf4bec5f019fd1e4f6e8d13e8cee3e9efd6a9
--- /dev/null
+++ b/dumux/porousmediumflow/2p/saturationreconstruction.hh
@@ -0,0 +1,105 @@
+// -*- 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 TwoPModel
+ * \brief copydoc Dumux::TwoPScvSaturationReconstruction
+ */
+#ifndef DUMUX_2P_SCV_SATURATION_RECONSTRUCTION_HH
+#define DUMUX_2P_SCV_SATURATION_RECONSTRUCTION_HH
+
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TwoPModel
+ * \brief Class that computes the non-wetting saturation in an scv from the saturation
+ * at the global degree of freedom. This is only necessary in conjunction with
+ * the box scheme where the degrees of freedom lie on material interfaces. There
+ * the non-wetting phase saturation is generally discontinuous.
+ */
+template<DiscretizationMethod M, bool enableReconstruction>
+class TwoPScvSaturationReconstruction
+{
+public:
+ /*!
+ * \brief Compute the non-wetting phase saturation.
+ * \note In the general case, we don't do anything. We do
+ * stuff only in the specialization for the box scheme below.
+ *
+ * \param SpatialParams Class encapsulating the spatial parameters
+ * \param element The finite element the scv is embedded in
+ * \param scv The sub-control volume for which the saturation is computed
+ * \param elemSol The solution at all dofs inside this element
+ * \param Sn The non-wetting phase saturation at the global dof
+ */
+ template<class SpatialParams, class Element, class Scv, class ElemSol>
+ static typename ElemSol::PrimaryVariables::value_type
+ reconstructSn(const SpatialParams& spatialParams,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol,
+ typename ElemSol::PrimaryVariables::value_type Sn)
+ { return Sn; }
+};
+
+//! Specialization for the box scheme with the interface solver enabled
+template<>
+class TwoPScvSaturationReconstruction<DiscretizationMethod::box, /*enableReconstruction*/true>
+{
+public:
+ /*!
+ * \brief Compute
+ *
+ * \param SpatialParams Class encapsulating the spatial parameters
+ * \param element The finite element the scv is embedded in
+ * \param scv The sub-control volume for which the saturation is computed
+ * \param elemSol The solution at all dofs inside this element
+ * \param Sn The non-wetting phase saturation at the global dof
+ */
+ template<class SpatialParams, class Element, class Scv, class ElemSol>
+ static typename ElemSol::PrimaryVariables::value_type
+ reconstructSn(const SpatialParams& spatialParams,
+ const Element& element,
+ const Scv& scv,
+ const ElemSol& elemSol,
+ typename ElemSol::PrimaryVariables::value_type Sn)
+ {
+ // if this dof doesn't lie on a material interface, simply return Sn
+ if (!spatialParams.materialInterfaceParams().isOnMaterialInterface(scv))
+ return Sn;
+
+ using MaterialLaw = typename SpatialParams::MaterialLaw;
+ // compute capillary pressure using material parameters associated with the dof
+ const auto& ifMaterialParams = spatialParams.materialInterfaceParams().getDofParams(scv);
+ const auto pc = MaterialLaw::pc(ifMaterialParams, /*Sw=*/1.0 - Sn);
+
+ // reconstruct by inverting the pc-sw curve
+ const auto& materialLawParams = spatialParams.materialLawParams(element, scv, elemSol);
+ const auto pcMin = MaterialLaw::endPointPc(materialLawParams);
+
+ if (pc < pcMin && pcMin > 0.0) return 0.0;
+ else return 1.0 - MaterialLaw::sw(materialLawParams, pc);
+ }
+};
+
+}
+
+#endif
diff --git a/dumux/porousmediumflow/2p/volumevariables.hh b/dumux/porousmediumflow/2p/volumevariables.hh
index 2675a44511882cf48f51703e0dfd3ca6cb60987c..f24e5fdaa94250364b93e70bc8fdda50d2272763 100644
--- a/dumux/porousmediumflow/2p/volumevariables.hh
+++ b/dumux/porousmediumflow/2p/volumevariables.hh
@@ -44,7 +44,6 @@ class TwoPVolumeVariables
using PermeabilityType = typename Traits::PermeabilityType;
using ModelTraits = typename Traits::ModelTraits;
using Indices = typename ModelTraits::Indices;
- static constexpr auto formulation = ModelTraits::priVarFormulation();
using Scalar = typename Traits::PrimaryVariables::value_type;
using FS = typename Traits::FluidSystem;
@@ -57,6 +56,7 @@ class TwoPVolumeVariables
phase1Idx = FS::phase1Idx
};
+ static constexpr auto formulation = ModelTraits::priVarFormulation();
public:
//! export type of fluid system
@@ -130,28 +130,44 @@ public:
const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
if (formulation == TwoPFormulation::p0s1)
- {
- fluidState.setSaturation(phase1Idx, priVars[saturationIdx]);
- fluidState.setSaturation(phase0Idx, 1 - priVars[saturationIdx]);
- }
- else if (formulation == TwoPFormulation::p1s0)
- {
- fluidState.setSaturation(phase0Idx, priVars[saturationIdx]);
- fluidState.setSaturation(phase1Idx, 1 - priVars[saturationIdx]);
- }
-
- pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
- if (formulation == TwoPFormulation::p0s1)
{
fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
- fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
- : priVars[pressureIdx] - pc_);
+ if (wPhaseIdx == phase1Idx)
+ {
+ fluidState.setSaturation(phase1Idx, priVars[saturationIdx]);
+ fluidState.setSaturation(phase0Idx, 1 - priVars[saturationIdx]);
+ pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
+ fluidState.setPressure(phase1Idx, priVars[pressureIdx] - pc_);
+ }
+ else
+ {
+ const auto Sn = Traits::SaturationReconstruction::reconstructSn(problem.spatialParams(), element,
+ scv, elemSol, priVars[saturationIdx]);
+ fluidState.setSaturation(phase1Idx, Sn);
+ fluidState.setSaturation(phase0Idx, 1 - Sn);
+ pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
+ fluidState.setPressure(phase1Idx, priVars[pressureIdx] + pc_);
+ }
}
else if (formulation == TwoPFormulation::p1s0)
{
fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
- fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
- : priVars[pressureIdx] + pc_);
+ if (wPhaseIdx == phase1Idx)
+ {
+ const auto Sn = Traits::SaturationReconstruction::reconstructSn(problem.spatialParams(), element,
+ scv, elemSol, priVars[saturationIdx]);
+ fluidState.setSaturation(phase0Idx, Sn);
+ fluidState.setSaturation(phase1Idx, 1 - Sn);
+ pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
+ fluidState.setPressure(phase0Idx, priVars[pressureIdx] + pc_);
+ }
+ else
+ {
+ fluidState.setSaturation(phase1Idx, priVars[saturationIdx]);
+ fluidState.setSaturation(phase0Idx, 1 - priVars[saturationIdx]);
+ pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
+ fluidState.setPressure(phase0Idx, priVars[pressureIdx] - pc_);
+ }
}
typename FluidSystem::ParameterCache paramCache;
diff --git a/dumux/porousmediumflow/properties.hh b/dumux/porousmediumflow/properties.hh
index 6ef4ddc262c7c287a08373fe7b0a60e400192cfe..10a5765aee101f6051624a3214728c5341e3473b 100644
--- a/dumux/porousmediumflow/properties.hh
+++ b/dumux/porousmediumflow/properties.hh
@@ -79,6 +79,9 @@ SET_TYPE_PROP(PorousMediumFlow, VelocityOutput, PorousMediumFlowVelocityOutput<T
SET_TYPE_PROP(PorousMediumFlow, PrimaryVariableSwitch, NoPrimaryVariableSwitch);
SET_BOOL_PROP(PorousMediumFlow, EnableThermalNonEquilibrium, false);
+
+//! Per default, we disable the box interface solver
+SET_BOOL_PROP(PorousMediumFlow, EnableBoxInterfaceSolver, false);
} // namespace Properties
} // namespace Dumux