Commit d587cb1e authored by Ned Coltman's avatar Ned Coltman
Browse files

[mm][buckleyLeverett] adapt to new material laws

parent 2d76ea3e
......@@ -21,8 +21,7 @@
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/2p/sequential/properties.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh>
/*!
* \file
......@@ -36,33 +35,6 @@ namespace Dumux {
* the Buckley-Leverett problem
*/
template<typename Scalar, typename Law>
struct CheckMaterialLaw
{
static bool isLinear()
{
return false;
}
};
template<typename Scalar>
struct CheckMaterialLaw<Scalar, LinearMaterial<Scalar> >
{
static bool isLinear()
{
return true;
}
};
template<typename Scalar>
struct CheckMaterialLaw<Scalar, EffToAbsLaw< LinearMaterial<Scalar> > >
{
static bool isLinear()
{
return true;
}
};
template<class TypeTag> class BuckleyLeverettAnalytic
{
using Problem = GetPropType<TypeTag, Properties::Problem>;
......@@ -70,8 +42,6 @@ template<class TypeTag> class BuckleyLeverettAnalytic
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
using MaterialLaw = typename SpatialParams::MaterialLaw;
using MaterialLawParams = typename MaterialLaw::Params;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using FluidState = GetPropType<TypeTag, Properties::FluidState>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
......@@ -87,7 +57,6 @@ template<class TypeTag> class BuckleyLeverettAnalytic
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
public:
// functions needed for analytical solution
void initializeAnalytic()
......@@ -106,11 +75,11 @@ public:
*/
void prepareAnalytic()
{
const MaterialLawParams& materialLawParams(problem_.spatialParams().materialLawParams(dummyElement_));
swr_ = materialLawParams.swr();
snr_ = materialLawParams.snr();
porosity_ = problem_.spatialParams().porosity(dummyElement_);
const auto& dummyElement = *problem_.gridView().template begin<0>();
const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(dummyElement.geometry().center());
swr_ = fluidMatrixInteraction.effToAbsParams().swr();
snr_ = fluidMatrixInteraction.effToAbsParams().snr();
porosity_ = problem_.spatialParams().porosity(dummyElement);
time_ = 0;
satVec_ = swr_;
......@@ -128,12 +97,11 @@ public:
for (int i = 0; i < pointNum_; i++)
{
fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW;
fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW);
fractionalW_[i] = fluidMatrixInteraction.krw(satVec_[i])/viscosityW;
fractionalW_[i] /= (fractionalW_[i] + fluidMatrixInteraction.krn(satVec_[i])/viscosityNW);
}
dfwdsw_ = 0;
for (int i = 1; i < intervalNum_; i++)
{
dfwdsw_[i] = (fractionalW_[i + 1] - fractionalW_[i - 1]) / (satVec_[i + 1] - satVec_[i - 1]);
......@@ -177,12 +145,13 @@ public:
{
xf_[i] = vTot_ * time_ / porosity_ * dfwdsw_[i];
}
// position of maximum xf_
int xfmax = 0;
int xhelp = pointNum_ / 3;
int xhelpold = 0;
int xhelpoldold = 0;
int xfmax = 0;
// position of maximum xf_
int xhelp2 = 0;
for (int i = 0; i < pointNum_; i++)
{
if (xf_[i] > xf_[i + 1])
......@@ -197,13 +166,9 @@ public:
Scalar A1;
Scalar A2;
Scalar b;
int xhelp2 = 0;
while (a)
{
if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
break;
A1 = 0;
for (int i = 0; i <= xhelp - 1; i++)
......@@ -262,7 +227,7 @@ public:
int index = problem_.variables().index(*eIt);
// account for linear material law
if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
if constexpr (SpatialParams::pcSwCurveIsLinear())
{
if (globalPos[0] <= xf_[1])
{
......@@ -313,9 +278,7 @@ public:
}
BlockVector AnalyticSolution() const
{
return analyticSolution_;
}
{ return analyticSolution_; }
//Write saturation and pressure into file
template<class MultiWriter>
......@@ -332,9 +295,13 @@ public:
}
//! Construct an IMPES object.
BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7) :
problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), vTot_(totalVelocity), dummyElement_(
*(problem_.gridView().template begin<0> ()))
BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7)
: problem_(problem)
, analyticSolution_(0)
, error_(0)
, elementVolume_(0)
, size_(problem.gridView().size(0))
, vTot_(totalVelocity)
{
dummyGlobal_ = 0.0;
dummyGlobal_[0] = 1.0;
......@@ -365,7 +332,6 @@ private:
Dune::FieldVector<Scalar, pointNum_> dfwdsw_;
Dune::FieldVector<Scalar, pointNum_> xf_;
int dfwdswmax_;
const Element& dummyElement_;
GlobalPosition dummyGlobal_;
};
......
......@@ -48,10 +48,10 @@ void usage(const char *progName, const std::string &errorMsg)
"\t-SpatialParams.Porosity The porosity of the porous medium [-]\n"
"\t-SpatialParams.BrooksCoreyLambda The pore size distribution parameter for the \n"
"\t \t Brooks-Corey capillary pressure - saturation relationship [-]\n"
"\t-SpatialParams.BrooksCoreyEntryPressure The entry pressure for the \n"
"\t-SpatialParams.BrooksCoreyPcEntry The entry pressure for the \n"
"\t \t Brooks-Corey capillary pressure - saturation relationship [Pa]\n"
"\t-SpatialParams.ResidualSaturationWetting The residual saturation of the wetting phase [-]\n"
"\t-SpatialParams.ResidualSaturationNonwetting The residual saturation of the nonwetting phase [-]\n"
"\t-SpatialParams.Swr The residual saturation of the wetting phase [-]\n"
"\t-SpatialParams.Snr The residual saturation of the nonwetting phase [-]\n"
"\t-Fluid.DensityW The density of the wetting phase [kg/m^3]\n"
"\t-Fluid.DensityNW The density of the nonwetting phase [kg/m^3]\n"
"\t-Fluid.ViscosityW The dynamic viscosity of the wetting phase [kg/(ms)]\n"
......
......@@ -6,15 +6,12 @@ DtInitial = 1e3 # initial time step for the simulati
EnableGravity = false
[SpatialParams]
Permeability = 1.01936799e-14 # intrinsic permeability of the porous medium [m^2]
Porosity = 0.2 # porosity of the porous medium [-]
BrooksCoreyLambda = 4.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship [-]
BrooksCoreyEntryPressure = 0 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship [Pa]
ResidualSaturationWetting = 0.2 # residual saturation of the wetting phase [-]
ResidualSaturationNonwetting = 0.2 # residual saturation of the nonwetting phase [-]
BrooksCoreyPcEntry = 0 # entry pressure for the Brooks-Corey capillary pressure - saturation relationship [Pa]
Swr = 0.2 # residual saturation of the wetting phase [-]
Snr = 0.2 # residual saturation of the nonwetting phase [-]
[Fluid]
DensityW = 1e3 # density of the wetting phase [kg/m^3]
......
......@@ -90,8 +90,8 @@ public:
densityNonwetting_ = getParam<Scalar>("Fluid.DensityNW");
swr_ = getParam<Scalar>("SpatialParams.ResidualSaturationWetting");
snr_ = getParam<Scalar>("SpatialParams.ResidualSaturationNonwetting");
swr_ = getParam<Scalar>("SpatialParams.Swr");
snr_ = getParam<Scalar>("SpatialParams.Snr");
paraviewOutput_ = getParam<bool>("Output.paraviewOutput", true);
}
......
......@@ -23,9 +23,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/material/spatialparams/sequentialfv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
namespace Dumux {
......@@ -44,19 +42,12 @@ struct BuckleyLeverettSpatialParamsTypeTag {};
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::BuckleyLeverettSpatialParamsTypeTag> { using type = BuckleyLeverettSpatialParams<TypeTag>; };
// Set the material law
template<class TypeTag>
struct MaterialLaw<TypeTag, TTag::BuckleyLeverettSpatialParamsTypeTag>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using RawMaterialLaw = RegularizedBrooksCorey<Scalar>;
public:
using type = EffToAbsLaw<RawMaterialLaw>;
};
} // end namespace Properties
// forward declaration
class LinearMaterialDefault;
class LinearMaterial;
template<class TypeTag>
class BuckleyLeverettSpatialParams: public SequentialFVSpatialParams<TypeTag>
{
......@@ -71,50 +62,44 @@ class BuckleyLeverettSpatialParams: public SequentialFVSpatialParams<TypeTag>
using Element = typename Grid::Traits::template Codim<0>::Entity;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
using FieldMatrix = Dune::FieldMatrix<Scalar,dim,dim>;
using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>;
public:
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
Scalar intrinsicPermeability (const Element& element) const
BuckleyLeverettSpatialParams(const Problem& problem)
: ParentType(problem)
, pcKrSwCurve_("SpatialParams")
{
return constPermeability_;
Scalar permFactor = 1.0; //0.001/(1000*9.81);
constPermeability_ = getParam<Scalar>("SpatialParams.Permeability")*permFactor;
porosity_ = getParam<Scalar>("SpatialParams.Porosity");
}
static constexpr bool pcSwCurveIsLinear()
{ return (std::is_same_v<PcKrSwCurve, LinearMaterial> || std::is_same_v<PcKrSwCurve, LinearMaterialDefault>); }
Scalar intrinsicPermeability (const Element& element) const
{ return constPermeability_; }
Scalar porosity(const Element &element) const
{
return porosity_;
}
{ return porosity_; }
/*!
* \brief DOC ME!
* \param element DOC ME!
* \brief Returns the parameters for the material law at a given location
*
* \param globalPos The global coordinates for the given location
*/
// return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const Element &element) const
{
return materialLawParams_;
}
BuckleyLeverettSpatialParams(const Problem& problem)
:ParentType(problem)
{
Scalar permFactor = 1.0; //0.001/(1000*9.81);
constPermeability_ = getParam<double>("SpatialParams.Permeability")*permFactor;
materialLawParams_.setSwr( getParam<double>("SpatialParams.ResidualSaturationWetting") );
materialLawParams_.setSnr( getParam<double>("SpatialParams.ResidualSaturationNonwetting") );
//set Brooks-Corey parameters
materialLawParams_.setPe( getParam<double>("SpatialParams.BrooksCoreyEntryPressure") );
materialLawParams_.setLambda( getParam<double>("SpatialParams.BrooksCoreyLambda") );
porosity_ = getParam<double>("SpatialParams.Porosity");
}
auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{ return pcKrSwCurve_; }
private:
MaterialLawParams materialLawParams_;
const PcKrSwCurve pcKrSwCurve_;
Scalar constPermeability_;
Scalar porosity_;
Scalar swr_;
Scalar snr_;
Scalar pcEntry_;
Scalar lambda_;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment