Commit 8e37ca8b authored by Martin Utz's avatar Martin Utz
Browse files

[swe][frictionlaws] Use std::vector<FrictionLaw>

For each element a object of a child class of FrictionLaw is
instantiated and stored within a vector.
This makes it possible to change the friction law on runtime, use also friction
laws which have more than one parameter, use spatial variying friction
values and even to use muliple friction laws.
parent afef51b0
......@@ -66,7 +66,7 @@ private:
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = RoughChannelSpatialParams<FVGridGeometry, Scalar>;
using type = RoughChannelSpatialParams<FVGridGeometry, Scalar, TypeTag>;
};
template<class TypeTag>
......@@ -144,23 +144,6 @@ public:
bedSlope_ = getParam<Scalar>("Problem.BedSlope");
discharge_ = getParam<Scalar>("Problem.Discharge");
hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
const auto gravity = this->spatialParams().gravity();
if (getParam<std::string>("Problem.FrictionLaw") == ("Manning"))
{
frictionLaw_ = std::make_shared<FrictionLawManning<Scalar, VolumeVariables>>(gravity, constManningN_);
}
else if (getParam<std::string>("Problem.FrictionLaw") == ("Nikuradse"))
{
frictionLaw_ = std::make_shared<FrictionLawNikuradse<Scalar, VolumeVariables>>(0.003);
std::cout<<"\nWARNING: This test is meant to be run for the friction law after Manning.\n";
std::cout<<" You are running it with Nikuradse, although the friction values in\n";
std::cout<<" the input file and the analytic solution don't fit to Nikuradse!\n\n.";
}
else
{
std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are 'Manning' and 'Nikuradse'!";
}
}
//! Get the analytical water depth
......@@ -234,17 +217,16 @@ public:
* that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
NumEqVector source(const Element &element,
NumEqVector source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
const auto& globalPos = scv.center();
const auto& volVars = elemVolVars[scv];
NumEqVector source (0.0);
Dune::FieldVector<Scalar, 2> bottomShearStress = frictionLaw_->computeShearStress(volVars);
Dune::FieldVector<Scalar, 2> bottomShearStress = this->spatialParams().frictionLaw(element)->computeShearStress(volVars);
source += computeBottomFrictionSource(bottomShearStress);
return source;
......@@ -395,7 +377,6 @@ private:
Scalar constManningN_; // analytic solution is only available for const friction.
Scalar bedSlope_;
Scalar discharge_; // discharge at the inflow boundary
std::shared_ptr<FrictionLaw<Scalar, VolumeVariables>> frictionLaw_;
static constexpr Scalar eps_ = 1.0e-6;
std::string name_;
};
......
......@@ -26,6 +26,9 @@
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/common/parameters.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh>
namespace Dumux {
......@@ -34,18 +37,20 @@ namespace Dumux {
* \brief The spatial parameters class for the rough channel test.
*
*/
template<class FVGridGeometry, class Scalar>
template<class FVGridGeometry, class Scalar, class TypeTag>
class RoughChannelSpatialParams
: public FVSpatialParams<FVGridGeometry, Scalar,
RoughChannelSpatialParams<FVGridGeometry, Scalar>>
RoughChannelSpatialParams<FVGridGeometry, Scalar, TypeTag>>
{
using ThisType = RoughChannelSpatialParams<FVGridGeometry, Scalar>;
using ThisType = RoughChannelSpatialParams<FVGridGeometry, Scalar, TypeTag>;
using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
using GridView = typename FVGridGeometry::GridView;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
public:
RoughChannelSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
......@@ -54,6 +59,34 @@ public:
gravity_ = getParam<Scalar>("Problem.Gravity");
bedSlope_ = getParam<Scalar>("Problem.BedSlope");
frictionValue_ = getParam<Scalar>("Problem.FrictionValue");
frictionLawType_ = getParam<std::string>("Problem.FrictionLaw");
initFrictionLaw(fvGridGeometry);
}
/*!
* \brief Initialize FrictionLaws
*/
void initFrictionLaw(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
{
if (frictionLawType_ == "Manning")
{
for (const auto& element : elements(fvGridGeometry->gridView()))
{
frictionLaw_.push_back(std::make_shared<FrictionLawManning<Scalar, VolumeVariables>>(gravity_, frictionValue_));
}
}
if (frictionLawType_ == "Nikuradse")
{
for (const auto& element : elements(fvGridGeometry->gridView()))
{
frictionLaw_.push_back(std::make_shared<FrictionLawNikuradse<Scalar, VolumeVariables>>(frictionValue_));
}
}
else
{
std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are 'Manning' and 'Nikuradse'!";
}
}
/*! \brief Define the gravitation.
......@@ -90,6 +123,11 @@ public:
/*! \brief Define the friction value.
*
* The unit of the friction value depends on the used friction law.
* The Mannng friction value is just a coefficient and has no unit $[-]$
* In the Nikuradse law the friction value is the equivalent sand roughness
* with the unit $[m]$.
*
* \return friction value
*/
Scalar frictionValue() const
......@@ -97,6 +135,18 @@ public:
return frictionValue_;
}
/*! \brief Get the frictionLaw.
*
* Get the frictionLaw, which already includes the friction value.
*
* \return frictionLaw
*/
std::shared_ptr<FrictionLaw<Scalar, VolumeVariables>> frictionLaw(const Element element) const
{
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
return frictionLaw_[eIdx];
}
/*! \brief Define the bed surface
*
* \param element The current element
......@@ -115,6 +165,8 @@ private:
Scalar gravity_;
Scalar bedSlope_;
Scalar frictionValue_;
std::string frictionLawType_;
std::vector<std::shared_ptr<FrictionLaw<Scalar, VolumeVariables>>> frictionLaw_;
};
} // end namespace Dumux
......
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