Commit 49cbd33e authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[multidomain,Beavers-Joseph] Move evalBeaversJoseph to LOP which prevents a...

[multidomain,Beavers-Joseph] Move evalBeaversJoseph to LOP which prevents a reimplementation in the non-isothermal coupling local residual.

Move beaversJospehCoeff to SpatialParamters, where it belongs.
Changes in reference solution was triggered, comments see below.

2cstokes2p2c:
Checked storage.out (e.g. the evaporation rate):
  - changes were lass than 1 per mile.
Checked vertical velocity profiles in the middle of the domain:
  - v_0: no visible difference
  - v_1: main deviation (5%) in middle, where the profile is oscillating anyway and NOT at the interface
Checked horizontal velocity profiles:
  - v_0: small differences towards the outflow boundary
  - v_1: difference during the profile development

2cnizeroeq2p2cni (freeflow):
Checked vertical velocity profiles in the middle of the domain:
  - v_0: no visible difference
  - v_1: no visible difference
Checked horizontal velocity profiles:
  - v_0: no visible difference
  - v_1: small differences towards the outflow boundary

2cnizeroeq2p2cni (porousmedium):
Checked vertical velocity profiles in the middle of the domain:
  - v_0: no visible difference
  - v_1: very small deviation towards the interface (3%)
