OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Refactor the twophase laws

Browse Source 
They are now all generic capillary pressure laws, but with some
additional methods to make working with them easier. For this reason,
they have been moved up one directory to opm/material/fluidmatrixinteractions

Also, a unit test which ensures that all capillary pressure laws
conform to their respective APIs has been added.
This commit is contained in:
Andreas Lauser
2013-11-06 14:32:22 +01:00
parent 332c2f64b3
commit 7c0ff8309f
30 changed files with 2341 additions and 1639 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@@ -30,6 +30,7 @@
# originally generated with the command:
# find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
list (APPEND TEST_SOURCE_FILES
tests/test_fluidmatrixinteractions.cpp
tests/material/pengrobinson/test_pengrobinson.cpp
tests/material/ncpflash/test_ncpflash.cpp
tests/material/tabulation/test_tabulation.cpp

206
opm/material/fluidmatrixinteractions/2p/BrooksCorey.hpp
View File

@@ -1,206 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2008-2013 by Andreas Lauser *
* Copyright (C) 2010 by Philipp Nuske *
* *
* 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
* \copydoc Opm::BrooksCorey
*/
#ifndef OPM_BROOKS_COREY_HH
#define OPM_BROOKS_COREY_HH
#include "BrooksCoreyParams.hpp"
#include <algorithm>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Implementation of the Brooks-Corey capillary pressure <->
* saturation relation.
*
* This class provides the "raw" curves as static members and doesn't
* concern itself converting absolute to effective saturations and
* vice versa.
*
*\see BrooksCoreyParams
*/
template <class ScalarT, class ParamsT = BrooksCoreyParams<ScalarT> >
class BrooksCorey
{
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
/*!
* \brief The capillary pressure-saturation curve according to
* Brooks and Corey.
*
* The empirical Brooks-Corey capillary pressure-saturation
* function is defined as
* \f[
* p_C = p_e\overline{S}_w^{-1/\lambda}
* \f]
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar pC(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return params.pe()*std::pow(Sw, -1.0/params.lambda());
}
/*!
* \brief The saturation-capillary pressure curve according to
* Brooks & Corey.
*
* This is the inverse of the capillary pressure-saturation curve:
* \f[
\overline{S}_w = (\frac{p_C}{p_e})^{-\lambda}
\f]
*
* \param pc Capillary pressure \f$[Pa]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar Sw(const Params &params, Scalar pc)
{
assert(pc > 0); // if we don't assume that, std::pow will screw up!
return std::pow(pc/params.pe(), -params.lambda());
}
/*!
* \brief The partial derivative of the capillary
* pressure w.r.t. the effective saturation according to Brooks & Corey.
*
* This is equivalent to
* \f[
\frac{\partial p_C}{\partial \overline{S}_w} =
-\frac{p_e}{\lambda} \overline{S}_w^{-1/\lambda - 1}
\f]
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar dpC_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return - params.pe()/params.lambda() * std::pow(Sw, -1/params.lambda() - 1);
}
/*!
* \brief The partial derivative of the effective saturation with
* regard to the capillary pressure according to Brooks and
* Corey.
*
* \param pC Capillary pressure \f$[Pa]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
assert(pC > 0); // required for std::pow
return -params.lambda()/params.pe() * std::pow(pC/params.pe(), - params.lambda() - 1);
}
/*!
* \brief The relative permeability for the wetting phase of
* the medium implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar krw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return std::pow(Sw, 2.0/params.lambda() + 3);
}
/*!
* \brief The derivative of the relative permeability for the
* wetting phase with regard to the wetting saturation of the
* medium implied by the Brooks-Corey parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar dkrw_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return (2.0/params.lambda() + 3)*std::pow(Sw, 2.0/params.lambda() + 2);
}
/*!
* \brief The relative permeability for the non-wetting phase of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar krn(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
Scalar exponent = 2.0/params.lambda() + 1;
Scalar Sn = 1. - Sw;
return Sn*Sn*(1. - std::pow(Sw, exponent));
}
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar dkrn_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return
2.0*(Sw - 1)*(
1 +
std::pow(Sw,
2.0/params.lambda())*(
1.0/params.lambda() + 1.0/2 -
Sw*(1.0/params.lambda() + 1.0/2)
));
}
};
} // namespace Opm
#endif // BROOKS_COREY_HH

247
opm/material/fluidmatrixinteractions/2p/EffToAbsLaw.hpp
View File

@@ -1,247 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2010-2011 by Philipp Nuske *
* Copyright (C) 2009-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::EffToAbsLaw
*/
#ifndef OPM_EFF_TO_ABS_LAW_HH
#define OPM_EFF_TO_ABS_LAW_HH
#include "EffToAbsLawParams.hpp"
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief This material law takes a material law defined for effective
* saturations and converts it to a material law defined on absolute
* saturations.
*
* The idea: "material laws" (like VanGenuchten or BrooksCorey) are
* defined for effective saturations. The numeric calculations
* however are performed with absolute saturations. The EffToAbsLaw
* class gets the "material laws" actually used as well as the
* corresponding parameter container as template arguments.
*
* Subsequently, the desired function (pc, Sw... ) of the actually
* used "material laws" are called but with the saturations already
* converted from absolute to effective.
*
* This approach makes sure that in the "material laws" only effective
* saturations are considered, which makes sense, as these laws only
* deal with effective saturations. This also allows for changing the
* calculation of the effective saturations easily, as this is subject
* of discussion / may be problem specific.
*
* Additionally, handing over effective saturations to the "material
* laws" in stead of them calculating effective saturations prevents
* accidently "converting twice".
*
* This boils down to:
* - the actual material laws (linear, VanGenuchten...) do not need to
* deal with any kind of conversion
* - the definition of the material law in the spatial parameters is
* not really intuitive, but using it is: Hand in values, get back
* values, do not deal with conversion.
*/
template <class EffLawT, class AbsParamsT = EffToAbsLawParams<typename EffLawT::Params> >
class EffToAbsLaw
{
typedef EffLawT EffLaw;
public:
typedef AbsParamsT Params;
typedef typename EffLaw::Scalar Scalar;
/*!
* \brief The capillary pressure-saturation curve.
*
*
* \param Sw Absolute saturation of the wetting phase
* \f$\overline{S}_w\f$. It is converted to effective
* saturation and then handed over to the material law
* actually used for calculation.
* \param params A object that stores the appropriate coefficients
* for the respective law.
*
* \return Capillary pressure [Pa] calculated by specific
* constitutive relation (e.g. Brooks & Corey, van
* Genuchten, linear...)
*/
static Scalar pC(const Params &params, Scalar Sw)
{ return EffLaw::pC(params, SwToSwe(params, Sw));}
/*!
* \brief The saturation-capillary pressure curve.
*
* \param pC Capillary pressure \f$p_C\f$:
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*\return Absolute wetting phase saturation calculated as inverse of (EffLaw e.g. Brooks & Corey, van Genuchten, linear...) constitutive relation.
*
* \return The absolute saturation of the wetting phase \f$S_w\f$
*/
static Scalar Sw(const Params &params, Scalar pC)
{
return SweToSw_(params, EffLaw::Sw(params, pC));
}
/*!
* \brief Returns the partial derivative of the capillary
* pressure w.r.t the absolute saturation.
*
* In this case the chain rule needs to be applied:
\f[
p_c = p_c( \overline S_w (S_w))
\rightarrow p_c ^\prime = \frac{\partial p_c}{\partial \overline S_w} \frac{\partial \overline S_w}{\partial S_w}
\f]
* \param Sw Absolute saturation of the wetting phase \f$\overline{S}_w\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of \f$p_c\f$ w.r.t. effective saturation according to EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar dpC_dSw(const Params &params, Scalar Sw)
{
return EffLaw::dpC_dSw(params, SwToSwe(params, Sw) )*dSwe_dSw_(params);
}
/*!
* \brief Returns the partial derivative of the absolute
* saturation w.r.t. the capillary pressure.
*
* In this case the chain rule needs to be applied:
\f[
S_w = S_w(\overline{S}_w (p_c) )
\rightarrow S_w^\prime = \frac{\partial S_w}{\partial \overline S_w} \frac{\partial \overline S_w}{\partial p_c}
\f]
*
*
* \param pC Capillary pressure \f$p_C\f$:
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of effective saturation w.r.t. \f$p_c\f$ according to EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
return EffLaw::dSw_dpC(params, pC)*dSw_dSwe_(params);
}
/*!
* \brief The relative permeability for the wetting phase.
*
* \param Sw Absolute saturation of the wetting phase \f$\overline{S}_w\f$. It is converted to effective saturation
* and then handed over to the material law actually used for calculation.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the wetting phase calculated as implied by EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*
*/
static Scalar krw(const Params &params, Scalar Sw)
{
return EffLaw::krw(params, SwToSwe(params, Sw));
}
/*!
* \brief The relative permeability for the non-wetting phase.
*
* \param Sw Absolute saturation of the wetting phase \f${S}_w\f$. It is converted to effective saturation
* and then handed over to the material law actually used for calculation.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the non-wetting phase calculated as implied by EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar krn(const Params &params, Scalar Sw)
{
return EffLaw::krn(params, SwToSwe(params, Sw));
}
/*!
* \brief Convert an absolute wetting saturation to an effective one.
*
* \param Sw Absolute saturation of the wetting phase \f${S}_w\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Effective saturation of the wetting phase.
*/
static Scalar SwToSwe(const Params &params, Scalar Sw)
{
return (Sw - params.Swr())/(1 - params.Swr() - params.Snr());
}
/*!
* \brief Convert an absolute non-wetting saturation to an effective one.
*
* \param Sn Absolute saturation of the non-wetting phase \f${S}_n\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Effective saturation of the non-wetting phase.
*/
static Scalar SnToSne(const Params &params, Scalar Sn)
{
return (Sn - params.Snr())/(1 - params.Swr() - params.Snr());
}
private:
/*!
* \brief Convert an effective wetting saturation to an absolute one.
*
* \param Swe Effective saturation of the non-wetting phase \f$\overline{S}_n\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Absolute saturation of the non-wetting phase.
*/
static Scalar SweToSw_(const Params &params, Scalar Swe)
{
return Swe*(1 - params.Swr() - params.Snr()) + params.Swr();
}
/*!
* \brief Derivative of the effective saturation w.r.t. the absolute saturation.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the effective saturation w.r.t. the absolute saturation.
*/
static Scalar dSwe_dSw_(const Params &params)
{ return 1.0/(1 - params.Swr() - params.Snr()); }
/*!
* \brief Derivative of the absolute saturation w.r.t. the effective saturation.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the absolute saturation w.r.t. the effective saturation.
*/
static Scalar dSw_dSwe_(const Params &params)
{ return 1 - params.Swr() - params.Snr(); }
};
} // namespace Opm
#endif

172
opm/material/fluidmatrixinteractions/2p/LinearMaterial.hpp
View File

@@ -1,172 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2008-2012 by Andreas Lauser *
* Copyright (C) 2010 by Philipp Nuske *
* *
* 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
* \copydoc Opm::LinearMaterial
*/
#ifndef OPM_LINEAR_MATERIAL_HH
#define OPM_LINEAR_MATERIAL_HH
#include "LinearMaterialParams.hpp"
#include <algorithm>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Linear capillary pressure and relative permeability <->
* saturation relations
*
* The entry pressure is reached at \f$ \overline S_w = 1\f$, the maximum
* capillary pressure is observed at \f$ \overline S_w = 0\f$.
*
* For general info see EffToAbsLaw
*
* \see LinearMaterialParams
*/
template <class ScalarT, class ParamsT = LinearMaterialParams<ScalarT> >
class LinearMaterial
{
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
/*!
* \brief The linear capillary pressure-saturation curve.
*
* This material law is linear:
* \f[
p_C = (1 - \overline{S}_w) (p_{C,max} - p_{C,entry}) + p_{C,entry}
\f]
*
* \param Swe Effective saturation of the wetting phase
* \f$\overline{S}_w\f$ conversion from absolute
* saturation happened in EffToAbsLaw.
* \param params A object that provides the appropriate coefficients
* for the respective law.
*
* \return Capillary pressure calculated by linear constitutive
* relation.
*/
static Scalar pC(const Params &params, Scalar Swe)
{
return (1 - Swe)*(params.maxPC() - params.entryPC()) + params.entryPC();
}
/*!
* \brief The saturation-capillary pressure curve.
*
* This is the inverse of the capillary pressure-saturation curve:
* \f[
S_w = 1 - \frac{p_C - p_{C,entry}}{p_{C,max} - p_{C,entry}}
\f]
*
* \param pC Capillary pressure \f$p_C\f$
* \param params A container object that is populated with the
* appropriate coefficients for the respective law.
*
* \return Effective wetting phase saturation calculated as
* inverse of the linear constitutive relation.
*/
static Scalar Sw(const Params &params, Scalar pC)
{
return 1 - (pC - params.entryPC())/(params.maxPC() - params.entryPC());
}
/*!
* \brief Returns the partial derivative of the capillary
* pressure w.r.t. the effective saturation.
*
* This is equivalent to
* \f[
\frac{\partial p_C}{\partial \overline{S}_w} =
- (p_{C,max} - p_{C,min})
\f]
*
* \param Swe Effective saturation of the wetting phase
* \f$\overline{S}_w\f$ conversion from absolute
* saturation happened in EffToAbsLaw.
* \param params A object that provides the appropriate coefficients
* for the respective law.
*
* \return Partial derivative of \f$p_c\f$ w.r.t. effective
* saturation according to linear material relation.
*/
static Scalar dpC_dSw(const Params &params, Scalar Swe)
{
return - (params.maxPC() - params.entryPC());
}
/*!
* \brief Returns the partial derivative of the effective
* saturation to the capillary pressure.
*
* \param pC Capillary pressure \f$p_C\f$ [Pa]
* \param params A object that provides the appropriate coefficients
* for the respective law.
*
* \return Partial derivative of effective saturation
* w.r.t. \f$p_c\f$ according to linear relation.
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
return - 1/(params.maxPC() - params.entryPC());
}
/*!
* \brief The relative permeability for the wetting phase.
*
* \param params A object that provides the appropriate coefficients
* for the respective law.
* \param Swe Effective saturation of the wetting phase
* \f$\overline{S}_w\f$. The conversion from absolute
* saturation is done by EffToAbsLaw.
*
* \return Relative permeability of the wetting phase calculated
* as linear relation.
*/
static Scalar krw(const Params &params, Scalar Swe)
{
return std::max(std::min(Swe,1.0),0.0);
}
/*!
* \brief The relative permeability for the non-wetting phase.
*
*
* \param params A object that provides the appropriate coefficients
* for the respective law.
* \param Swe Effective saturation of the wetting phase
* \f$\overline{S}_w\f$. The conversion from absolute
* saturation is done by EffToAbsLaw.
*
* \return Relative permeability of the non-wetting phase
* calculated as linear relation.
*/
static Scalar krn(const Params &params, Scalar Swe)
{
Scalar Sne = 1 - Swe;
return std::max(std::min(Sne,1.0),0.0);
}
};
} // namespace Opm
#endif

