Commit 5240abd6 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[multidomain]

cleaned up the stokes-darcy models:
  - SET_TYPE_PROP task from Dumux Day
  - cleaned up input
  - removed commented lines of code
  - removed try/catch blocks

updated the coupled-freeflow reference results which was broken in one
of the previous commits

reviewed by gruenich



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14094 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 925d4ce3
......@@ -26,9 +26,6 @@
#include <dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh>
//#include <dumux/freeflow/stokesncni/stokesncnimodel.hh>
//#include <dumux/implicit/2p2cnicoupling/2p2cnicouplingmodel.hh>
namespace Dumux {
/*!
......
......@@ -133,9 +133,9 @@ class TwoCStokesTwoPTwoCLocalOperator :
blModel_ = GET_PARAM_FROM_GROUP(TypeTag, int, FreeFlow, BoundaryLayerModel);
massTransferModel_ = GET_PARAM_FROM_GROUP(TypeTag, int, FreeFlow, MassTransferModel);
if (blModel_ > 0)
if (blModel_ != 0)
std::cout << "Using boundary layer model " << blModel_ << std::endl;
if (massTransferModel_ > 0)
if (massTransferModel_ != 0)
std::cout << "Using mass transfer model " << massTransferModel_ << std::endl;
}
......@@ -551,7 +551,7 @@ class TwoCStokesTwoPTwoCLocalOperator :
const Scalar boundaryLayerOffset =
GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, BoundaryLayerOffset);
return 5*(globalPos[0]+boundaryLayerOffset) / sqrt(reynoldsX);
return 5*(globalPos[0]+boundaryLayerOffset) / std::sqrt(reynoldsX);
}
if (blModel_ == 2)
{
......@@ -562,7 +562,7 @@ class TwoCStokesTwoPTwoCLocalOperator :
const Scalar boundaryLayerOffset =
GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, BoundaryLayerOffset);
return 0.37*(globalPos[0]+boundaryLayerOffset) / pow(reynoldsX, 0.2);
return 0.37*(globalPos[0]+boundaryLayerOffset) / std::pow(reynoldsX, 0.2);
}
if (blModel_ == 3)
{
......@@ -573,9 +573,9 @@ class TwoCStokesTwoPTwoCLocalOperator :
const Scalar boundaryLayerOffset =
GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, BoundaryLayerOffset);
const Scalar cf = 2*pow(0.41*1.5/log(reynoldsX),2);
const Scalar cf = 2*std::pow(0.41*1.5/std::log(reynoldsX),2);
return 50*(globalPos[0]+boundaryLayerOffset)/(reynoldsX*sqrt(cf/2));
return 50*(globalPos[0]+boundaryLayerOffset)/(reynoldsX*std::sqrt(cf/2));
}
if (blModel_ == 9)
{
......@@ -603,33 +603,28 @@ class TwoCStokesTwoPTwoCLocalOperator :
const DimVector& globalPos,
const int vertInElem1, const int vertInElem2) const
{
if (massTransferModel_ == 1){
Scalar exponentMTC = 0;
if (ParameterTree::tree().hasKey("FreeFlow.ExponentMTC"))
exponentMTC = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, ExponentMTC);
else
std::cerr << "FreeFlow.ExponentMTC not defined\n";
return pow(cParams.elemVolVarsCur2[vertInElem2].saturation(wPhaseIdx2), exponentMTC);
if (massTransferModel_ == 1)
{
Scalar exponentMTC = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, ExponentMTC);
return std::pow(cParams.elemVolVarsCur2[vertInElem2].saturation(wPhaseIdx2), exponentMTC);
}
// Schlünder model (Schlünder, CES 1988)
if (massTransferModel_ == 2){
Scalar charPoreRadius = 0;
if (ParameterTree::tree().hasKey("PorousMedium.CharPoreDiameter"))
charPoreRadius = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, CharPoreDiameter);
else
std::cerr << "PorousMedium.PoreDiameter not defined\n";
else if (massTransferModel_ == 2)
{
Scalar charPoreRadius = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, CharPoreDiameter);
const Scalar blThickness = computeBoundaryLayerThickness(cParams, globalPos, vertInElem1);
const Scalar moistureContent = cParams.elemVolVarsCur2[vertInElem2].saturation(wPhaseIdx2) *
globalProblem_.sdProblem2().spatialParams().porosity(sdElement2, cParams.fvGeometry2, vertInElem2);
Scalar massTransferCoeff = 1. + 2./M_PI * charPoreRadius/blThickness
* sqrt(M_PI/(4*moistureContent)) * (sqrt(M_PI/(4*moistureContent)) - 1.);
* std::sqrt(M_PI/(4*moistureContent)) * (std::sqrt(M_PI/(4*moistureContent)) - 1.);
return 1./massTransferCoeff;
}
// Schlünder model (Schlünder, CES 1988) with variable char. pore diameter depending on Pc
if (massTransferModel_ == 3){
else if (massTransferModel_ == 3)
{
const Scalar surfaceTension = 72.85e-3; // source: Wikipedia
const Scalar charPoreRadius = 2*surfaceTension/cParams.elemVolVarsCur2[vertInElem2].capillaryPressure();
......@@ -638,19 +633,14 @@ class TwoCStokesTwoPTwoCLocalOperator :
globalProblem_.sdProblem2().spatialParams().porosity(sdElement2, cParams.fvGeometry2, vertInElem2);
Scalar massTransferCoeff = 1. + 2./M_PI * charPoreRadius/blThickness
* sqrt(M_PI/(4*moistureContent)) * (sqrt(M_PI/(4*moistureContent)) - 1.);
* std::sqrt(M_PI/(4*moistureContent)) * (std::sqrt(M_PI/(4*moistureContent)) - 1.);
return 1./massTransferCoeff;
}
// modified Schlünder model
if (massTransferModel_ == 4){
// const Scalar surfaceTension = 72.85e-3; // source: Wikipedia
// const Scalar charPoreRadius = 2*surfaceTension/cParams.elemVolVarsCur2[vertInElem2].capillaryPressure();
Scalar charPoreRadius = 0;
if (ParameterTree::tree().hasKey("PorousMedium.CharPoreRadius"))
charPoreRadius = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, CharPoreDiameter);
else
std::cerr << "PorousMedium.CharPoreRadius not defined\n";
else if (massTransferModel_ == 4)
{
Scalar charPoreRadius = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, CharPoreDiameter);
const Scalar blThickness = computeBoundaryLayerThickness(cParams, globalPos, vertInElem1);
const Scalar moistureContent = cParams.elemVolVarsCur2[vertInElem2].saturation(wPhaseIdx2) *
......
......@@ -62,16 +62,16 @@ SET_TYPE_PROP(TwoCStokesTwoPTwoC, JacobianAssembler, Dumux::MultiDomainAssembler
// Specif the used Newton controller
SET_TYPE_PROP(TwoCStokesTwoPTwoC, NewtonController, Dumux::TwoCStokesTwoPTwoCNewtonController<TypeTag>);
// set this to one here (must fit to the structure of the coupled matrix which has block length 1)
// Set this to one here (must fit to the structure of the coupled matrix which has block length 1)
SET_INT_PROP(TwoCStokesTwoPTwoC, NumEq, 1);
// Write the solutions of individual newton iterations?
// Specitfy if the solutions of individual newton iterations should be written
SET_BOOL_PROP(TwoCStokesTwoPTwoC, NewtonWriteConvergence, false);
//! 0 = none, 1 = Blasius, 2 and 3 = turbulent BL, 9 = constant thickness
// 0 = none, 1 = Blasius, 2 and 3 = turbulent BL, 9 = constant thickness
SET_INT_PROP(TwoCStokesTwoPTwoC, FreeFlowBoundaryLayerModel, 0);
//! 0 = none, 1 = power law, 2 = Schluender model mass transfer model
// 0 = none, 1 = power law, 2 = Schluender model mass transfer model
SET_INT_PROP(TwoCStokesTwoPTwoC, FreeFlowMassTransferModel, 0);
} // end namespace properties
......
......@@ -86,16 +86,6 @@ public:
TwoPTwoCNICouplingLocalResidual()
{ };
// void evalBoundary_()
// {
// ParentType::evalBoundary_();
//
// // evaluate Dirichlet-like coupling conditions
// for (int idx = 0; idx < this->fvGeometry_().numVertices; idx++)
// if (boundaryHasCoupling_(this->bcTypes_(idx)))
// evalCouplingVertex_(idx);
// }
/*!
* \brief Implementation of the boundary evaluation
*/
......
......@@ -56,7 +56,6 @@ namespace Dumux
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
......
......@@ -40,7 +40,7 @@
#include <dumux/multidomain/common/multidomainproblem.hh>
#include <dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh>
#include <dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnipropertydefaults.hh>
//#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/seqsolverbackend.hh>
#ifdef HAVE_PARDISO
#include <dumux/linear/pardisobackend.hh>
......@@ -60,17 +60,13 @@ namespace Properties
NEW_TYPE_TAG(TwoCNIStokesTwoPTwoCNIProblem, INHERITS_FROM(TwoCNIStokesTwoPTwoCNI));
// Set the grid type
SET_PROP(TwoCNIStokesTwoPTwoCNIProblem, Grid)
{
public:
#ifdef HAVE_UG
typedef typename Dune::UGGrid<2> type;
SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, Grid, Dune::UGGrid<2>);
#elif HAVE_ALUGRID || HAVE_DUNE_ALUGRID
typedef typename Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming> type;
SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>;
#else
typedef typename Dune::YaspGrid<2> type;
SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, Grid, Dune::YaspGrid<2>);
#endif
};
// Set the global problem
SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, Problem, TwoCNIStokesTwoPTwoCNIProblem<TypeTag>);
......@@ -91,11 +87,9 @@ SET_TYPE_PROP(TwoPTwoCNISubProblem, MultiDomainTypeTag, TTAG(TwoCNIStokesTwoPTwo
SET_TYPE_PROP(Stokes2cniSubProblem, OtherSubDomainTypeTag, TTAG(TwoPTwoCNISubProblem));
SET_TYPE_PROP(TwoPTwoCNISubProblem, OtherSubDomainTypeTag, TTAG(Stokes2cniSubProblem));
// Set the same spatial parameters for both sub-problems
SET_PROP(Stokes2cniSubProblem, SpatialParams)
{ typedef Dumux::TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag> type; };
SET_PROP(TwoPTwoCNISubProblem, SpatialParams)
{ typedef Dumux::TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag> type; };
// Set the spatial parameters used for the problems
SET_TYPE_PROP(Stokes2cniSubProblem, SpatialParams, Dumux::TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag>);
SET_TYPE_PROP(TwoPTwoCNISubProblem, SpatialParams, Dumux::TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag>);
// Set the fluid system
SET_PROP(TwoCNIStokesTwoPTwoCNIProblem, FluidSystem)
......@@ -228,29 +222,17 @@ public:
TimeManager &timeManager)
: ParentType(mdGrid, timeManager)
{
try
{
// define location of the interface
interface_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
noDarcyX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, NoDarcyX);
episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
dtInit_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
// define output options
freqRestart_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqRestart);
freqOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqOutput);
freqFluxOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqFluxOutput);
freqVaporFluxOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqVaporFluxOutput);
}
catch (Dumux::ParameterException &e) {
std::cerr << e << ". Abort!\n";
exit(1) ;
}
catch (...) {
std::cerr << "Unknown exception thrown!\n";
exit(1);
}
interfacePosY_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
noDarcyX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX);
episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
dtInit_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
// define output options
freqRestart_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqRestart);
freqOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqOutput);
freqFluxOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqFluxOutput);
freqVaporFluxOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqVaporFluxOutput);
stokes2cni_ = this->sdID1();
twoPtwoCNI_ = this->sdID2();
......@@ -286,9 +268,6 @@ public:
{
ParentType::init();
// plot the Pc-Sw curves, if requested
this->sdProblem2().spatialParams().plotMaterialLaw();
std::cout << "Writing flux data at interface\n";
if (this->timeManager().time() == 0)
{
......@@ -332,7 +311,7 @@ public:
GlobalPosition globalPos = elementIt->geometry().center();
if (globalPos[1] > interface_)
if (globalPos[1] > interfacePosY_)
mdGrid.addToSubDomain(stokes2cni_,*elementIt);
else
if(globalPos[0] > noDarcyX_)
......@@ -708,7 +687,6 @@ public:
return false;
}
/*!
* \brief Returns true if the current solution should be written to
* disk (i.e. as a VTK file)
......@@ -750,7 +728,6 @@ public:
|| this->timeManager().episodeWillBeOver())
return true;
return false;
}
/*!
......@@ -779,7 +756,7 @@ private:
unsigned freqFluxOutput_;
unsigned freqVaporFluxOutput_;
Scalar interface_;
Scalar interfacePosY_;
Scalar noDarcyX_;
Scalar episodeLength_;
Scalar initializationTime_;
......@@ -808,6 +785,6 @@ private:
std::ofstream fluxFile_;
};
} //end namespace
} // namespace Dumux
#endif // DUMUX_2CNISTOKES2P2CNIPROBLEM_HH
......@@ -33,11 +33,9 @@
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
//#include <dumux/io/plotfluidmatrixlaw.hh>
namespace Dumux
{
//forward declaration
template<class TypeTag>
class TwoCNIStokesTwoPTwoCNISpatialParams;
......@@ -49,27 +47,19 @@ NEW_TYPE_TAG(TwoCNIStokesTwoPTwoCNISpatialParams);
// Set the spatial parameters
SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNISpatialParams, SpatialParams,
TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag>);
TwoCNIStokesTwoPTwoCNISpatialParams<TypeTag>);
// Set the material Law
SET_PROP(TwoCNIStokesTwoPTwoCNISpatialParams, MaterialLaw)
{
private:
// define the material law which is parameterized by effective
// saturations
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef RegularizedVanGenuchten<Scalar> EffMaterialLaw;
public:
// define the material law parameterized by absolute saturations
typedef EffToAbsLaw<EffMaterialLaw> type;
};
// Set the material law parametrized by absolute saturations
SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNISpatialParams,
MaterialLaw,
EffToAbsLaw<RegularizedVanGenuchten<typename GET_PROP_TYPE(TypeTag, Scalar)>>);
// EffToAbsLaw<RegularizedBrooksCorey<typename GET_PROP_TYPE(TypeTag, Scalar)> >);
}
/*!
* \ingroup TwoPTwoCNiModel
* \ingroup StokesniModel
* \ingroup ImplicitTestProblems
* \ingroup MultidomainProblems
* \brief Definition of the spatial parameters for
* the coupling of an non-isothermal two-component Stokes
* and an non-isothermal two-phase two-component Darcy model.
......@@ -116,58 +106,18 @@ public:
* \param gridView The GridView which is used by the problem
*/
TwoCNIStokesTwoPTwoCNISpatialParams(const GridView& gridView)
: ParentType(gridView),
permeability_(gridView.size(dim), 0.0),
vanGenuchtenAlpha_(gridView.size(dim), 0.0),
indexSet_(gridView.indexSet())
: ParentType(gridView)
{
try
{
soilType_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, SpatialParams, SoilType);
xMaterialInterface_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, MaterialInterfaceX);
// porosities
coarsePorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Porosity1);
mediumPorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Porosity2);
finePorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Porosity3);
// intrinsic permeabilities
coarsePermeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Permeability1);
mediumPermeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Permeability2);
finePermeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Permeability3);
// thermal conductivity of the solid material
coarseLambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, LambdaSolid1);
mediumLambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, LambdaSolid2);
fineLambdaSolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, LambdaSolid3);
if (soilType_ != 0)
{
// residual saturations
coarseParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Swr1));
coarseParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Snr1));
mediumParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2));
mediumParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
fineParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Swr3));
fineParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Snr3));
// parameters for the vanGenuchten law
coarseParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, VgAlpha1));
coarseParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, VgN1));
mediumParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, VgAlpha2));
mediumParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, VgN2));
fineParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, VgAlpha3));
fineParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, VgN3));
}
}
catch (Dumux::ParameterException &e) {
std::cerr << e << ". Abort!\n";
exit(1) ;
}
catch (...) {
std::cerr << "Unknown exception thrown!\n";
exit(1);
}
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);
// residual saturations
params_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
params_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Snr));
// parameters for the vanGenuchten law
params_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, VgAlpha));
params_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, VgN));
}
/*!
......@@ -187,18 +137,7 @@ public:
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
if (checkSoilType(globalPos) == 1)
return coarsePermeability_;
if (checkSoilType(globalPos) == 2)
return mediumPermeability_;
if (checkSoilType(globalPos) == 3)
return finePermeability_;
if (checkSoilType(globalPos) == 4)
return finePermeability_;
else
return mediumPermeability_;
return permeability_;
}
/*!
......@@ -212,19 +151,9 @@ public:
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
if (checkSoilType(globalPos) == 1)
return coarsePorosity_;
if (checkSoilType(globalPos) == 2)
return mediumPorosity_;
if (checkSoilType(globalPos) == 3)
return finePorosity_;
else
return mediumPorosity_;
return porosity_;
}
/*!
* \brief Returns the parameter object for the material law
*
......@@ -236,20 +165,9 @@ public:
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
if (checkSoilType(globalPos)==1)
return coarseParams_;
if (checkSoilType(globalPos)==2)
return mediumParams_;
if (checkSoilType(globalPos)==3)
return fineParams_;
// if (checkSoilType(globalPos)==4)
// return leverettJParams_;
else
return mediumParams_;
return params_;
}
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
*
......@@ -263,7 +181,7 @@ public:
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 790; // specific heat capacity of granite [J / (kg K)]
return 790;
}
/*!
......@@ -295,117 +213,16 @@ public:
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
const GlobalPosition &pos = element.geometry().corner(scvIdx);
if (checkSoilType(pos) == 1)
return coarseLambdaSolid_;
if (checkSoilType(pos) == 2)
return mediumLambdaSolid_;
if (checkSoilType(pos) == 3)
return fineLambdaSolid_;
else
return mediumLambdaSolid_;
}
/*!
* \brief Returns the index of the used soil type
*
* The soil, can be chosen as runtime parameter:
* 1: coarse,
* 2: medium,
* 3: fine,
* 4: LeverettJ (x < xMaterialInterface)
*
* \param pos The global position
*/
const unsigned checkSoilType(const GlobalPosition &pos) const
{
return soilType_;
}
/*!
* \brief This is called from the coupled problem and creates
* a gnuplot output of the Pc-Sw curve
*
* If this function should be used, uncomment the lines between
* the curly brackets.
*/
void plotMaterialLaw()
{
// if (soilType_ == 0)
// {
// std::cout << "Material law plot not possible for heterogeneous media!\n";
// return;
// }
// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.Coarse, PlotMaterialLaw1))
// {
// PlotFluidMatrixLaw<MaterialLaw> coarseFluidMatrixLaw_;
// coarseFluidMatrixLaw_.plotpC(coarseParams_,
// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Swr1),
// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Coarse, Snr1),
// "pcSw_coarse");
// }
// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.Medium, PlotMaterialLaw2))
// {
// PlotFluidMatrixLaw<MaterialLaw> mediumFluidMatrixLaw_;
// mediumFluidMatrixLaw_.plotpC(mediumParams_,
// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2),
// "pcSw_medium");
// PlotFluidMatrixLaw<MaterialLaw> mediumFluidMatrixLaw2_;
// mediumFluidMatrixLaw_.plotkrw(mediumParams_,
// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
// PlotFluidMatrixLaw<MaterialLaw> mediumFluidMatrixLaw3_;
// mediumFluidMatrixLaw_.plotkrn(mediumParams_,
// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2));
// }
// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.Fine, PlotMaterialLaw3))
// {
// PlotFluidMatrixLaw<MaterialLaw> fineFluidMatrixLaw_;
// fineFluidMatrixLaw_.plotpC(fineParams_,
// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Swr3),
// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Fine, Snr3),
// "pcSw_fine");
// }
// if (GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams.LeverettJ, PlotMaterialLawJ))
// {
// PlotFluidMatrixLaw<MaterialLaw> leverettJFluidMatrixLaw_;
// leverettJFluidMatrixLaw_.plotpC(leverettJParams_,
// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Swr2),
// 1.0 - GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams.Medium, Snr2),
// "pcSw_leverett");
// }
return lambdaSolid_;
}
private:
unsigned soilType_;