parent 2691a7a7
......@@ -81,7 +81,6 @@ protected:
};
typedef typename GridView::ctype CoordScalar;
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> DimVector;
......@@ -179,13 +178,6 @@ public:
averagedNormal += boundaryFaceNormal;
}
// TODO: move scope below to coupling localoperator
// Beavers-Joseph condition at the coupling boundary/interface
if(bcTypes.hasCoupling())
{
evalBeaversJoseph_(&intersection, scvIdx, boundaryFaceIdx, boundaryVars);
}
// TODO: move scope below to coupling localoperator/ BUG (potentially): sollte das nicht dirichlet sein?
// set velocity normal to the interface
if (bcTypes.isCouplingNeumann(momentumYIdx))
......@@ -330,6 +322,7 @@ protected:
}
template <class IntersectionIterator>
DUNE_DEPRECATED_MSG("evalBeaversJoseph_ is deprecated. Its functionality is now included in the LOP function evalCoupling21().")
void evalBeaversJoseph_(const IntersectionIterator &isIt,
const int scvIdx,
const int boundaryFaceIdx,
......
......@@ -139,7 +139,8 @@ class TwoCStokesTwoPTwoCLocalOperator :
typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) Stokes2cTypeTag;
typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) TwoPTwoCTypeTag;
typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(Stokes2cTypeTag, SpatialParams) SpatialParams;
typedef typename GET_PROP_TYPE(Stokes2cTypeTag, ElementVolumeVariables) ElementVolumeVariables1;
typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, ElementVolumeVariables) ElementVolumeVariables2;
......@@ -165,7 +166,10 @@ class TwoCStokesTwoPTwoCLocalOperator :
typedef typename GET_PROP_TYPE(Stokes2cTypeTag, Indices) Stokes2cIndices;
typedef typename GET_PROP_TYPE(TwoPTwoCTypeTag, Indices) TwoPTwoCIndices;
enum { dim = MDGrid::dimension };
enum {
dim = MDGrid::dimension,
dimWorld = MDGrid::dimensionworld
};
// FREE FLOW
enum { numEq1 = GET_PROP_VALUE(Stokes2cTypeTag, NumEq) };
......@@ -199,6 +203,9 @@ class TwoCStokesTwoPTwoCLocalOperator :
nPhaseIdx2 = TwoPTwoCIndices::nPhaseIdx //!< Index for the gas phase
};
typedef typename MDGrid::ctype CoordScalar;
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dim> DimVector;
......@@ -231,6 +238,7 @@ class TwoCStokesTwoPTwoCLocalOperator :
* Based on them, a coupling residual is calculated and added at the
* respective positions in the matrix.
*
* TODO: rename _s _n parameters
* \param intersectionGeometry the geometry of the intersection
* \param lfsu_s local basis for the trial space of the Stokes domain
* \param unknowns1 the unknowns vector of the Stokes element (formatted according to PDELab)
......@@ -616,6 +624,54 @@ class TwoCStokesTwoPTwoCLocalOperator :
}
}
SpatialParams spatialParams = globalProblem_.sdProblem1().spatialParams();
const GlobalPosition& globalPos = cParams.fvGeometry1.subContVol[vertInElem1].global;
Scalar beaversJosephCoeff = spatialParams.beaversJosephCoeffAtPos(globalPos);
assert(beaversJosephCoeff > 0);
const Scalar Kxx = spatialParams.intrinsicPermeability(sdElement1, cParams.fvGeometry1,
vertInElem1);
beaversJosephCoeff /= std::sqrt(Kxx);
const DimVector& elementUnitNormal = boundaryVars1.face().normal;
// TODO revise comment
// Implementation as Neumann condition: (v.n)n
if (cParams.boundaryTypes1.isCouplingDirichlet(momentumXIdx1))
{
const Scalar normalComp = boundaryVars1.velocity()*elementUnitNormal;
DimVector normalV = elementUnitNormal;
normalV *= normalComp; // v*n*n
// v_tau = v - (v.n)n
const DimVector tangentialV = boundaryVars1.velocity() - normalV;
const Scalar boundaryFaceArea = boundaryVars1.face().area;
for (int dimIdx=0; dimIdx < dim; ++dimIdx)
{
couplingRes1.accumulate(lfsu_s.child(momentumXIdx1), vertInElem1,
beaversJosephCoeff
* boundaryFaceArea
* tangentialV[dimIdx]
* (boundaryVars1.dynamicViscosity()
+ boundaryVars1.dynamicEddyViscosity()));
}
}
// TODO revise comment
// Implementation as Dirichlet condition
// tangential component: vx = sqrt K /alpha * (grad v n(unity))t
if (cParams.boundaryTypes1.isCouplingNeumann(momentumXIdx1))
{
DimVector tangentialVelGrad(0);
boundaryVars1.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
tangentialVelGrad /= -beaversJosephCoeff; // was - before
// TODO: 3 lines below not implemtented because curPrimaryVars_ is protected
// it could be that this part of code has never been checked
// couplingRes1.accumulate(lfsu_s.child(momentumXIdx1), vertInElem1,
// this->residual_[vertInElem1][momentumXIdx1] =
// tangentialVelGrad[momentumXIdx1] - globalProblem_.localResidual1().curPriVars_(vertInElem1)[momentumXIdx1]);
}
//coupling residual is added to "real" residual
//here each node is passed twice, hence only half of the dirichlet condition has to be set
if (cParams.boundaryTypes1.isCouplingDirichlet(transportEqIdx1))
......
......@@ -177,13 +177,6 @@ public:
averagedNormal += boundaryFaceNormal;
}
// TODO: move scope below to coupling localoperator
// Beavers-Joseph condition at the coupling boundary/interface
if(bcTypes.hasCoupling() || bcTypes.hasCouplingMortar())
{
evalBeaversJoseph_(&intersection, scvIdx, boundaryFaceIdx, boundaryVars);
}
// TODO: move scope below to coupling localoperator/ BUG (potentially): sollte das nicht dirichlet sein?
// set velocity normal to the interface
if (bcTypes.isCouplingNeumann(momentumYIdx))
......@@ -312,6 +305,7 @@ protected:
}
template <class IntersectionIterator>
DUNE_DEPRECATED_MSG("evalBeaversJoseph_ is deprecated. Its functionality is now included in the LOP function evalCoupling21().")
void evalBeaversJoseph_(const IntersectionIterator &isIt,
const int scvIdx,
const int boundaryFaceIdx,
......
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -74,26 +74,15 @@ class TwoCNIStokesTwoPTwoCNISpatialParams : public ImplicitSpatialParams<TypeTag
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> DimVector;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef std::vector<Scalar> PermeabilityType;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
typedef std::vector<MaterialLawParams> MaterialLawParamsVector;
public:
/*!
* \brief Spatial parameters for the
* coupling of an isothermal two-component Stokes
......@@ -107,6 +96,7 @@ public:
porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Porosity);
permeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Permeability);
lambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LambdaSolid);
alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
// residual saturations
params_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
......@@ -206,10 +196,23 @@ public:
return lambdaSolid_;
}
/*!
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* \param globalPos The global position
*
* \return Beavers-Joseph coefficient
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return alphaBJ_;
}
private:
Scalar permeability_;
Scalar porosity_;
Scalar lambdaSolid_;
Scalar alphaBJ_;
MaterialLawParams params_;
};
......
......@@ -96,7 +96,6 @@ class Stokes2cniSubProblem : public StokesProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
// soil parameters for beavers & joseph
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
enum {
......@@ -176,7 +175,6 @@ public:
sinusTPeriod_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperaturePeriod);
useDirichletBC_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, FreeFlow, UseDirichletBC);
alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
}
......@@ -347,34 +345,6 @@ public:
}
}
/*!
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* \param globalPos The global position
*
* \return Beavers-Joseph coefficient
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return alphaBJ_;
}
/*!
* \brief Returns the intrinsic permeability tensor \f$[m^2]\f$
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
*/
Scalar permeability(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return spatialParams_.intrinsicPermeability(element,
fvGeometry,
scvIdx);
}
// \}
/*!
......@@ -519,7 +489,6 @@ private:
const Scalar height_() const
{ return bBoxMax_[1] - bBoxMin_[1]; }
// spatial parameters
SpatialParams spatialParams_;
static constexpr Scalar eps_ = 1e-8;
......@@ -542,7 +511,6 @@ private:
Scalar sinusTPeriod_;
bool useDirichletBC_;
Scalar alphaBJ_;
Scalar runUpDistanceX_;
Scalar initializationTime_;
......
......@@ -74,16 +74,15 @@ class TwoCNIZeroEqTwoPTwoCNISpatialParams : public ImplicitSpatialParams<TypeTag
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
public:
/*!
* \brief Spatial parameters for the
* coupling of a non-isothermal two-component ZeroEq
......@@ -97,6 +96,7 @@ public:
permeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Permeability);
porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Porosity);
thermalConductivitySolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, ThermalConductivitySolid);
alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
spatialParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
spatialParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Snr));
......@@ -192,10 +192,24 @@ public:
return thermalConductivitySolid_;
}
/*!
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* \param globalPos The global position
*
* \return Beavers-Joseph coefficient
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return alphaBJ_;
}
private:
Scalar permeability_;
Scalar porosity_;
Scalar thermalConductivitySolid_;
Scalar alphaBJ_;
MaterialLawParams spatialParams_;
};
} // end namespace
......
......@@ -270,32 +270,6 @@ public:
values = 0.;
}
/*!
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* \param globalPos The global position
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return alphaBJ_;
}
/*!
* \brief Evaluate the intrinsic permeability at the corner of a given element
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scvIdx The local subcontrolvolume index
*/
Scalar permeability(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return spatialParams_.intrinsicPermeability(element,
fvGeometry,
scvIdx);
}
//! \copydoc Dumux::ImplicitProblem::sourceAtPos()
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
......@@ -324,6 +298,14 @@ public:
|| (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)));
}
/*!
* \brief Returns the spatial parameters object.
*/
SpatialParams &spatialParams()
{ return spatialParams_; }
const SpatialParams &spatialParams() const
{ return spatialParams_; }
/*!
* \brief Auxiliary function used for the mortar coupling, if mortar coupling,
* this should return true
......
......@@ -77,26 +77,15 @@ class TwoCStokesTwoPTwoCSpatialParams : public ImplicitSpatialParams<TypeTag>
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> DimVector;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef std::vector<Scalar> PermeabilityType;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
typedef std::vector<MaterialLawParams> MaterialLawParamsVector;
public:
/*!
* \brief Spatial parameters for the
* coupling of an isothermal two-component Stokes
......@@ -109,7 +98,7 @@ public:
{
porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Porosity);
permeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Permeability);
lambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LambdaSolid);
alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
// residual saturations
params_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
......@@ -162,57 +151,21 @@ public:
}
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
*/
Scalar solidHeatCapacity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 790;
}
/*!
* \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* This is only required for non-isothermal models.
* \param globalPos The global position
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
*/
Scalar solidDensity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 2700; // density of granite [kg/m^3]
}
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
* \return Beavers-Joseph coefficient
*/
Scalar solidThermalConductivity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return lambdaSolid_;
return alphaBJ_;
}
private:
Scalar permeability_;
Scalar porosity_;
Scalar lambdaSolid_;
Scalar alphaBJ_;
MaterialLawParams params_;
};
......
......@@ -93,7 +93,6 @@ class Stokes2cSubProblem : public StokesProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
// soil parameters for beavers & joseph
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
......@@ -171,7 +170,6 @@ public:
sinusTAmplitude_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperatureAmplitude);
sinusTPeriod_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperaturePeriod);
alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
}
......@@ -316,33 +314,6 @@ public:
}
}
/*!
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* \param globalPos The global position
*
* \return Beavers-Joseph coefficient
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return alphaBJ_;
}
/*!
* \brief Returns the intrinsic permeability tensor \f$[m^2]\f$
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
*/
Scalar permeability(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return spatialParams_.intrinsicPermeability(element,
fvGeometry,
scvIdx);
}
// \}
/*!
......@@ -483,7 +454,6 @@ private:
const Scalar height_() const
{ return bBoxMax_[1] - bBoxMin_[1]; }
// spatial parameters
SpatialParams spatialParams_;
static constexpr Scalar eps_ = 1e-8;
......@@ -504,8 +474,6 @@ private:
Scalar sinusTAmplitude_;
Scalar sinusTPeriod_;
Scalar alphaBJ_;
Scalar runUpDistanceX_;
Scalar initializationTime_;
};
......
......@@ -74,16 +74,15 @@ class TwoCZeroEqTwoPTwoCSpatialParams : public ImplicitSpatialParams<TypeTag>
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
public:
/*!
* \brief Spatial parameters for the
* coupling of an isothermal two-component ZeroEq
......@@ -96,6 +95,7 @@ public:
{
permeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Permeability);
porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Porosity);
alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
spatialParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
spatialParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Snr));
......@@ -145,9 +145,22 @@ public:
return spatialParams_;
}
/*!
* \brief Evaluate the Beavers-Joseph coefficient at given position
*
* \param globalPos The global position
*
* \return Beavers-Joseph coefficient
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
return alphaBJ_;
}
private:
Scalar permeability_;
Scalar porosity_;
Scalar alphaBJ_;
MaterialLawParams spatialParams_;
};