92
opm/material/fluidmatrixinteractions/2p/LinearMaterialParams.hpp
View File

@@ -1,92 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2009-2012 by Andreas Lauser *
* Copyright (C) 2010 by Philipp Nuske *
* *
* 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
* \copydoc Opm::LinearMaterialParams
*/
#ifndef LINEAR_MATERIAL_PARAMS_HH
#define LINEAR_MATERIAL_PARAMS_HH
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Reference implementation of the parameter object for the
* linear material law.
*/
template<class ScalarT>
class LinearMaterialParams
{
public:
typedef ScalarT Scalar;
LinearMaterialParams()
{
setEntryPC(0);
setMaxPC(0);
}
LinearMaterialParams(Scalar entryPC, Scalar maxPC)
{
setEntryPC(entryPC);
setMaxPC(maxPC);
}
/*!
* \brief Return the entry pressure for the linear material law.
*
* The entry pressure is reached at \f$\overline S_w = 1\f$
*/
Scalar entryPC() const
{ return entryPC_; }
/*!
* \brief Set the entry pressure for the linear material law.
*
* The entry pressure is reached at \f$ \overline S_w = 1\f$
*/
void setEntryPC(Scalar v)
{ entryPC_ = v; }
/*!
* \brief Return the maximum capillary pressure for the linear material law.
*
* The maximum capillary pressure is reached at \f$ \overline S_w = 0\f$
*/
Scalar maxPC() const
{ return maxPC_; }
/*!
* \brief Set the maximum capillary pressure for the linear material law.
*
* The maximum capillary pressure is reached at \f$ \overline S_w = 0\f$
*/
void setMaxPC(Scalar v)
{ maxPC_ = v; }
private:
Scalar entryPC_;
Scalar maxPC_;
};
} // namespace Opm
#endif

197
opm/material/fluidmatrixinteractions/2p/RegularizedLinearMaterial.hpp
View File

@@ -1,197 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2010 by Markus Wolff *
* Copyright (C) 2010 by Philipp Nuske *
* Copyright (C) 2008-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::RegularizedLinearMaterial
*/
#ifndef OPM_REGULARIZED_LINEAR_MATERIAL_HH
#define OPM_REGULARIZED_LINEAR_MATERIAL_HH
#include "LinearMaterial.hpp"
#include "RegularizedLinearMaterialParams.hpp"
#include <opm/core/utility/Spline.hpp>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
* \brief Implements a linear saturation-capillary pressure relation
*
* The entry pressure is reached at \f$\overline S_w = 1\f$, the maximum
* capillary pressure is observed at \f$\overline S_w = 0\f$.
*
* It may seem strange to regularize a linear material, here comes the rationale:
*
* The relative permeabilities are 0 or 1 outside of the range of effective saturation.
* However, the transition between the linearly changing and the constant part is not smooth but with a kink.
* The Newton scheme does not like that. Therefore a smooth transition is accomplished by interpolating these
* regions with a spline.
*
* An example of the regularization of the relative permeability is shown below:
* \image html regularizedLinearKr.png
*
*
* \see LinearMaterial
*/
template <class ScalarT, class ParamsT = RegularizedLinearMaterialParams<ScalarT> >
class RegularizedLinearMaterial
{
typedef Opm::LinearMaterial<ScalarT, ParamsT> LinearMaterial;
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
/*!
* \brief The linear capillary pressure-saturation curve.
*
* This material law is linear:
* \f[
p_C = (1 - \overline{S}_w) (p_{C,max} - p_{C,entry}) + p_{C,entry}
\f]
*
* \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$ conversion from absolute saturation happened in EffToAbsLaw.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar pC(const Params &params, Scalar Swe)
{
return LinearMaterial::pC(params, Swe);
}
/*!
* \brief The saturation-capillary pressure curve.
*
* This is the inverse of the capillary pressure-saturation curve:
* \f[
S_w = 1 - \frac{p_C - p_{C,entry}}{p_{C,max} - p_{C,entry}}
\f]
*
* \param pC Capillary pressure \f$p_C\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return The effective saturation of the wetting phase \f$\overline{S}_w\f$
*/
static Scalar Sw(const Params &params, Scalar pC)
{
return LinearMaterial::Sw(params, pC);
}
/*!
* \brief Returns the partial derivative of the capillary
* pressure to the effective saturation.
*
* This is equivalent to
* \f[
\frac{\partial p_C}{\partial \overline{S}_w} =
- (p_{C,max} - p_{C,min})
\f]
* \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$ conversion from absolute saturation happened in EffToAbsLaw.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar dpC_dSw(const Params &params, Scalar Swe)
{
return LinearMaterial::dpC_dSw(params, Swe);
}
/*!
* \brief Returns the partial derivative of the effective
* saturation to the capillary pressure.
*
* \param pC Capillary pressure \f$p_C\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
return LinearMaterial::dSw_dpC(params, pC);
}
/*!
* \brief The relative permeability for the wetting phase.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$ conversion from absolute saturation happened in EffToAbsLaw.
*/
static Scalar krw(const Params &params, Scalar Swe)
{
return relperm_(params, Swe);
}
/*!
* \brief The relative permeability for the non-wetting phase.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$ conversion from absolute saturation happened in EffToAbsLaw.
*/
static Scalar krn(const Params &params, Scalar Swe)
{
Scalar Sne = 1 - Swe;
return relperm_(params, Sne);
}
private:
static Scalar relperm_(const Params &params, Scalar S)
{
const Scalar lowS = params.krLowS();
const Scalar highS = params.krHighS();
const Scalar m = (1 - ((1 - highS) + lowS)/2 ) / (1 - (1 - highS) - lowS);
// check whether the saturation is unpyhsical
if (S >= 1.0)
return 1.0;
else if (S <= 0.0)
return 0;
// check wether the permeability needs to be regularized
else if (S < lowS) {
typedef Opm::Spline<Scalar> Spline;
Spline sp(0, lowS,
0, lowS/2,
0, m);
return sp.eval(S);
}
else if (S > highS) {
typedef Opm::Spline<Scalar> Spline;
Spline sp(highS, 1,
1 - (1 - highS)/2, 1,
m, 0);
return sp.eval(S);
}
// straight line for S \in [lowS, highS]
return lowS/2 + m*(S - lowS);
}
};
} // namespace Opm
#endif

68
opm/material/fluidmatrixinteractions/2p/RegularizedLinearMaterialParams.hpp
View File

@@ -1,68 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2010 by Philipp Nuske *
* *
* 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
* \copydoc Opm::RegularizedLinearMaterialParams
*/
#ifndef REGULARIZED_LINEAR_PARAMS_HH
#define REGULARIZED_LINEAR_PARAMS_HH
#include "LinearMaterialParams.hpp"
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Parameters that are necessary for the \em regularization of
* the linear constitutive relations.
*
*/
template<class ScalarT>
class RegularizedLinearMaterialParams : public LinearMaterialParams<ScalarT>
{
public:
typedef ScalarT Scalar;
RegularizedLinearMaterialParams()
{}
/*!
* \brief Return the threshold saturation respective phase below
* which the relative permeability gets regularized.
*
* This is just 5%. If you need a different value, write your own
* parameter class.
*/
Scalar krLowS() const
{ return 0.05; }
/*!
* \brief Return the threshold saturation of the respective phase
* above which the relative permeability gets regularized.
*
* This is just 95%. If you need a different value, write your own
* parameter class.
*/
Scalar krHighS() const
{ return 0.95; }
};
} // namespace Opm
#endif

109
opm/material/fluidmatrixinteractions/2p/RegularizedVanGenuchtenParams.hpp
View File

@@ -1,109 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2008-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::RegularizedVanGenuchtenParams
*/
#ifndef OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
#define OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
#include "VanGenuchtenParams.hpp"
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Parameters that are necessary for the \em regularization of
* VanGenuchten "material law".
*
*/
template<class ScalarT>
class RegularizedVanGenuchtenParams : public VanGenuchtenParams<ScalarT>
{
public:
typedef ScalarT Scalar;
typedef VanGenuchtenParams<Scalar> Parent;
RegularizedVanGenuchtenParams()
: pCLowSw_(0.01)
, pCHighSw_(0.99)
{}
RegularizedVanGenuchtenParams(Scalar vgAlpha,
Scalar vgN)
: Parent(vgAlpha, vgN)
, pCLowSw_(0.01)
, pCHighSw_(0.99)
{}
/*!
* \brief Return the threshold saturation below which the
* capillary pressure is regularized.
*/
Scalar pCLowSw() const
{ return pCLowSw_; }
/*!
* \brief Set the threshold saturation below which the capillary
* pressure is regularized.
*/
void setPCLowSw(Scalar value)
{ pCLowSw_ = value; }
/*!
* \brief Return the threshold saturation below which the
* capillary pressure is regularized.
*/
Scalar pCHighSw() const
{ return pCHighSw_; }
/*!
* \brief Set the threshold saturation below which the capillary
* pressure is regularized.
*/
void setPCHighSw(Scalar value)
{ pCHighSw_ = value; }
/*!
* \brief Threshold saturation below which the relative
* permeability of the non-wetting phase gets regulatized.
*
* This is just 10%. If you need a different value, overload this
* class.
*/
Scalar krnLowSw() const
{ return 0.10; }
/*!
* \brief Threshold saturation above which the relative
* permeability of the wetting phase gets regulatized.
*
* This is just 90%. If you need a different value, overload this
* class.
*/
Scalar krwHighSw() const
{ return 0.90; }
private:
Scalar pCLowSw_;
Scalar pCHighSw_;
};
} // namespace Opm
#endif

86
opm/material/fluidmatrixinteractions/2pAdapter.hpp
View File

@@ -1,86 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2011-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::TwoPAdapter
*/
#ifndef OPM_MP_2P_ADAPTER_HH
#define OPM_MP_2P_ADAPTER_HH
#include <algorithm>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
* \brief Provides an adapter class to use two-phase material laws
* with the generic M-phase API.
*/
template <int wPhaseIdxT, class TwoPLaw >
class TwoPAdapter
{
public:
typedef typename TwoPLaw::Params Params;
enum { numPhases = 2 };
enum { wPhaseIdx = wPhaseIdxT };
enum { nPhaseIdx = (wPhaseIdx == 0)?1:0 };
/*!
* \brief The capillary pressure-saturation curve.
*/
template <class ContainerT, class FluidState>
static void capillaryPressures(ContainerT &values,
const Params &params,
const FluidState &fluidState)
{
// wetting phase does not get anything added
values[wPhaseIdx] = 0;
// non-wetting phase gets the capillary pressure added
values[nPhaseIdx] = TwoPLaw::pC(params, fluidState.saturation(wPhaseIdx));
}
/*!
* \brief The inverse of the capillary pressure-saturation curve.
*/
template <class ContainerT, class FluidState>
static void saturations(ContainerT &values,
const Params &params,
const FluidState &fluidState)
{
// wetting phase does not get anything added
values[wPhaseIdx] = TwoPLaw::Sw(params,
fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx));
values[nPhaseIdx] = 1.0 - values[wPhaseIdx];
}
/*!
* \brief The relative permeability of all phases.
*/
template <class ContainerT, class FluidState>
static void relativePermeabilities(ContainerT &values,
const Params &params,
const FluidState &fluidState)
{
values[wPhaseIdx] = TwoPLaw::krw(params, fluidState.saturation(wPhaseIdx));
values[nPhaseIdx] = TwoPLaw::krn(params, fluidState.saturation(wPhaseIdx));
}
};
} // namespace Opm
#endif

334
opm/material/fluidmatrixinteractions/BrooksCorey.hpp Normal file
View File

@@ -0,0 +1,334 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2008-2013 by Andreas Lauser *
* Copyright (C) 2010 by Philipp Nuske *
* *
* 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
* \copydoc Opm::BrooksCorey
*/
#ifndef OPM_BROOKS_COREY_HH
#define OPM_BROOKS_COREY_HH
#include "BrooksCoreyParams.hpp"
#include <algorithm>
#include <cassert>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Implementation of the Brooks-Corey capillary pressure <->
* saturation relation.
*
* This class provides the "raw" curves as static members and doesn't
* concern itself converting absolute to effective saturations and
* vice versa.
*
*\see BrooksCoreyParams
*/
template <class TraitsT, class ParamsT = BrooksCoreyParams<TraitsT> >
class BrooksCorey : public TraitsT
{
public:
typedef TraitsT Traits;
typedef ParamsT Params;
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases to which this material law applies.
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The Brooks-Corey capillary pressure law only applies "
"to the case of two fluid phases");
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
static_assert(Traits::numPhases == 2,
"The number of fluid phases must be two if you want to use "
"this material law!");
/*!
* \brief The capillary pressure-saturation curves according to van Genuchten.
*
* Van Genuchten's empirical capillary pressure <-> saturation
* function is given by
* \f[
* p_{c,wn} = p_n - p_w = ({S_w}^{-1/m} - 1)^{1/n}/\alpha
* \f]
*
* \param values A random access container which stores the
* relative pressure of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the capillary pressure
* ought to be calculated
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcwn(params, fs);
}
/*!
* \brief Calculate the saturations of the phases starting from
* their pressure differences.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
}
/*!
* \brief The relative permeability-saturation curves according to van Genuchten.
*
* \param values A random access container which stores the
* relative permeability of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the relative permeabilities
* ought to be calculated
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
}
/*!
* \brief The capillary pressure-saturation curve according to
* Brooks and Corey.
*
* The empirical Brooks-Corey capillary pressure-saturation
* function is defined as
* \f[
* p_C = p_e\overline{S}_w^{-1/\lambda}
* \f]
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
return twoPhaseSatPcwn(params, Sw);
}
static Scalar twoPhaseSatPcwn(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return params.entryPressure()*std::pow(Sw, -1.0/params.lambda());
}
/*!
* \brief The saturation-capillary pressure curve according to
* Brooks & Corey.
*
* This is the inverse of the capillary pressure-saturation curve:
* \f[
\overline{S}_w = (\frac{p_C}{p_e})^{-\lambda}
\f]
*
* \param pc Capillary pressure \f$[Pa]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pc = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
return twoPhaseSatSw(params, pc);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pc)
{
assert(pc > 0); // if we don't assume that, std::pow will screw up!
return std::pow(pc/params.entryPressure(), -params.lambda());
}
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pc)
{ return 1 - twoPhaseSatSw(params, pc); }
/*!
* \brief The partial derivative of the capillary
* pressure w.r.t. the effective saturation according to Brooks & Corey.
*
* This is equivalent to
* \f[
\frac{\partial p_C}{\partial \overline{S}_w} =
-\frac{p_e}{\lambda} \overline{S}_w^{-1/\lambda - 1}
\f]
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar dpcwn_dSw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
return twoPhaseSatDpcwn_dSw(params, Sw);
}
static Scalar twoPhaseSatDpcwn_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return - params.entryPressure()/params.lambda() * std::pow(Sw, -1/params.lambda() - 1);
}
/*!
* \brief The partial derivative of the effective saturation with
* regard to the capillary pressure according to Brooks and
* Corey.
*
* \param pcwn Capillary pressure \f$[Pa]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar dSw_dpcwn(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
return twoPhaseSatDSw_dpcwn(params, Sw);
}
static Scalar twoPhaseSatDSw_dpcwn(const Params &params, Scalar Sw)
{
assert(pcwn > 0); // required for std::pow
return -params.lambda()/params.entryPressure() * std::pow(pcwn/params.entryPressure(), - params.lambda() - 1);
}
/*!
* \brief The relative permeability for the wetting phase of
* the medium implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return std::pow(Sw, 2.0/params.lambda() + 3);
}
/*!
* \brief The derivative of the relative permeability for the
* wetting phase with regard to the wetting saturation of the
* medium implied by the Brooks-Corey parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar dkrw_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return (2.0/params.lambda() + 3)*std::pow(Sw, 2.0/params.lambda() + 2);
}
/*!
* \brief The relative permeability for the non-wetting phase of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
Scalar exponent = 2.0/params.lambda() + 1;
Scalar Sn = 1. - Sw;
return Sn*Sn*(1. - std::pow(Sw, exponent));
}
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar dkrn_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
Scalar alpha =
1.0/params.lambda() + 1.0/2 -
Sw*(1.0/params.lambda() + 1.0/2);
return 2.0*(Sw - 1)*(1 + std::pow(Sw, 2.0/params.lambda())*alpha);
}
};
} // namespace Opm
#endif // BROOKS_COREY_HH

34
opm/material/fluidmatrixinteractions/2p/BrooksCoreyParams.hpp → opm/material/fluidmatrixinteractions/BrooksCoreyParams.hpp
View File

@@ -36,33 +36,39 @@ namespace Opm {
*
*\see BrooksCorey
*/
template <class ScalarT>
template <class TraitsT>
class BrooksCoreyParams
{
typedef typename TraitsT::Scalar Scalar;
public:
typedef ScalarT Scalar;
typedef TraitsT Traits;
BrooksCoreyParams()
{
Valgrind::SetUndefined(*this);
}
{ Valgrind::SetUndefined(*this); }
BrooksCoreyParams(Scalar pe, Scalar lambda)
: pe_(pe), lambda_(lambda)
{
}
BrooksCoreyParams(Scalar entryPressure, Scalar lambda)
: entryPressure_(entryPressure), lambda_(lambda)
{ }
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{ }
/*!
* \brief Returns the entry pressure [Pa]
*/
Scalar pe() const
{ return pe_; }
Scalar entryPressure() const
{ return entryPressure_; }
/*!
* \brief Set the entry pressure [Pa]
*/
void setPe(Scalar v)
{ pe_ = v; }
void setEntryPressure(Scalar v)
{ entryPressure_ = v; }
/*!
@@ -78,7 +84,7 @@ public:
{ lambda_ = v; }
private:
Scalar pe_;
Scalar entryPressure_;
Scalar lambda_;
};
} // namespace Opm

431
opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp Normal file
View File

@@ -0,0 +1,431 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2010-2011 by Philipp Nuske *
* Copyright (C) 2009-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::EffToAbsLaw
*/
#ifndef OPM_EFF_TO_ABS_LAW_HH
#define OPM_EFF_TO_ABS_LAW_HH
#include "EffToAbsLawParams.hpp"
#include <opm/material/fluidstates/SaturationOverlayFluidState.hpp>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief This material law takes a material law defined for effective
* saturations and converts it to a material law defined on absolute
* saturations.
*
* The idea: "material laws" (like VanGenuchten or BrooksCorey) are
* defined for effective saturations. The numeric calculations
* however are performed with absolute saturations. The EffToAbsLaw
* class gets the "material laws" actually used as well as the
* corresponding parameter container as template arguments.
*
* Subsequently, the desired function (pc, Sw... ) of the actually
* used "material laws" are called but with the saturations already
* converted from absolute to effective.
*
* This approach makes sure that in the "material laws" only effective
* saturations are considered, which makes sense, as these laws only
* deal with effective saturations. This also allows for changing the
* calculation of the effective saturations easily, as this is subject
* of discussion / may be problem specific.
*
* Additionally, handing over effective saturations to the "material
* laws" in stead of them calculating effective saturations prevents
* accidently "converting twice".
*
* This boils down to:
* - the actual material laws (linear, VanGenuchten...) do not need to
* deal with any kind of conversion
* - the definition of the material law in the spatial parameters is
* not really intuitive, but using it is: Hand in values, get back
* values, do not deal with conversion.
*/
template <class EffLawT, class ParamsT = EffToAbsLawParams<typename EffLawT::Params, EffLawT::numPhases> >
class EffToAbsLaw : public EffLawT::Traits
{
typedef EffLawT EffLaw;
public:
typedef typename EffLaw::Traits Traits;
typedef ParamsT Params;
typedef typename EffLaw::Scalar Scalar;
//! The number of fluid phases
static const int numPhases = EffLaw::numPhases;
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = EffLaw::implementsTwoPhaseApi;
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = EffLaw::implementsTwoPhaseSatApi;
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = EffLaw::isSaturationDependent;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = EffLaw::isPressureDependent;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = EffLaw::isTemperatureDependent;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = EffLaw::isCompositionDependent;
/*!
* \brief The capillary pressure-saturation curves depending on absolute saturations.
*
* \param values A random access container which stores the
* relative pressure of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the capillary pressure
* ought to be calculated
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<Scalar, FluidState> OverlayFluidState;
OverlayFluidState overlayFs(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
overlayFs.setSaturation(phaseIdx,
effectiveSaturation(params,
fs.saturation(phaseIdx),
phaseIdx));
}
EffLaw::capillaryPressures(values, params, overlayFs);
}
/*!
* \brief The relative permeability-saturation curves depending on absolute saturations.
*
* \param values A random access container which stores the
* relative permeability of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the relative permeabilities
* ought to be calculated
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<Scalar, FluidState> OverlayFluidState;
OverlayFluidState overlayFs(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
overlayFs.setSaturation(phaseIdx,
effectiveSaturation(params,
fs.saturation(phaseIdx),
phaseIdx));
}
EffLaw::relativePermeabilities(values, params, overlayFs);
}
/*!
* \brief The capillary pressure-saturation curve.
*
*
* \param Sw Absolute saturation of the wetting phase
* \f$\overline{S}_w\f$. It is converted to effective
* saturation and then handed over to the material law
* actually used for calculation.
* \param params A object that stores the appropriate coefficients
* for the respective law.
*
* \return Capillary pressure [Pa] calculated by specific
* constitutive relation (e.g. Brooks & Corey, van
* Genuchten, linear...)
*/
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<Scalar, FluidState> OverlayFluidState;
static_assert(FluidState::numPhases == numPhases,
"The fluid state and the material law must exhibit the same "
"number of phases!");
OverlayFluidState overlayFs(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
overlayFs.setSaturation(phaseIdx,
effectiveSaturation(params,
fs.saturation(phaseIdx),
phaseIdx));
}
return EffLaw::pcwn(params, overlayFs);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatPcwn(const Params &params, Scalar SwAbs)
{
Scalar SwEff = effectiveSaturation(params, SwAbs, Traits::wPhaseIdx);
return EffLaw::twoPhaseSatPcwn(params, SwEff);
}
/*!
* \brief The saturation-capillary pressure curves.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
EffLaw::saturations(values, params, fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
values[phaseIdx] = absoluteSaturation(params, values[phaseIdx], phaseIdx);
}
}
/*!
* \brief Calculate wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sw(params, fs), Traits::wPhaseIdx); }
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatSw(const Params &params, Scalar Sw)
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSw(params, Sw), Traits::wPhaseIdx); }
/*!
* \brief Calculate non-wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sn(params, fs), Traits::nPhaseIdx); }
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatSn(const Params &params, Scalar Sw)
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSn(params, Sw), Traits::nPhaseIdx); }
/*!
* \brief Calculate gas phase saturation given that the rest of
* the fluid state has been initialized
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
Sg(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sg(params, fs), Traits::gPhaseIdx); }
/*!
* \brief Returns the partial derivative of the capillary
* pressure w.r.t the absolute saturation.
*
* In this case the chain rule needs to be applied:
\f[
p_c = p_c( \overline S_w (S_w))
\rightarrow p_c ^\prime = \frac{\partial p_c}{\partial \overline S_w} \frac{\partial \overline S_w}{\partial S_w}
\f]
* \param Sw Absolute saturation of the wetting phase \f$\overline{S}_w\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of \f$p_c\f$ w.r.t. effective saturation according to EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar dpC_dSw(const Params &params, Scalar Sw)
{
return EffLaw::dpC_dSw(params, SabsToSeff(params, Sw) )*dSwe_dSw_(params);
}
/*!
* \brief Returns the partial derivative of the absolute
* saturation w.r.t. the capillary pressure.
*
* In this case the chain rule needs to be applied:
\f[
S_w = S_w(\overline{S}_w (p_c) )
\rightarrow S_w^\prime = \frac{\partial S_w}{\partial \overline S_w} \frac{\partial \overline S_w}{\partial p_c}
\f]
*
*
* \param pC Capillary pressure \f$p_C\f$:
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of effective saturation w.r.t. \f$p_c\f$ according to EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
return EffLaw::dSw_dpC(params, pC)*dSw_dSwe_(params);
}
/*!
* \brief The relative permeability for the wetting phase.
*
* \param Sw Absolute saturation of the wetting phase \f$\overline{S}_w\f$. It is converted to effective saturation
* and then handed over to the material law actually used for calculation.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the wetting phase calculated as implied by EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<Scalar, FluidState> OverlayFluidState;
static_assert(FluidState::numPhases == numPhases,
"The fluid state and the material law must exhibit the same "
"number of phases!");
OverlayFluidState overlayFs(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
overlayFs.setSaturation(phaseIdx,
effectiveSaturation(params,
fs.saturation(phaseIdx),
phaseIdx));
}
return EffLaw::krw(params, overlayFs);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return EffLaw::twoPhaseSatKrw(params, effectiveSaturation(params, Sw, Traits::nPhaseIdx)); }
/*!
* \brief The relative permeability of the non-wetting phase.
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<Scalar, FluidState> OverlayFluidState;
static_assert(FluidState::numPhases == numPhases,
"The fluid state and the material law must exhibit the same "
"number of phases!");
OverlayFluidState overlayFs(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
overlayFs.setSaturation(phaseIdx,
effectiveSaturation(params,
fs.saturation(phaseIdx),
phaseIdx));
}
return EffLaw::krn(params, overlayFs);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return EffLaw::twoPhaseSatKrn(params, effectiveSaturation(params, Sw, Traits::nPhaseIdx)); }
/*!
* \brief The relative permability of the gas phase
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
krg(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<Scalar, FluidState> OverlayFluidState;
static_assert(FluidState::numPhases == numPhases,
"The fluid state and the material law must exhibit the same "
"number of phases!");
OverlayFluidState overlayFs(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
overlayFs.setSaturation(phaseIdx,
effectiveSaturation(params,
fs.saturation(phaseIdx),
phaseIdx));
}
return EffLaw::krg(params, overlayFs);
}
/*!
* \brief Convert an absolute saturation to an effective one.
*/
static Scalar effectiveSaturation(const Params &params, Scalar S, int phaseIdx)
{ return (S - params.residualSaturation(phaseIdx))/(1 - params.sumResidualSaturations()); }
/*!
* \brief Convert an effective saturation to an absolute one.
*/
static Scalar absoluteSaturation(const Params &params, Scalar S, int phaseIdx)
{ return S*(1 - params.sumResidualSaturations()) + params.residualSaturation(phaseIdx); }
private:
/*!
* \brief Convert an effective wetting saturation to an absolute one.
*
* \param Swe Effective saturation of the non-wetting phase \f$\overline{S}_n\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Absolute saturation of the non-wetting phase.
*/
static Scalar SweToSw_(const Params &params, Scalar Swe)
{
return Swe*(1 - params.Swr() - params.Snr()) + params.Swr();
}
/*!
* \brief Derivative of the effective saturation w.r.t. the absolute saturation.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the effective saturation w.r.t. the absolute saturation.
*/
static Scalar dSwe_dSw_(const Params &params)
{ return 1.0/(1 - params.Swr() - params.Snr()); }
/*!
* \brief Derivative of the absolute saturation w.r.t. the effective saturation.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the absolute saturation w.r.t. the effective saturation.
*/
static Scalar dSw_dSwe_(const Params &params)
{ return 1 - params.Swr() - params.Snr(); }
};
} // namespace Opm
#endif

48
opm/material/fluidmatrixinteractions/2p/EffToAbsLawParams.hpp → opm/material/fluidmatrixinteractions/EffToAbsLawParams.hpp
View File

