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

[mm][mcwhorter] adapt to new material laws

parent d587cb1e
......@@ -60,8 +60,6 @@ class McWhorterAnalytic
using ElementIterator = typename GridView::template Codim<0>::Iterator;
public:
using MaterialLaw = typename SpatialParams::MaterialLaw;
using MaterialLawParams = typename MaterialLaw::Params;
// functions needed for analytical solution
void initializeAnalytic()
......@@ -97,11 +95,12 @@ public:
void prepareAnalytic()
{
const MaterialLawParams& materialLawParams( problem_.spatialParams().materialLawParams(dummyElement_) );
swr_ = materialLawParams.swr();
snr_ = materialLawParams.snr();
porosity_ = problem_.spatialParams().porosity(dummyElement_);
permeability_ = problem_.spatialParams().intrinsicPermeability(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);
permeability_ = problem_.spatialParams().intrinsicPermeability(dummyElement);
PrimaryVariables initVec;
problem_.initialAtPos(initVec, dummyGlobal_);
sInit_ = initVec[saturationIdx];
......@@ -126,8 +125,8 @@ public:
// get fractional flow function vector
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);
}
// get capillary pressure derivatives
......@@ -135,7 +134,7 @@ public:
for (int i=0; i<pointNum_; i++)
{
dpcdsw_[i] = MaterialLaw::dpc_dsw(materialLawParams, satVec_[i]);
dpcdsw_[i] = fluidMatrixInteraction.dpc_dsw(satVec_[i]);
}
// set initial fW
......@@ -156,7 +155,7 @@ public:
// diffusivity function
for (int i=0; i<pointNum_; i++)
{
d_[i] = fractionalW_[i]*(-dpcdsw_[i])*(MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW)*permeability_;
d_[i] = fractionalW_[i]*(-dpcdsw_[i])*(fluidMatrixInteraction.krn(satVec_[i])/viscosityNW)*permeability_;
}
// gk_: part of fractional flow function
......@@ -302,9 +301,13 @@ public:
prepareAnalytic();
}
McWhorterAnalytic(Problem& problem, Scalar tolAnalytic = 1e-14) :
problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), dummyElement_(
*(problem_.gridView().template begin<0> ())), tolAnalytic_(tolAnalytic)
McWhorterAnalytic(Problem& problem, Scalar tolAnalytic = 1e-14)
: problem_(problem)
, analyticSolution_(0)
, error_(0)
, elementVolume_(0)
, size_(problem.gridView().size(0))
, tolAnalytic_(tolAnalytic)
{
dummyGlobal_ = 0.0;
dummyGlobal_[0] = 1.0;
......@@ -318,8 +321,6 @@ private:
BlockVector elementVolume_;
int size_;
const Element& dummyElement_;
GlobalPosition dummyGlobal_;
Scalar tolAnalytic_;
......
......@@ -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"
......
......@@ -9,9 +9,9 @@ EnableGravity = false
Permeability = 1.01936799e-13 # intrinsic permeability of the porous medium [m^2]
Porosity = 0.2 # porosity of the porous medium [-]
BrooksCoreyLambda = 3.0 # pore size distribution parameter for the Brooks-Corey capillary pressure - saturation relationship [-]
BrooksCoreyEntryPressure = 8000 # 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 = 8000 # 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]
......
......@@ -79,8 +79,8 @@ public:
PseudoOil<Scalar>::setDensity( getParam<Scalar>("Fluid.DensityW") );
PseudoH2O<Scalar>::setDensity( 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);
}
......
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