Skip to content
Snippets Groups Projects

Adapt to new material laws

Merged Ned Coltman requested to merge feature/newMaterialLaws into fix/new-matlaw-warnings
All threads resolved!
Files
36
@@ -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_;
};
Loading