@@ -31,44 +31,58 @@ namespace Opm {
* class to convert material laws from effective to absolute
* saturations.
*/
template <class EffLawParamsT>
template <class EffLawParamsT, int numPhases>
class EffToAbsLawParams : public EffLawParamsT
{
typedef EffLawParamsT EffLawParams;
typedef typename EffLawParams::Traits::Scalar Scalar;
public:
typedef typename EffLawParams::Scalar Scalar;
typedef typename EffLawParams::Traits Traits;
EffToAbsLawParams()
: EffLawParams()
{ Swr_ = Snr_ = 0; }
{
sumResidualSaturations_ = 0.0;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
residualSaturation_[phaseIdx] = 0.0;
}
/*!
* \brief Return the residual wetting saturation.
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
Scalar Swr() const
{ return Swr_; }
void finalize()
{
EffLawParams::finalize();
}
/*!
* \brief Set the residual wetting saturation.
* \brief Return the residual saturation of a phase.
*/
void setSwr(Scalar v)
{ Swr_ = v; }
Scalar residualSaturation(int phaseIdx) const
{ return residualSaturation_[phaseIdx]; }
/*!
* \brief Return the residual non-wetting saturation.
* \brief Return the sum of the residual saturations.
*/
Scalar Snr() const
{ return Snr_; }
Scalar sumResidualSaturations() const
{ return sumResidualSaturations_; }
/*!
* \brief Set the residual non-wetting saturation.
* \brief Set the residual saturation of a phase.
*/
void setSnr(Scalar v)
{ Snr_ = v; }
void setResidualSaturation(int phaseIdx, Scalar value)
{
sumResidualSaturations_ -= residualSaturation_[phaseIdx];
sumResidualSaturations_ += value;
residualSaturation_[phaseIdx] = value;
}
private:
Scalar Swr_;
Scalar Snr_;
Scalar residualSaturation_[numPhases];
Scalar sumResidualSaturations_;
};
} // namespace Opm

275
opm/material/fluidmatrixinteractions/LinearMaterial.hpp Normal file
View File

@@ -0,0 +1,275 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2011-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::LinearMaterial
*/
#ifndef OPM_LINEAR_MATERIAL_HH
#define OPM_LINEAR_MATERIAL_HH
#include "LinearMaterialParams.hpp"
#include <opm/material/Valgrind.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <type_traits>
namespace Opm {
/*!
* \ingroup material
*
* \brief Implements a linear saturation-capillary pressure relation
*
* Implements a linear saturation-capillary pressure relation for
* M-phase fluid systems.
*
* \sa LinearMaterialParams
*/
template <class TraitsT, class ParamsT = LinearMaterialParams<TraitsT> >
class LinearMaterial : public TraitsT
{
public:
typedef TraitsT Traits;
typedef ParamsT Params;
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases
static const int numPhases = Traits::numPhases;
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = (numPhases == 2);
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = (numPhases == 2);
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
/*!
* \brief The linear capillary pressure-saturation curve.
*
* This material law is linear:
* \f[
p_C = (1 - \overline{S}_w) (p_{C,max} - p_{C,entry}) + p_{C,entry}
\f]
*
* \param values Container for the return values
* \param params Parameters
* \param state The fluid state
*/
template <class ContainerT, class FluidState>
static void capillaryPressures(ContainerT &values,
const Params &params,
const FluidState &state)
{
for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
Scalar S = state.saturation(phaseIdx);
Valgrind::CheckDefined(S);
values[phaseIdx] =
S*params.pcMaxSat(phaseIdx) +
(1.0 - S)*params.pcMinSat(phaseIdx);
}
}
/*!
* \brief The inverse of the capillary pressure
*/
template <class ContainerT, class FluidState>
static void saturations(ContainerT &values,
const Params &params,
const FluidState &state)
{
OPM_THROW(std::runtime_error, "Not implemented: MpLinearMaterial::saturations()");
}
/*!
* \brief The relative permeability of all phases.
*/
template <class ContainerT, class FluidState>
static void relativePermeabilities(ContainerT &values,
const Params &params,
const FluidState &state)
{
for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
Scalar S = state.saturation(phaseIdx);
Valgrind::CheckDefined(S);
values[phaseIdx] = std::max(std::min(S,1.0),0.0);
}
}
/*!
* \brief The difference between the pressures of the non-wetting and wetting phase.
*/
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{
Scalar S = fs.saturation(Traits::wPhaseIdx);
Valgrind::CheckDefined(S);
Scalar wPhasePressure =
S*params.pcMaxSat(Traits::wPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::wPhaseIdx);
S = fs.saturation(Traits::nPhaseIdx);
Valgrind::CheckDefined(S);
Scalar nPhasePressure =
S*params.pcMaxSat(Traits::nPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nPhaseIdx);
return nPhasePressure - wPhasePressure;
}
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatPcwn(const Params &params, Scalar Sw)
{
Scalar wPhasePressure =
Sw*params.pcMaxSat(Traits::wPhaseIdx) +
(1.0 - Sw)*params.pcMinSat(Traits::wPhaseIdx);
Scalar nPhasePressure =
(1.0 - Sw)*params.pcMaxSat(Traits::nPhaseIdx) +
Sw*params.pcMinSat(Traits::nPhaseIdx);
return nPhasePressure - wPhasePressure;
}
/*!
* \brief Calculate wetting phase saturation given that the rest
* of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{ OPM_THROW(std::runtime_error, "Not implemented: Sw()"); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatSw(const Params &params, Scalar Sw)
{ OPM_THROW(std::runtime_error, "Not implemented: twoPhaseSatSw()"); }
/*!
* \brief Calculate non-wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ OPM_THROW(std::runtime_error, "Not implemented: Sn()"); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatSn(const Params &params, Scalar Sw)
{ OPM_THROW(std::runtime_error, "Not implemented: twoPhaseSatSn()"); }
/*!
* \brief Calculate gas phase saturation given that the rest of
* the fluid state has been initialized
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
Sg(const Params &params, const FluidState &fs)
{ OPM_THROW(std::runtime_error, "Not implemented: Sg()"); }
/*!
* \brief The relative permability of the wetting phase
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wPhaseIdx))); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, Sw)); }
/*!
* \brief The relative permability of the liquid non-wetting phase
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nPhaseIdx))); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, 1 - Sw)); }
/*!
* \brief The relative permability of the gas phase
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
krg(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gPhaseIdx))); }
/*!
* \brief The difference between the pressures of the gas and the non-wetting phase.
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
pcng(const Params &params, const FluidState &fs)
{
Scalar S = fs.saturation(Traits::nPhaseIdx);
Valgrind::CheckDefined(S);
Scalar nPhasePressure =
S*params.pcMaxSat(Traits::nPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nPhaseIdx);
S = fs.saturation(Traits::gPhaseIdx);
Valgrind::CheckDefined(S);
Scalar gPhasePressure =
S*params.pcMaxSat(Traits::gPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::gPhaseIdx);
return gPhasePressure - nPhasePressure;
}
};
} // namespace Opm
#endif

57
opm/material/fluidmatrixinteractions/MpLinearMaterialParams.hpp → opm/material/fluidmatrixinteractions/LinearMaterialParams.hpp
View File

@@ -18,10 +18,10 @@
*****************************************************************************/
/*!
* \file
* \copydoc Opm::MpLinearMaterialParams
* \copydoc Opm::LinearMaterialParams
*/
#ifndef MP_LINEAR_MATERIAL_PARAMS_HH
#define MP_LINEAR_MATERIAL_PARAMS_HH
#ifndef OPM_LINEAR_MATERIAL_PARAMS_HH
#define OPM_LINEAR_MATERIAL_PARAMS_HH
namespace Opm {
@@ -29,66 +29,69 @@ namespace Opm {
* \brief Reference implementation of params for the linear M-phase
* material material.
*/
template<int numPhasesV, class ScalarT>
class MpLinearMaterialParams
template<class TraitsT>
class LinearMaterialParams
{
public:
typedef ScalarT Scalar;
enum { numPhases = numPhasesV };
enum { numPhases = TraitsT::numPhases };
typedef typename TraitsT::Scalar Scalar;
public:
typedef TraitsT Traits;
/*!
* \brief The default constructor.
*
* We set the capillary pressure to zero, if not specified otherwise.
*/
MpLinearMaterialParams()
LinearMaterialParams()
{
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
setPcMinSat(phaseIdx, 0.0);
setPcMaxSat(phaseIdx, 0.0);
setResidSat(phaseIdx, 0.0);
}
}
/*!
* \brief Return the capillary pressure for a phase \f$\alpha\f$ at \f$S_\alpha=0\f$.
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{ }
/*!
* \brief Return the relative phase pressure at the minimum saturation of a phase.
*
* This means \f$p_{c\alpha}\f$ at \f$S_\alpha=0\f$.
*/
Scalar pcMinSat(int phaseIdx) const
{ return pcMinSat_[phaseIdx]; }
/*!
* \brief Set the capillary pressure for a phase \f$\alpha\f$ at \f$S_\alpha=0\f$.
* \brief Set the relative phase pressure at the minimum saturation of a phase.
*
* This means \f$p_{c\alpha}\f$ at \f$S_\alpha=0\f$.
*/
void setPcMinSat(int phaseIdx, Scalar val)
{ pcMinSat_[phaseIdx] = val; }
/*!
* \brief Return the capillary pressure for a phase \f$\alpha\f$ at \f$S_\alpha=1\f$.
* \brief Return the relative phase pressure at the maximum saturation of a phase.
*
* This means \f$p_{c\alpha}\f$ at \f$S_\alpha=1\f$.
*/
Scalar pcMaxSat(int phaseIdx) const
{ return pcMaxSat_[phaseIdx]; }
/*!
* \brief Set the capillary pressure for a phase \f$\alpha\f$ at \f$S_\alpha=1\f$.
* \brief Set the relative phase pressure at the maximum saturation of a phase.
*
* This means \f$p_{c\alpha}\f$ at \f$S_\alpha=1\f$.
*/
void setPcMaxSat(int phaseIdx, Scalar val)
{ pcMaxSat_[phaseIdx] = val; }
/*!
* \brief Return the residual saturation for a phase \f$\alpha\f$.
*/
Scalar residSat(int phaseIdx) const
{ return residSat_[phaseIdx]; }
/*!
* \brief Set the residual saturation for a phase \f$\alpha\f$.
*/
void setResidSat(int phaseIdx, Scalar val)
{ residSat_[phaseIdx] = val; }
private:
Scalar residSat_[numPhases];
Scalar pcMaxSat_[numPhases];
Scalar pcMinSat_[numPhases];
};

103
opm/material/fluidmatrixinteractions/MaterialTraits.hpp Normal file
View File

@@ -0,0 +1,103 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2013 by Andreas Lauser *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief This file contains helper classes for the material laws.
*
* These classes specify the information which connects fluid systems
* and the fluid-matrix interaction laws. This includes stuff like the
* index of the wetting and non-wetting phase, etc.
*/
#ifndef OPM_MATERIAL_TRAITS_HH
#define OPM_MATERIAL_TRAITS_HH
namespace Opm {
/*!
* \ingroup material
*
* \brief A generic traits class for two-phase material laws.
*/
template <class ScalarT, int wettingPhaseIdxV, int nonWettingPhaseIdxV>
class TwoPhaseMaterialTraits
{
public:
//! The type used for scalar floating point values
typedef ScalarT Scalar;
//! The number of fluid phases
static const int numPhases = 2;
//! The index of the wetting phase
static const int wPhaseIdx = wettingPhaseIdxV;
//! The index of the non-wetting phase
static const int nPhaseIdx = nonWettingPhaseIdxV;
// some safety checks...
static_assert(0 <= wPhaseIdx && wPhaseIdx < numPhases,
"wPhaseIdx is out of range");
static_assert(0 <= nPhaseIdx && nPhaseIdx < numPhases,
"nPhaseIdx is out of range");
static_assert(wPhaseIdx != nPhaseIdx,
"wPhaseIdx and nPhaseIdx must be different");
};
/*!
* \ingroup material
*
* \brief A generic traits class for three-phase material laws.
*/
template <class ScalarT, int wettingPhaseIdxV, int nonWettingPhaseIdxV, int gasPhaseIdxV>
class ThreePhaseMaterialTraits
{
public:
//! The type used for scalar floating point values
typedef ScalarT Scalar;
//! The number of fluid phases
static const int numPhases = 3;
//! The index of the wetting liquid phase
static const int wPhaseIdx = wettingPhaseIdxV;
//! The index of the non-wetting liquid phase
static const int nPhaseIdx = nonWettingPhaseIdxV;
//! The index of the gas phase (i.e., the least wetting phase)
static const int gPhaseIdx = gasPhaseIdxV;
// some safety checks...
static_assert(0 <= wPhaseIdx && wPhaseIdx < numPhases,
"wPhaseIdx is out of range");
static_assert(0 <= nPhaseIdx && nPhaseIdx < numPhases,
"nPhaseIdx is out of range");
static_assert(0 <= gPhaseIdx && gPhaseIdx < numPhases,
"gPhaseIdx is out of range");
static_assert(wPhaseIdx != nPhaseIdx,
"wPhaseIdx and nPhaseIdx must be different");
static_assert(wPhaseIdx != gPhaseIdx,
"wPhaseIdx and gPhaseIdx must be different");
static_assert(nPhaseIdx != gPhaseIdx,
"nPhaseIdx and gPhaseIdx must be different");
};
} // namespace Opm
#endif

123
opm/material/fluidmatrixinteractions/MpLinearMaterial.hpp
View File

@@ -1,123 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2011-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::MpLinearMaterial
*/
#ifndef OPM_MP_LINEAR_MATERIAL_HH
#define OPM_MP_LINEAR_MATERIAL_HH
#include "MpLinearMaterialParams.hpp"
#include <opm/material/Valgrind.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
namespace Opm {
/*!
* \ingroup material
*
* \brief Implements a linear saturation-capillary pressure relation
*
* Implements a linear saturation-capillary pressure relation for
* M-phase fluid systems.
*
* \sa MpLinearMaterialParams
*/
template <int numPhasesV, class ScalarT, class ParamsT = MpLinearMaterialParams<numPhasesV, ScalarT> >
class MpLinearMaterial
{
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
enum { numPhases = numPhasesV };
/*!
* \brief The linear capillary pressure-saturation curve.
*
* This material law is linear:
* \f[
p_C = (1 - \overline{S}_w) (p_{C,max} - p_{C,entry}) + p_{C,entry}
\f]
*
* \param values Container for the return values
* \param params Parameters
* \param state The fluid state
*/
template <class ContainerT, class FluidState>
static void capillaryPressures(ContainerT &values,
const Params &params,
const FluidState &state)
{
Scalar sumResidSat = 0;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
sumResidSat += params.residSat(phaseIdx);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar Sabs = state.saturation(phaseIdx);
Valgrind::CheckDefined(Sabs);
Scalar S =
(Sabs - params.residSat(phaseIdx))
/ (1.0 - sumResidSat + params.residSat(phaseIdx));
values[phaseIdx] =
S*params.pcMaxSat(phaseIdx) +
(1.0 - S)*params.pcMinSat(phaseIdx);
}
}
/*!
* \brief The inverse of the capillary pressure
*/
template <class ContainerT, class FluidState>
static void saturations(ContainerT &values,
const Params &params,
const FluidState &state)
{
OPM_THROW(std::runtime_error, "Not implemented: MpLinearMaterial::saturations()");
}
/*!
* \brief The relative permeability of all phases.
*/
template <class ContainerT, class FluidState>
static void relativePermeabilities(ContainerT &values,
const Params &params,
const FluidState &state)
{
Scalar sumResidSat = 0;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
sumResidSat += params.residSat(phaseIdx);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar Sabs = state.saturation(phaseIdx);
Scalar S = (Sabs - params.residSat(phaseIdx))
/ (1.0 - sumResidSat + params.residSat(phaseIdx));
values[phaseIdx] = std::max(std::min(S,1.0),0.0);
}
}
};
} // namespace Opm
#endif

156
opm/material/fluidmatrixinteractions/2p/ParkerLenhard.hpp → opm/material/fluidmatrixinteractions/ParkerLenhard.hpp
View File

@@ -25,7 +25,7 @@
#include "ParkerLenhardParams.hpp"
#include <opm/material/fluidmatrixinteractions/2p/VanGenuchten.hpp>
#include <opm/material/fluidmatrixinteractions/VanGenuchten.hpp>
#include <algorithm>
#include <cassert>
@@ -58,13 +58,13 @@ public:
this, // next
-1, // loop number
Swr, // Sw
1e12, // pC
1e12, // pcwn
Swr, // SwMic
Swr); // SwMdc
next_ = NULL;
Sw_ = 1.0;
pC_ = 0.0;
pcwn_ = 0.0;
SwMic_ = 1.0;
SwMdc_ = 1.0;
}
@@ -74,7 +74,7 @@ protected:
PLScanningCurve *next,
int loopN,
Scalar Sw,
Scalar pC,
Scalar pcwn,
Scalar SwMic,
Scalar SwMdc)
{
@@ -82,7 +82,7 @@ protected:
next_ = next;
loopNum_ = loopN;
Sw_ = Sw;
pC_ = pC;
pcwn_ = pcwn;
SwMic_ = SwMic;
SwMdc_ = SwMdc;
}
@@ -124,7 +124,7 @@ public:
* deleted and thus forgotten.
*/
void setNext(Scalar Sw,
Scalar pC,
Scalar pcwn,
Scalar SwMic,
Scalar SwMdc)
{
@@ -136,7 +136,7 @@ public:
NULL, // next
loopNum() + 1,
Sw,
pC,
pcwn,
SwMic,
SwMdc);
}
@@ -195,8 +195,8 @@ public:
/*!
* \brief Capillary pressure at the last reversal point.
*/
Scalar pC() const
{ return pC_; }
Scalar pcwn() const
{ return pcwn_; }
/*!
* \brief Apparent saturation of the last reversal point on
@@ -219,7 +219,7 @@ private:
int loopNum_;
Scalar Sw_;
Scalar pC_;
Scalar pcwn_;
Scalar SwMdc_;
Scalar SwMic_;
@@ -231,12 +231,47 @@ private:
* p_c-Sw hysteresis model. This class adheres to the twophase
* capillary pressure API.
*/
template <class ScalarT, class ParamsT = ParkerLenhardParams<ScalarT> >
class ParkerLenhard
template <class TraitsT, class ParamsT = ParkerLenhardParams<TraitsT> >
class ParkerLenhard : public TraitsT
{
public:
typedef TraitsT Traits;
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The Parker-Lenhard capillary pressure law only "
"applies to the case of two fluid phases");
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
static_assert(Traits::numPhases == 2,
"The number of fluid phases must be two if you want to use "
"this material law!");
private:
typedef typename ParamsT::VanGenuchten VanGenuchten;
@@ -260,8 +295,11 @@ public:
* \brief Set the current absolute saturation for the
* current timestep
*/
static void update(Params &params, Scalar Sw)
template <class FluidState>
static void update(Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
if (Sw > 1 - 1e-5) {
// if the absolute saturation is almost 1,
// it means that we're back to the beginning
@@ -276,9 +314,9 @@ public:
// calculate the apparent saturation on the MIC and MDC
// which yield the same capillary pressure as the
// Sw at the current scanning curve
Scalar pc = pC(params, Sw);
Scalar Sw_mic = VanGenuchten::Sw(params.micParams(), pc);
Scalar Sw_mdc = VanGenuchten::Sw(params.mdcParams(), pc);
Scalar pc = pcwn(params, fs);
Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
curve->setNext(Sw, pc, Sw_mic, Sw_mdc);
if (!curve->next())
@@ -294,10 +332,44 @@ public:
}
/*!
* \brief Returns the capillary pressure dependend on
* the absolute__ saturation of the wetting phase.
* \brief Returns the capillary pressure dependening on
* the phase saturations.
*/
static Scalar pC(const Params &params, Scalar Sw)
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcwn(params, fs);
}
/*!
* \brief Returns the capillary pressure dependening on
* the phase saturations.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::saturations()"); }
/*!
* \brief Returns the relative permeabilities of the phases
* dependening on the phase saturations.
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
}
/*!
* \brief Returns the capillary pressure dependend on
* the phase saturations.
*/
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcwn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatPcwn(const Params &params, Scalar Sw)
{
// calculate the current apparent saturation
ScanningCurve *sc = findScanningCurve_(params, Sw);
@@ -308,7 +380,7 @@ public:
// if the apparent saturation exceeds the 'legal' limits,
// we also the underlying material law decide what to do.
if (Sw_app > 1) {
return 0.0; // VanGenuchten::pC(params.mdcParams(), Sw_app);
return 0.0; // VanGenuchten::pcwn(params.mdcParams(), Sw_app);
}
// put the apparent saturation into the main imbibition or
@@ -320,38 +392,68 @@ public:
Scalar SwMic =
pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
return VanGenuchten::pC(params.micParams(), SwMic);
return VanGenuchten::twoPhaseSatPcwn(params.micParams(), SwMic);
}
else { // sc->isDrain()
Scalar SwMdc =
pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
return VanGenuchten::pC(params.mdcParams(), SwMdc);
return VanGenuchten::twoPhaseSatPcwn(params.mdcParams(), SwMdc);
}
}
/*!
* \brief Calculate the wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::Sw()"); }
static Scalar twoPhaseSatSw(const Params &params, Scalar pc)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pc)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
/*!
* \brief The relative permeability for the wetting phase of
* the medium.
*/
static Scalar krw(const Params &params, Scalar Sw)
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
// for the effective permeability we only use Land's law and
// the relative permeability of the main drainage curve.
Scalar Sw_app = absoluteToApparentSw_(params, Sw);
return VanGenuchten::krw(params.mdcParams(), Sw_app);
return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
}
/*!
* \brief The relative permeability for the non-wetting phase
* of the params.
*/
static Scalar krn(const Params &params, Scalar Sw)
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
// for the effective permeability we only use Land's law and
// the relative permeability of the main drainage curve.
Scalar Sw_app = absoluteToApparentSw_(params, Sw);
return VanGenuchten::krn(params.mdcParams(), Sw_app);
return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
}
/*!

15
opm/material/fluidmatrixinteractions/2p/ParkerLenhardParams.hpp → opm/material/fluidmatrixinteractions/ParkerLenhardParams.hpp
View File

@@ -23,7 +23,7 @@
#ifndef OPM_PARKER_LENHARD_PARAMS_HH
#define OPM_PARKER_LENHARD_PARAMS_HH
#include <opm/material/fluidmatrixinteractions/2p/RegularizedVanGenuchten.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
namespace Opm
{
@@ -35,12 +35,12 @@ class PLScanningCurve;
* \brief Default parameter class for the Parker-Lenhard hysteresis
* model.
*/
template <class ScalarT>
template <class TraitsT>
class ParkerLenhardParams
{
public:
typedef ScalarT Scalar;
typedef Opm::RegularizedVanGenuchten<Scalar> VanGenuchten;
typedef typename TraitsT::Scalar Scalar;
typedef Opm::RegularizedVanGenuchten<TraitsT> VanGenuchten;
typedef typename VanGenuchten::Params VanGenuchtenParams;
typedef PLScanningCurve<Scalar> ScanningCurve;
@@ -62,6 +62,13 @@ public:
~ParkerLenhardParams()
{ delete mdc_; }
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{ }
/*!
* \brief Returns the parameters of the main imbibition curve (which uses
* the van Genuchten capillary pressure model).

196
opm/material/fluidmatrixinteractions/2p/RegularizedBrooksCorey.hpp → opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp
View File

@@ -58,14 +58,94 @@ namespace Opm {
*
* \see BrooksCorey
*/
template <class ScalarT, class ParamsT = RegularizedBrooksCoreyParams<ScalarT> >
class RegularizedBrooksCorey
template <class TraitsT, class ParamsT = RegularizedBrooksCoreyParams<TraitsT> >
class RegularizedBrooksCorey : public TraitsT
{
typedef Opm::BrooksCorey<ScalarT, ParamsT> BrooksCorey;
typedef Opm::BrooksCorey<TraitsT, ParamsT> BrooksCorey;
public:
typedef TraitsT Traits;
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The regularized Brooks-Corey capillary pressure law only "
"applies to the case of two fluid phases");
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
static_assert(Traits::numPhases == 2,
"The number of fluid phases must be two if you want to use "
"this material law!");
/*!
* \brief The capillary pressure-saturation curves depending on absolute saturations.
*
* \param values A random access container which stores the
* relative pressure of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the capillary pressure
* ought to be calculated
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcwn(params, fs);
}
/*!
* \brief Calculate the saturations of the phases starting from
* their pressure differences.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
}
/*!
* \brief The relative permeability-saturation curves depending on absolute saturations.
*
* \param values A random access container which stores the
* relative permeability of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the relative permeabilities
* ought to be calculated
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
}
/*!
* \brief A regularized Brooks-Corey capillary pressure-saturation
@@ -89,9 +169,13 @@ public:
* slope of the unregularized Brooks-Corey curve at \f$S_w =
* 1\f$
*
* \sa BrooksCorey::pC
* \sa BrooksCorey::pcwn
*/
static Scalar pC(const Params &params, Scalar Sw)
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcwn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatPcwn(const Params &params, Scalar Sw)
{
const Scalar Sthres = params.thresholdSw();
@@ -102,28 +186,35 @@ public:
// saturation moving to the right direction if it
// temporarily is in an 'illegal' range.
if (Sw <= Sthres) {
Scalar m = BrooksCorey::dpC_dSw(params, Sthres);
Scalar pC_SwLow = BrooksCorey::pC(params, Sthres);
return pC_SwLow + m*(Sw - Sthres);
Scalar m = BrooksCorey::twoPhaseSatDpcwn_dSw(params, Sthres);
Scalar pcwn_SwLow = BrooksCorey::twoPhaseSatPcwn(params, Sthres);
return pcwn_SwLow + m*(Sw - Sthres);
}
else if (Sw > 1.0) {
Scalar m = BrooksCorey::dpC_dSw(params, 1.0);
Scalar pC_SwHigh = BrooksCorey::pC(params, 1.0);
return pC_SwHigh + m*(Sw - 1.0);
Scalar m = BrooksCorey::twoPhaseSatDpcwn_dSw(params, 1.0);
Scalar pcwn_SwHigh = params.entryPressure();
return pcwn_SwHigh + m*(Sw - 1.0);
}
// if the effective saturation is in an 'reasonable'
// range, we use the real Brooks-Corey law...
return BrooksCorey::pC(params, Sw);
return BrooksCorey::twoPhaseSatPcwn(params, Sw);
}
/*!
* \brief A regularized Brooks-Corey saturation-capillary pressure
* curve.
*
* This is the inverse of the pC() method.
* This is the inverse of the pcwn() method.
*/
static Scalar Sw(const Params &params, Scalar pC)
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pcwn = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
return twoPhaseSatSw(params, pcwn);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pcwn)
{
const Scalar Sthres = params.thresholdSw();
@@ -133,8 +224,8 @@ public:
// smaller than the entry pressure, make sure that we will
// regularize.
Scalar Sw = 1.5;
if (pC >= params.pe())
Sw = BrooksCorey::Sw(params, pC);
if (pcwn >= params.entryPressure())
Sw = BrooksCorey::twoPhaseSatSw(params, pcwn);
// make sure that the capilary pressure observes a
// derivative != 0 for 'illegal' saturations. This is
@@ -143,48 +234,59 @@ public:
// saturation moving to the right direction if it
// temporarily is in an 'illegal' range.
if (Sw <= Sthres) {
// invert the low saturation regularization of pC()
Scalar m = BrooksCorey::dpC_dSw(params, Sthres);
Scalar pC_SwLow = BrooksCorey::pC(params, Sthres);
return Sthres + (pC - pC_SwLow)/m;
// invert the low saturation regularization of pcwn()
Scalar m = BrooksCorey::twoPhaseSatDpcwn_dSw(params, Sthres);
Scalar pcwn_SwLow = BrooksCorey::twoPhaseSatPcwn(params, Sthres);
return Sthres + (pcwn - pcwn_SwLow)/m;
}
else if (Sw > 1.0) {
Scalar m = BrooksCorey::dpC_dSw(params, 1.0);
Scalar pC_SwHigh = BrooksCorey::pC(params, 1.0);
return 1.0 + (pC - pC_SwHigh)/m;;
Scalar m = BrooksCorey::twoPhaseSatDpcwn_dSw(params, 1.0);
Scalar pcwn_SwHigh = BrooksCorey::twoPhaseSatPcwn(params, 1.0);
return 1.0 + (pcwn - pcwn_SwHigh)/m;;
}
return BrooksCorey::Sw(params, pC);
return BrooksCorey::twoPhaseSatSw(params, pcwn);
}
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pcwn)
{ return 1 - twoPhaseSatSw(params, pcwn); }
/*!
* \brief The derivative of the regularized Brooks-Corey capillary
* pressure-saturation curve.
*/
static Scalar dpC_dSw(const Params &params, Scalar Sw)
static Scalar dpcwn_dSw(const Params &params, Scalar Sw)
{
const Scalar Sthres = params.thresholdSw();
// derivative of the regualarization
if (Sw <= Sthres) {
// calculate the slope of the straight line used in pC()
Scalar m = BrooksCorey::dpC_dSw(params, Sthres);
// calculate the slope of the straight line used in pcwn()
Scalar m = BrooksCorey::dpcwn_dSw(params, Sthres);
return m;
}
else if (Sw > 1.0) {
// calculate the slope of the straight line used in pC()
Scalar m = BrooksCorey::dpC_dSw(params, 1.0);
// calculate the slope of the straight line used in pcwn()
Scalar m = BrooksCorey::dpcwn_dSw(params, 1.0);
return m;
}
return BrooksCorey::dpC_dSw(params, Sw);
return BrooksCorey::dpcwn_dSw(params, Sw);
}
/*!
* \brief The derivative of the regularized Brooks-Corey
* saturation-capillary pressure curve.
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
static Scalar dSw_dpcwn(const Params &params, Scalar pcwn)
{
const Scalar Sthres = params.thresholdSw();
@@ -192,23 +294,23 @@ public:
// saturation in the non-regularized version of the
// Brooks-Corey law
Scalar Sw;
if (pC < params.pe())
if (pcwn < params.entryPressure())
Sw = 1.5; // make sure we regularize (see below)
else
Sw = BrooksCorey::Sw(params, pC);
Sw = BrooksCorey::Sw(params, pcwn);
// derivative of the regularization
if (Sw <= Sthres) {
// calculate the slope of the straight line used in pC()
Scalar m = BrooksCorey::dpC_dSw(params, Sthres);
// calculate the slope of the straight line used in pcwn()
Scalar m = BrooksCorey::dpcwn_dSw(params, Sthres);
return 1/m;
}
else if (Sw > 1.0) {
// calculate the slope of the straight line used in pC()
Scalar m = BrooksCorey::dpC_dSw(params, 1.0);
// calculate the slope of the straight line used in pcwn()
Scalar m = BrooksCorey::dpcwn_dSw(params, 1.0);
return 1/m;
}
return 1.0/BrooksCorey::dpC_dSw(params, Sw);
return 1.0/BrooksCorey::dpcwn_dSw(params, Sw);
}
/*!
@@ -225,14 +327,18 @@ public:
*
* \sa BrooksCorey::krw
*/
static Scalar krw(const Params &params, Scalar Sw)
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
if (Sw <= 0.0)
return 0.0;
else if (Sw >= 1.0)
return 1.0;
return BrooksCorey::krw(params, Sw);
return BrooksCorey::twoPhaseSatKrw(params, Sw);
}
/*!
@@ -249,14 +355,18 @@ public:
*
* \sa BrooksCorey::krn
*/
static Scalar krn(const Params &params, Scalar Sw)
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
if (Sw >= 1.0)
return 0.0;
else if (Sw <= 0.0)
return 1.0;
return BrooksCorey::krn(params, Sw);
return BrooksCorey::twoPhaseSatKrn(params, Sw);
}
};
} // namespace Opm

22
opm/material/fluidmatrixinteractions/2p/RegularizedBrooksCoreyParams.hpp → opm/material/fluidmatrixinteractions/RegularizedBrooksCoreyParams.hpp
View File

@@ -32,24 +32,34 @@ namespace Opm {
* \brief Parameters that are necessary for the \em regularization of
* the Brooks-Corey capillary pressure model.
*/
template <class ScalarT>
class RegularizedBrooksCoreyParams : public Opm::BrooksCoreyParams<ScalarT>
template <class TraitsT>
class RegularizedBrooksCoreyParams : public Opm::BrooksCoreyParams<TraitsT>
{
typedef Opm::BrooksCoreyParams<ScalarT> BrooksCoreyParams;
typedef Opm::BrooksCoreyParams<TraitsT> BrooksCoreyParams;
typedef typename TraitsT::Scalar Scalar;
public:
typedef ScalarT Scalar;
typedef TraitsT Traits;
RegularizedBrooksCoreyParams()
: BrooksCoreyParams()
, SwThres_(1e-2)
{ }
RegularizedBrooksCoreyParams(Scalar pe, Scalar lambda)
: BrooksCoreyParams(pe, lambda)
RegularizedBrooksCoreyParams(Scalar entryPressure, Scalar lambda)
: BrooksCoreyParams(entryPressure, lambda)
, SwThres_(1e-2)
{ }
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{
BrooksCoreyParams::finalize();
}
/*!
* \brief Return the threshold saturation below which the capillary pressure
* is regularized.

316
opm/material/fluidmatrixinteractions/2p/RegularizedVanGenuchten.hpp → opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp
View File

@@ -65,14 +65,79 @@ namespace Opm {
*
* \see VanGenuchten
*/
template <class ScalarT, class ParamsT = RegularizedVanGenuchtenParams<ScalarT> >
class RegularizedVanGenuchten
template <class TraitsT, class ParamsT = RegularizedVanGenuchtenParams<TraitsT> >
class RegularizedVanGenuchten : public TraitsT
{
typedef Opm::VanGenuchten<ScalarT, ParamsT> VanGenuchten;
typedef Opm::VanGenuchten<TraitsT, ParamsT> VanGenuchten;
public:
typedef TraitsT Traits;
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The regularized van Genuchten capillary pressure law only "
"applies to the case of two fluid phases");
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
/*!
* \brief Calculate the pressure difference of the phases in the
* most generic way.
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcwn(params, fs);
}
/*!
* \brief Calculate the saturations of the phases starting from
* their pressure differences.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
}
/*!
* \brief Returns the relative permeabilities of the phases
* dependening on the phase saturations.
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
}
/*!
* \brief A regularized van Genuchten capillary pressure-saturation
@@ -84,45 +149,47 @@ public:
*
* For not-regularized part:
*
\copydetails VanGenuchten::pC()
* \copydetails VanGenuchten::pC()
*/
static Scalar pC(const Params &params, Scalar Swe)
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcwn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatPcwn(const Params &params, Scalar Sw)
{
// retrieve the low and the high threshold saturations for the
// unregularized capillary pressure curve from the parameters
const Scalar SwThLow = params.pCLowSw();
const Scalar SwThHigh = params.pCHighSw();
const Scalar SwThLow = params.pcwnLowSw();
const Scalar SwThHigh = params.pcwnHighSw();
// make sure that the capillary pressure observes a derivative
// != 0 for 'illegal' saturations. This is favourable for the
// newton solver (if the derivative is calculated numerically)
// in order to get the saturation moving to the right
// direction if it temporarily is in an 'illegal' range.
if (Swe < SwThLow) {
return VanGenuchten::pC(params, SwThLow) + mLow_(params)*(Swe - SwThLow);
if (Sw < SwThLow) {
return params.pcwnLow() + params.pcwnSlopeLow()*(Sw - SwThLow);
}
else if (Swe > SwThHigh)
else if (Sw > SwThHigh)
{
Scalar yTh = VanGenuchten::pC(params, SwThHigh);
Scalar yTh = params.pcwnHigh();
Scalar m1 = (0.0 - yTh)/(1.0 - SwThHigh)*2;
if (Swe < 1.0) {
// use spline between threshold Swe and 1.0
Scalar mTh = VanGenuchten::dpC_dSw(params, SwThHigh);
Spline<Scalar> sp(SwThHigh, 1.0, // x0, x1
yTh, 0, // y0, y1
mTh, m1); // m0, m1
return sp.eval(Swe);
if (Sw < 1.0) {
// use spline between threshold Sw and 1.0
const Spline<Scalar> &sp = params.pcwnHighSpline();
return sp.eval(Sw);
}
else {
// straight line for Swe > 1.0
return m1*(Swe - 1.0) + 0.0;
// straight line for Sw > 1.0
return m1*(Sw - 1.0) + 0.0;
}
}
// if the effective saturation is in an 'reasonable'
// range, we use the real van genuchten law...
return VanGenuchten::pC(params, Swe);
return VanGenuchten::twoPhaseSatPcwn(params, Sw);
}
/*!
@@ -139,76 +206,96 @@ public:
\copydetails VanGenuchten::Sw()
*
*/
static Scalar Sw(const Params &params, Scalar pC)
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
return twoPhaseSatSw(params, pC);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
{
// retrieve the low and the high threshold saturations for the
// unregularized capillary pressure curve from the parameters
const Scalar SwThLow = params.pCLowSw();
const Scalar SwThHigh = params.pCHighSw();
const Scalar SwThLow = params.pcwnLowSw();
const Scalar SwThHigh = params.pcwnHighSw();
// calculate the saturation which corrosponds to the
// saturation in the non-regularized verision of van
// Genuchten's law
Scalar Sw;
if (pC <= 0) {
// invert straight line for Swe > 1.0
Scalar yTh = VanGenuchten::pC(params, SwThHigh);
Scalar m1 = (0.0 - yTh)/(1.0 - SwThHigh)*2;
// invert straight line for Sw > 1.0
Scalar m1 = params.pcwnSlopeHigh();
return pC/m1 + 1.0;
}
else
Sw = VanGenuchten::Sw(params, pC);
Sw = VanGenuchten::twoPhaseSatSw(params, pC);
// invert the regularization if necessary
if (Sw <= SwThLow) {
// invert the low saturation regularization of pC()
Scalar pC_SwLow = VanGenuchten::pC(params, SwThLow);
return (pC - pC_SwLow)/mLow_(params) + SwThLow;
Scalar pC_SwLow = VanGenuchten::twoPhaseSatPcwn(params, SwThLow);
return (pC - pC_SwLow)/params.pcwnSlopeLow() + SwThLow;
}
else if (Sw > SwThHigh)
else if (SwThHigh < Sw /* && Sw < 1.0*/)
{
Scalar yTh = VanGenuchten::pC(params, SwThHigh);
Scalar m1 = (0.0 - yTh)/(1.0 - SwThHigh)*2;
// invert spline between threshold saturation and 1.0
const Spline<Scalar>& spline = params.pcwnHighSpline();
// invert spline between threshold Swe and 1.0
Scalar mTh = VanGenuchten::dpC_dSw(params, SwThHigh);
Spline<Scalar> sp(SwThHigh, 1.0, // x0, x1
yTh, 0, // m0, m1
mTh, m1); // m0, m1
return sp.intersectInterval(SwThHigh, 1.0,
0, 0, 0, pC);
return spline.intersectInterval(/*x0=*/SwThHigh, /*x1=*/1.0,
/*a=*/0, /*b=*/0, /*c=*/0, /*d=*/pC);
}
return Sw;
}
/*!
* \brief A regularized version of the partial derivative
* of the \f$p_c(\overline S_w)\f$ w.r.t. effective saturation
* according to van Genuchten.
*
* regularized part:
* - low saturation: use the slope of the regularization point (i.e. no kink).
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line and use that slope (yes, there is a kink :-( ).
*
* For not-regularized part:
*
\copydetails VanGenuchten::dpC_dSw()
*
*/
static Scalar dpC_dSw(const Params &params, Scalar Swe)
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pc)
{ return 1 - twoPhaseSatSw(params, pc); }
/*!
* \brief A regularized version of the partial derivative
* of the \f$p_c(\overline S_w)\f$ w.r.t. effective saturation
* according to van Genuchten.
*
* regularized part:
* - low saturation: use the slope of the regularization point (i.e. no kink).
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line and use that slope (yes, there is a kink :-( ).
*
* For not-regularized part:
*
* \copydetails VanGenuchten::dpC_dSw()
*
*/
template <class FluidState>
static Scalar dpcwn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDpcwn_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatDpcwn_dSw(const Params &params, Scalar Sw)
{
// derivative of the regualarization
if (Swe < params.pCLowSw()) {
if (Sw < params.pcwnLowSw()) {
// the slope of the straight line used in pC()
return mLow_(params);
return params.pcwnSlopeLow();
}
else if (Swe > params.pCHighSw()) {
// the slope of the straight line used in pC()
return mHigh_(params);
else if (params.pcwnHighSw() <= Sw) {
if (Sw < 1)
return params.pcwnHighSpline().evalDerivative(Sw);
else
// the slope of the straight line used for the right
// side of the capillary pressure function
return params.mHigh();
}
return VanGenuchten::dpC_dSw(params, Swe);
return VanGenuchten::dpcwn_dSw(params, Sw);
}
/*!
@@ -221,10 +308,13 @@ public:
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line and use that slope (yes, there is a kink :-( ).
*
* For not-regularized part:
\copydetails VanGenuchten::dSw_dpC()
* \copydetails VanGenuchten::dSw_dpC()
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
template <class FluidState>
static Scalar dSw_dpC(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
// calculate the saturation which corrosponds to the
// saturation in the non-regularized verision of van
// Genuchten's law
@@ -232,19 +322,22 @@ public:
if (pC < 0)
Sw = 1.5; // make sure we regularize below
else
Sw = VanGenuchten::Sw(params, pC);
Sw = VanGenuchten::Sw_raw(params, pC);
// derivative of the regularization
if (Sw < params.pCLowSw()) {
if (Sw < params.pcwnLowSw()) {
// same as in dpC_dSw() but inverted
return 1/mLow_(params);
return 1/params.pcwnSlopeLow();
}
if (Sw > params.pCHighSw()) {
if (Sw > params.pcwnHighSw()) {
if (Sw < 1)
return 1/params.pcwnHighSpline().evalDerivative(Sw);
// same as in dpC_dSw() but inverted
return 1/mHigh_(params);
return 1/params.pcwnSlopHigh();
}
return VanGenuchten::dSw_dpC(params, pC);
return VanGenuchten::dSw_dpnw(params, fs);
}
/*!
@@ -261,26 +354,19 @@ public:
* For not-regularized part:
\copydetails VanGenuchten::krw()
*/
static Scalar krw(const Params &params, Scalar Swe)
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
// retrieve the high threshold saturation for the
// unregularized relative permeability curve of the wetting
// phase from the parameters
const Scalar SwThHigh = params.krwHighSw();
if (Swe < 0)
// regularize
if (Sw < 0)
return 0;
else if (Swe > 1)
else if (Sw > 1)
return 1;
else if (Swe > SwThHigh) {
typedef Opm::Spline<Scalar> Spline;
Spline sp(SwThHigh, 1.0, // x1, x2
VanGenuchten::krw(params, SwThHigh), 1.0, // y1, y2
VanGenuchten::dkrw_dSw(params, SwThHigh), 0); // m1, m2
return sp.eval(Swe);
}
return VanGenuchten::krw(params, Swe);
return VanGenuchten::twoPhaseSatKrw(params, Sw);
}
/*!
@@ -297,61 +383,19 @@ public:
\copydetails VanGenuchten::krn()
*
*/
static Scalar krn(const Params &params, Scalar Swe)
{
// retrieve the low threshold saturation for the unregularized
// relative permeability curve of the non-wetting phase from
// the parameters
const Scalar SwThLow = params.krnLowSw();
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wPhaseIdx)); }
if (Swe <= 0)
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
// regularize
if (Sw <= 0)
return 1;
else if (Swe >= 1)
else if (Sw >= 1)
return 0;
else if (Swe < SwThLow) {
typedef Opm::Spline<Scalar> Spline;
Spline sp(0.0, SwThLow, // x1, x2
1.0, VanGenuchten::krn(params, SwThLow), // y1, y2
0.0, VanGenuchten::dkrn_dSw(params, SwThLow)); // m1, m2
return sp.eval(Swe);
}
return VanGenuchten::krn(params, Swe);
}
private:
// the slope of the straight line used to regularize saturations
// below the minimum saturation
/*!
* \brief The slope of the straight line used to regularize
* saturations below the minimum saturation.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar mLow_(const Params &params)
{
const Scalar SwThLow = params.pCLowSw();
return VanGenuchten::dpC_dSw(params, SwThLow);
}
/*!
* \brief The slope of the straight line used to regularize
* saturations above the minimum saturation.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar mHigh_(const Params &params)
{
const Scalar SwThHigh = params.pCHighSw();
Scalar pC_SwHigh = VanGenuchten::pC(params, SwThHigh);
return (0 - pC_SwHigh)/(1.0 - SwThHigh);
return VanGenuchten::twoPhaseSatKrn(params, Sw);
}
};

161
opm/material/fluidmatrixinteractions/RegularizedVanGenuchtenParams.hpp Normal file
View File

@@ -0,0 +1,161 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2008-2012 by Andreas Lauser *
* *
* 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
* \copydoc Opm::RegularizedVanGenuchtenParams
*/
#ifndef OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
#define OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
#include "VanGenuchtenParams.hpp"
#include <opm/core/utility/Spline.hpp>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Parameters that are necessary for the \em regularization of
* VanGenuchten "material law".
*
*/
template<class TraitsT>
class RegularizedVanGenuchtenParams : public VanGenuchtenParams<TraitsT>
{
typedef typename TraitsT::Scalar Scalar;
typedef VanGenuchtenParams<TraitsT> Parent;
typedef Opm::VanGenuchten<TraitsT> VanGenuchten;
public:
typedef TraitsT Traits;
RegularizedVanGenuchtenParams()
: pcwnLowSw_(0.01)
, pcwnHighSw_(0.99)
{}
RegularizedVanGenuchtenParams(Scalar vgAlpha, Scalar vgN)
: Parent(vgAlpha, vgN)
, pcwnLowSw_(0.01)
, pcwnHighSw_(0.99)
{}
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{
Parent::finalize();
pcwnLow_ = VanGenuchten::twoPhaseSatPcwn(*this, pcwnLowSw_);
pcwnSlopeLow_ = VanGenuchten::twoPhaseSatDpcwn_dSw(*this, pcwnLowSw_);
pcwnHigh_ = VanGenuchten::twoPhaseSatPcwn(*this, pcwnHighSw_);
pcwnSlopeHigh_ = 2*(0.0 - pcwnHigh_)/(1.0 - pcwnHighSw_);
Scalar mThreshold = VanGenuchten::twoPhaseSatDpcwn_dSw(*this, pcwnHighSw_);
pcwnHighSpline_.set(pcwnHighSw_, 1.0, // x0, x1
pcwnHigh_, 0, // y0, y1
mThreshold, pcwnSlopeHigh_); // m0, m1
}
/*!
* \brief Return the threshold saturation below which the
* capillary pressure is regularized.
*/
Scalar pcwnLowSw() const
{ return pcwnLowSw_; }
/*!
* \brief Return the capillary pressure at the low threshold
* saturation of the wetting phase.
*/
Scalar pcwnLow() const
{ return pcwnLow_; }
/*!
* \brief Return the slope capillary pressure curve if Sw is
* smaller or equal to the low threshold saturation.
*
* For this case, we extrapolate the curve using a straight line.
*/
Scalar pcwnSlopeLow() const
{ return pcwnSlopeLow_; }
/*!
* \brief Set the threshold saturation below which the capillary
* pressure is regularized.
*/
void setPCLowSw(Scalar value)
{ pcwnLowSw_ = value; }
/*!
* \brief Return the threshold saturation below which the
* capillary pressure is regularized.
*/
Scalar pcwnHighSw() const
{ return pcwnHighSw_; }
/*!
* \brief Return the capillary pressure at the high threshold
* saturation of the wetting phase.
*/
Scalar pcwnHigh() const
{ return pcwnHigh_; }
/*!
* \brief Return the spline curve which ought to be used between
* the upper threshold saturation and 1.
*/
const Spline<Scalar> &pcwnHighSpline() const
{ return pcwnHighSpline_; }
/*!
* \brief Return the slope capillary pressure curve if Sw is
* larger or equal to 1.
*
* For this case, we extrapolate the curve using a straight line.
*/
Scalar pcwnSlopeHigh() const
{ return pcwnSlopeHigh_; }
/*!
* \brief Set the threshold saturation below which the capillary
* pressure is regularized.
*/
void setPCHighSw(Scalar value)
{ pcwnHighSw_ = value; }
private:
Scalar pcwnLowSw_;
Scalar pcwnHighSw_;
Scalar pcwnLow_;
Scalar pcwnHigh_;
Scalar pcwnSlopeLow_;
Scalar pcwnSlopeHigh_;
Spline<Scalar> pcwnHighSpline_;
};
} // namespace Opm
#endif

130
opm/material/fluidmatrixinteractions/2p/VanGenuchten.hpp → opm/material/fluidmatrixinteractions/VanGenuchten.hpp
View File

@@ -47,27 +47,48 @@ namespace Opm {
*
* \see VanGenuchtenParams
*/
template <class ScalarT,
int wPhaseIdxV,
int nPhaseIdxV,
class ParamsT = VanGenuchtenParams<ScalarT> >
class VanGenuchten
template <class TraitsT, class ParamsT = VanGenuchtenParams<TraitsT> >
class VanGenuchten : public TraitsT
{
public:
//! The traits class for this material law
typedef TraitsT Traits;
//! The type of the parameter objects for this law
typedef ParamsT Params;
//! The type of the scalar values for this law
typedef typename Params::Scalar Scalar;
typedef typename Traits::Scalar Scalar;
//! The number of fluid phases to which this capillary pressure law applies
enum { numPhases = 2 };
//! The number of fluid phases
static const int numPhases = Traits::numPhases;
static_assert(numPhases == 2,
"The van Genuchten capillary pressure law only "
"applies to the case of two fluid phases");
//! The index of the wetting phase
enum { wPhaseIdx = wPhaseIdxV };
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
//! The index of the non-wetting phase
enum { nPhaseIdx = nPhaseIdxV };
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
/*!
* \brief The capillary pressure-saturation curves according to van Genuchten.
@@ -88,8 +109,19 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[wPhaseIdx] = 0.0; // reference phase
values[nPhaseIdx] = pcwn(params, fs);
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcwn(params, fs);
}
/*!
* \brief Calculate the saturations of the phases starting from
* their pressure differences.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
}
/*!
@@ -105,8 +137,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[wPhaseIdx] = krw(params, fs);
values[nPhaseIdx] = krn(params, fs);
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
}
/*!
@@ -126,11 +158,29 @@ public:
template <class FluidState>
static Scalar pcwn(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
assert(0 <= Sw && Sw <= 1);
return std::pow(std::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha();
return twoPhaseSatPcwn(params, Sw);
}
/*!
* \brief The saturation-capillary pressure curve according to van
* Genuchten using a material law specific API.
*
* The advantage of this model is that it is simpler to use
* because the baggage of the fluid state API does not need to be
* carried along. The disavantage of this is, that it is very
* specific to the van Genuchten law (i.e., depends only on the
* wetting phase saturation, assumes two fluid phases, etc)
*
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param Sw The effective wetting phase saturation
*/
static Scalar twoPhaseSatPcwn(const Params &params, Scalar Sw)
{ return std::pow(std::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha(); }
/*!
* \brief The saturation-capillary pressure curve according to van Genuchten.
*
@@ -146,12 +196,28 @@ public:
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(nPhaseIdx) - fs.pressure(wPhaseIdx);
Scalar pC = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
return twoPhaseSatSw(params, pC);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
{
assert(pC >= 0);
return std::pow(std::pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM());
}
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pC)
{ return 1 - twoPhaseSatSw(params, pC); }
/*!
* \brief The partial derivative of the capillary pressure with
* regard to the saturation according to van Genuchten.
@@ -168,9 +234,10 @@ public:
*/
template <class FluidState>
static Scalar dpcwn_dSw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(wPhaseIdx);
{ return twoPhaseSatDpcwn_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatDpcwn_dSw(const Params &params, Scalar Sw)
{
assert(0 < Sw && Sw < 1);
Scalar powSw = std::pow(Sw, -1/params.vgM());
@@ -191,9 +258,10 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(wPhaseIdx);
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
Scalar r = 1. - std::pow(1 - std::pow(Sw, 1/params.vgM()), params.vgM());
@@ -213,9 +281,10 @@ public:
*/
template <class FluidState>
static Scalar dkrw_dSw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(wPhaseIdx);
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatDkrw_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
const Scalar x = 1 - std::pow(Sw, 1.0/params.vgM());
@@ -233,10 +302,12 @@ public:
* \param fs The fluid state for which the derivative
* ought to be calculated
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(wPhaseIdx);
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return
@@ -257,9 +328,10 @@ public:
*/
template <class FluidState>
static Scalar dkrn_dSw(const Params &params, const FluidState &fs)
{
Scalar fs = fs.saturation(wPhaseIdx);
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
static Scalar twoPhaseSatDkrn_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
const Scalar x = std::pow(Sw, 1.0/params.vgM());

13
opm/material/fluidmatrixinteractions/2p/VanGenuchtenParams.hpp → opm/material/fluidmatrixinteractions/VanGenuchtenParams.hpp
View File

@@ -34,11 +34,13 @@ namespace Opm {
* parameter automatically calculates the other. I.e. they cannot be
* set independently.
*/
template<class ScalarT>
template<class TraitsT>
class VanGenuchtenParams
{
typedef typename TraitsT::Scalar Scalar;
public:
typedef ScalarT Scalar;
typedef TraitsT Traits;
VanGenuchtenParams()
{}
@@ -49,6 +51,13 @@ public:
setVgN(vgN);
}
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{ }
/*!
* \brief Return the \f$\alpha\f$ shape parameter of van Genuchten's
* curve.

15
opm/material/fluidstates/ImmiscibleFluidState.hpp
View File

@@ -29,8 +29,11 @@
#include "ModularFluidState.hpp"
#include <opm/material/Valgrind.hpp>
#include <algorithm>
#include <string.h>
namespace Opm {
/*!
@@ -55,6 +58,12 @@ class ImmiscibleFluidState<Scalar, FluidSystem, true>
FluidStateExplicitViscosityModule<Scalar, FluidSystem, ImmiscibleFluidState<Scalar, FluidSystem, true> >,
FluidStateExplicitEnthalpyModule<Scalar, FluidSystem, ImmiscibleFluidState<Scalar, FluidSystem, true> > >
{
public:
ImmiscibleFluidState()
{}
ImmiscibleFluidState(const ImmiscibleFluidState &fs)
{ memcpy(this, &fs, sizeof(fs)); }
};
// specialization for the enthalpy disabled case
@@ -71,6 +80,12 @@ class ImmiscibleFluidState<Scalar, FluidSystem, false>
FluidStateExplicitViscosityModule<Scalar, FluidSystem, ImmiscibleFluidState<Scalar, FluidSystem, false> >,
FluidStateNullEnthalpyModule<Scalar, FluidSystem, ImmiscibleFluidState<Scalar, FluidSystem, false> > >
{
public:
ImmiscibleFluidState()
{}
ImmiscibleFluidState(const ImmiscibleFluidState &fs)
{ memcpy(this, &fs, sizeof(fs)); }
};
} // namespace Opm

30
tests/material/immiscibleflash/test_immiscibleflash.cpp
View File

@@ -36,12 +36,10 @@
#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
#include <opm/material/fluidmatrixinteractions/MpLinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/2pAdapter.hpp>
#include <opm/material/fluidmatrixinteractions/2p/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/2p/RegularizedLinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/2p/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/2p/EffToAbsLaw.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
template <class Scalar, class FluidState>
void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
@@ -159,9 +157,9 @@ int main()
enum { H2OIdx = FluidSystem::H2OIdx };
enum { N2Idx = FluidSystem::N2Idx };
typedef Opm::RegularizedBrooksCorey<Scalar> EffMaterialLaw;
typedef Opm::EffToAbsLaw<EffMaterialLaw> TwoPMaterialLaw;
typedef Opm::TwoPAdapter<lPhaseIdx, TwoPMaterialLaw> MaterialLaw;
typedef Opm::TwoPhaseMaterialTraits<Scalar, lPhaseIdx, gPhaseIdx> MaterialLawTraits;
typedef Opm::RegularizedBrooksCorey<MaterialLawTraits> EffMaterialLaw;
typedef Opm::EffToAbsLaw<EffMaterialLaw> MaterialLaw;
typedef MaterialLaw::Params MaterialLawParams;
Scalar T = 273.15 + 25;
@@ -179,10 +177,11 @@ int main()
// set the parameters for the capillary pressure law
MaterialLawParams matParams;
matParams.setSwr(0.0);
matParams.setSnr(0.0);
matParams.setPe(0);
matParams.setResidualSaturation(MaterialLaw::wPhaseIdx, 0.0);
matParams.setResidualSaturation(MaterialLaw::nPhaseIdx, 0.0);
matParams.setEntryPressure(0);
matParams.setLambda(2.0);
matParams.finalize();
ImmiscibleFluidState fsRef;
@@ -242,10 +241,11 @@ int main()
std::cout << "testing two-phase with capillary pressure\n";
MaterialLawParams matParams2;
matParams2.setSwr(0.0);
matParams2.setSnr(0.0);
matParams2.setPe(1e3);
matParams2.setResidualSaturation(MaterialLaw::wPhaseIdx, 0.0);
matParams2.setResidualSaturation(MaterialLaw::nPhaseIdx, 0.0);
matParams2.setEntryPressure(1e3);
matParams2.setLambda(2.0);
matParams2.finalize();
// set liquid saturation
fsRef.setSaturation(lPhaseIdx, 0.5);

31
tests/material/ncpflash/test_ncpflash.cpp
View File

@@ -36,12 +36,10 @@
#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
#include <opm/material/fluidmatrixinteractions/MpLinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/2pAdapter.hpp>
#include <opm/material/fluidmatrixinteractions/2p/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/2p/RegularizedLinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/2p/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/2p/EffToAbsLaw.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
template <class Scalar, class FluidState>
void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
@@ -162,9 +160,9 @@ int main()
enum { H2OIdx = FluidSystem::H2OIdx };
enum { N2Idx = FluidSystem::N2Idx };
typedef Opm::RegularizedBrooksCorey<Scalar> EffMaterialLaw;
typedef Opm::EffToAbsLaw<EffMaterialLaw> TwoPMaterialLaw;
typedef Opm::TwoPAdapter<lPhaseIdx, TwoPMaterialLaw> MaterialLaw;
typedef Opm::TwoPhaseMaterialTraits<Scalar, lPhaseIdx, gPhaseIdx> MaterialTraits;
typedef Opm::RegularizedBrooksCorey<MaterialTraits> EffMaterialLaw;
typedef Opm::EffToAbsLaw<EffMaterialLaw> MaterialLaw;
typedef MaterialLaw::Params MaterialLawParams;
Scalar T = 273.15 + 25;
@@ -182,10 +180,11 @@ int main()
// set the parameters for the capillary pressure law
MaterialLawParams matParams;
matParams.setSwr(0.0);
matParams.setSnr(0.0);
matParams.setPe(0);
matParams.setResidualSaturation(MaterialLaw::wPhaseIdx, 0.0);
matParams.setResidualSaturation(MaterialLaw::nPhaseIdx, 0.0);
matParams.setEntryPressure(0);
matParams.setLambda(2.0);
matParams.finalize();
CompositionalFluidState fsRef;
@@ -260,12 +259,12 @@ int main()
////////////////
// with capillary pressure
////////////////
MaterialLawParams matParams2;
matParams2.setSwr(0.0);
matParams2.setSnr(0.0);
matParams2.setPe(1e3);
matParams2.setResidualSaturation(MaterialLaw::wPhaseIdx, 0.0);
matParams2.setResidualSaturation(MaterialLaw::nPhaseIdx, 0.0);
matParams2.setEntryPressure(1e3);
matParams2.setLambda(2.0);
matParams2.finalize();
// set gas saturation
fsRef.setSaturation(gPhaseIdx, 0.5);

11
tests/material/pengrobinson/test_pengrobinson.cpp
View File

@@ -28,7 +28,8 @@
#include <opm/material/constraintsolvers/NcpFlash.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/fluidsystems/Spe5FluidSystem.hpp>
#include <opm/material/fluidmatrixinteractions/MpLinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
template <class FluidSystem, class FluidState>
void guessInitial(FluidState &fluidState,
@@ -70,7 +71,8 @@ Scalar bringOilToSurface(FluidState &surfaceFluidState, Scalar alpha, const Flui
};
typedef Opm::NcpFlash<Scalar, FluidSystem> Flash;
typedef Opm::MpLinearMaterial<numPhases, Scalar> MaterialLaw;
typedef Opm::ThreePhaseMaterialTraits<Scalar, wPhaseIdx, oPhaseIdx, gPhaseIdx> MaterialTraits;
typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
@@ -82,6 +84,7 @@ Scalar bringOilToSurface(FluidState &surfaceFluidState, Scalar alpha, const Flui
matParams.setPcMinSat(phaseIdx, 0.0);
matParams.setPcMaxSat(phaseIdx, 0.0);
}
matParams.finalize();
// retieve the global volumetric component molarities
surfaceFluidState.setTemperature(273.15 + 20);
@@ -171,7 +174,8 @@ int main(int argc, char** argv)
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Opm::CompositionalFluidState<Scalar, FluidSystem> FluidState;
typedef Opm::MpLinearMaterial<numPhases, Scalar> MaterialLaw;
typedef Opm::ThreePhaseMaterialTraits<Scalar, wPhaseIdx, oPhaseIdx, gPhaseIdx> MaterialTraits;
typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
typedef MaterialLaw::Params MaterialLawParams;
typedef FluidSystem::ParameterCache ParameterCache;
@@ -192,6 +196,7 @@ int main(int argc, char** argv)
matParams.setPcMinSat(phaseIdx, 0.0);
matParams.setPcMaxSat(phaseIdx, 0.0);
}
matParams.finalize();
////////////
// Create a fluid state

301
tests/test_fluidmatrixinteractions.cpp Normal file
View File

@@ -0,0 +1,301 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2013 by Andreas Lauser *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief This test makes sure that the API for fluid-matrix
* interactions is observed by all capillary pressure / relperm
* laws.
*/
#include "config.h"
// include all capillary pressure laws
#include <opm/material/fluidmatrixinteractions/BrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/ParkerLenhard.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/VanGenuchten.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
// include the helper classes to construct traits
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
// include some fluid states
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
// include some fluid systems
#include <opm/material/fluidsystems/2pImmiscibleFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
// include some components
#include <opm/material/components/SimpleH2O.hpp>
#include <opm/material/components/N2.hpp>
#include <opm/core/utility/Unused.hpp>
// include the MPI header if available
#if HAVE_MPI
#include <mpi.h>
#endif // HAVE_MPI
// class to call MPI_Init() on construction and MPI_Finalize() in the
// destructor
class MyMpiHelper
{
public:
MyMpiHelper(int &argc, char **&argv)
{
#if HAVE_MPI
MPI_Init(&argc, &argv);
#endif // HAVE_MPI
};
~MyMpiHelper()
{
#if HAVE_MPI
MPI_Finalize();
#endif // HAVE_MPI
};
};
// this function makes sure that a capillary pressure law adheres to
// the generic programming interface for such laws. This API _must_ be
// implemented by all capillary pressure laws. If there are no _very_
// strong reasons to do otherwise, numerical models should only use on
// this API.
template <class MaterialLaw, class FluidState>
void testGenericApi()
{
while (0) {
// ensure the presence of the required values
static const int numPhases = MaterialLaw::numPhases;
// check for the presence of the is*Dependent values
OPM_UNUSED static const bool isSaturationDependent = MaterialLaw::isSaturationDependent;
OPM_UNUSED static const bool isPressureDependent = MaterialLaw::isPressureDependent;
OPM_UNUSED static const bool isTemperatureDependent = MaterialLaw::isTemperatureDependent;
OPM_UNUSED static const bool isCompositionDependent = MaterialLaw::isCompositionDependent;
// Make sure that the Traits, Params and Scalar typedefs are
// exported by the material law
typedef typename MaterialLaw::Params Params;
typedef typename MaterialLaw::Traits Traits;
typedef typename MaterialLaw::Scalar Scalar;
typedef typename MaterialLaw::Traits::Scalar TraitsScalar;
static_assert(std::is_same<Scalar, TraitsScalar>::value,
"The traits and the material law must use the same type as scalar value");
static_assert(numPhases == Traits::numPhases,
"The traits and the material law must use the number of fluid phases");
// check the API of the parameter class. setting the actual
// parameter values is implementation specific. But all
// parameters must be default and copy constructible as well
// as exhibit the finalize() method!
Params params;
params.finalize();
const Params paramsConst(params);
// test the generic methods which need to be implemented by
// all material laws
const FluidState fs;
double destValues[numPhases];
MaterialLaw::capillaryPressures(destValues, paramsConst, fs);
MaterialLaw::saturations(destValues, paramsConst, fs);
MaterialLaw::relativePermeabilities(destValues, paramsConst, fs);
}
}
// this function makes ensures that a pressure law adheres to the
// covenience programming interface for two-phase material laws. The
// main purpose of this interface is to simplify the implementation of
// nested material laws.
template <class MaterialLaw, class FluidState>
void testTwoPhaseApi()
{
typedef typename MaterialLaw::Scalar Scalar;
while (0) {
static const int numPhases = MaterialLaw::numPhases;
static_assert(numPhases == 2,
"The number of fluid phases for a twophase "
"capillary pressure law must be 2");
static_assert(MaterialLaw::implementsTwoPhaseApi,
"This material law is expected to implement "
"the two-phase API!");
OPM_UNUSED static const int wPhaseIdx = MaterialLaw::wPhaseIdx;
OPM_UNUSED static const int nPhaseIdx = MaterialLaw::nPhaseIdx;
// make sure the two-phase specific methods are present
const FluidState fs;
const typename MaterialLaw::Params params;
Scalar v;
v = MaterialLaw::pcwn(params, fs);
v = MaterialLaw::Sw(params, fs);
v = MaterialLaw::Sn(params, fs);
v = MaterialLaw::krw(params, fs);
v = MaterialLaw::krn(params, fs);
}
}
template <class MaterialLaw, class FluidState>
void testTwoPhaseSatApi()
{
typedef typename MaterialLaw::Scalar Scalar;
while (0) {
static_assert(MaterialLaw::implementsTwoPhaseSatApi,
"This material law is expected to implement "
"the two-phase saturation only API!");
OPM_UNUSED static const int numPhases = MaterialLaw::numPhases;
// make sure the two-phase specific methods are present
const typename MaterialLaw::Params params;
Scalar Sw;
Scalar v;
v = MaterialLaw::twoPhaseSatPcwn(params, Sw);
v = MaterialLaw::twoPhaseSatSw(params, Sw);
v = MaterialLaw::twoPhaseSatSn(params, Sw);
v = MaterialLaw::twoPhaseSatKrw(params, Sw);
v = MaterialLaw::twoPhaseSatKrn(params, Sw);
}
}
template <class MaterialLaw, class FluidState>
void testThreePhaseApi()
{
typedef typename MaterialLaw::Scalar Scalar;
while (0) {
static const int numPhases = MaterialLaw::numPhases;
static_assert(numPhases == 3,
"The number of fluid phases for a threephase "
"capillary pressure law must be 3");
OPM_UNUSED static const int wPhaseIdx = MaterialLaw::wPhaseIdx;
OPM_UNUSED static const int nPhaseIdx = MaterialLaw::nPhaseIdx;
OPM_UNUSED static const int gPhaseIdx = MaterialLaw::gPhaseIdx;
// make sure the two-phase specific methods are present
const FluidState fs;
const typename MaterialLaw::Params params;
Scalar v;
v = MaterialLaw::pcwn(params, fs);
v = MaterialLaw::Sw(params, fs);
v = MaterialLaw::Sn(params, fs);
v = MaterialLaw::krw(params, fs);
v = MaterialLaw::krn(params, fs);
v = MaterialLaw::krg(params, fs);
}
}
template <class MaterialLaw>
void testThreePhaseSatApi()
{
}
int main(int argc, char **argv)
{
typedef double Scalar;
typedef Opm::SimpleH2O<Scalar> H2O;
typedef Opm::N2<Scalar> N2;
typedef Opm::LiquidPhase<Scalar, H2O> Liquid;
typedef Opm::GasPhase<Scalar, N2> Gas;
typedef Opm::FluidSystems::TwoPImmiscible<Scalar, Liquid, Gas> TwoPFluidSystem;
typedef Opm::FluidSystems::BlackOil<Scalar> ThreePFluidSystem;
typedef Opm::TwoPhaseMaterialTraits<Scalar,
TwoPFluidSystem::wPhaseIdx,
TwoPFluidSystem::nPhaseIdx> TwoPhaseTraits;
typedef Opm::ThreePhaseMaterialTraits<Scalar,
ThreePFluidSystem::wPhaseIdx,
ThreePFluidSystem::oPhaseIdx,
ThreePFluidSystem::gPhaseIdx> ThreePhaseTraits;
typedef Opm::ImmiscibleFluidState<Scalar, TwoPFluidSystem> TwoPhaseFluidState;
typedef Opm::ImmiscibleFluidState<Scalar, ThreePFluidSystem> ThreePhaseFluidState;
MyMpiHelper mpiHelper(argc, argv);
// test conformance to the capillary pressure APIs
{
typedef Opm::BrooksCorey<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::LinearMaterial<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
typedef Opm::EffToAbsLaw<MaterialLaw> TwoPAbsLaw;
testGenericApi<TwoPAbsLaw, TwoPhaseFluidState>();
testTwoPhaseApi<TwoPAbsLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<TwoPAbsLaw, TwoPhaseFluidState>();
typedef Opm::LinearMaterial<ThreePhaseTraits> ThreePMaterialLaw;
testGenericApi<ThreePMaterialLaw, ThreePhaseFluidState>();
testThreePhaseApi<ThreePMaterialLaw, ThreePhaseFluidState>();
//testThreePhaseSatApi<ThreePMaterialLaw, ThreePhaseFluidState>();
typedef Opm::EffToAbsLaw<ThreePMaterialLaw> ThreePAbsLaw;
testGenericApi<ThreePAbsLaw, ThreePhaseFluidState>();
testThreePhaseApi<ThreePAbsLaw, ThreePhaseFluidState>();
//testThreePhaseSatApi<ThreePAbsLaw, ThreePhaseFluidState>();
}
{
typedef Opm::ParkerLenhard<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::VanGenuchten<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::RegularizedBrooksCorey<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::RegularizedVanGenuchten<TwoPhaseTraits> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
return 0;
}