Commit f563b6ae authored by Martin Utz's avatar Martin Utz
Browse files

[swe][frictionlaws] Enable test for Nikuradse

One can run the test now also with Nikuradse. But be aware that the
friction value in the input file is a Manning values.
Addintional changes
- The names of the variables used for the calculation of the
  boundaryfluxes are changed, because they were not really clear. Also the
  documentation of their file is improved.
- Some smaller modification were also applyied, which were demanded
  by the discussion about the merge request.
parent 28812f9f
......@@ -19,8 +19,13 @@
/*!
* \file
* \ingroup ShallowWaterModel
* \brief Compute boundary conditions (cell state) via Riemann invariants
* \brief Compute boundary conditions for the Riemann Solver
*
* The boundary conditions are given at the the outer face of the
* the boundary cells. In this form the boundary condition can't be
* processed by the riemann Solver, because it needs two cell states, one at
* each side of a face. Therefore they Riemann invariants are used to
* calculate a virtual outer state.
*/
#ifndef DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#define DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
......@@ -32,49 +37,69 @@ namespace Dumux {
namespace ShallowWater {
/*!
* \brief compute the cell state for fixed water depth boundary.
* \brief Compute the outer cell state for fixed water depth boundary.
*
* \param waterDepthBoundary Discharge per meter at the boundary face [m^2/s]
* \param waterDepthInside Water depth in the inner cell [m]
* \param velocityXInside Velocity in x-direction in the inner cell [m/s]
* \param velocityYInside Velocity in y-direction in the inner cell [m/s]
* \param gravity Gravity constant [m^3/s]
* \param nxy Normal vector of the boundary face
*
* \return cellStateOutside The outer cell state
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar, 3> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const Scalar waterDepthInside,
const Scalar velocityXInside,
const Scalar velocityYInside,
const Scalar gravity,
const GlobalPosition& nxy)
{
std::array<Scalar, 3> cellStateRight;
cellStateRight[0] = waterDepthBoundary;
std::array<Scalar, 3> cellStateOutside;
cellStateOutside[0] = waterDepthBoundary;
using std::sqrt;
const auto uboundIn = nxy[0] * velocityXLeft + nxy[1] * velocityYLeft;
const auto uboundQut = uboundIn + 2.0 * sqrt(9.81 * waterDepthLeft) - 2.0 * sqrt(9.81 * cellStateRight[0]);
const auto uboundIn = nxy[0] * velocityXInside + nxy[1] * velocityYInside;
const auto uboundQut = uboundIn + 2.0 * sqrt(gravity * waterDepthInside) - 2.0 * sqrt(gravity * cellStateOutside[0]);
cellStateRight[1] = (nxy[0] * uboundQut); // we only use the normal part
cellStateRight[2] = (nxy[1] * uboundQut); // we only use the normal part
cellStateOutside[1] = (nxy[0] * uboundQut); // we only use the normal part
cellStateOutside[2] = (nxy[1] * uboundQut); // we only use the normal part
return cellStateRight;
return cellStateOutside;
}
/*!
* \brief compute the cell state for a fixed discharge boundary.
* \brief Compute the outer cell state for a fixed discharge boundary.
*
* \param dischargeBoundary Discharge per meter at the boundary face [m^2/s]
* \param waterDepthInside Water depth in the inner cell [m]
* \param velocityXInside Velocity in x-direction in the inner cell [m/s]
* \param velocityYInside Velocity in y-direction in the inner cell [m/s]
* \param gravity Gravity constant [m^3/s]
* \param nxy Normal vector of the boundary face
*
* \return cellStateOutside The outer cell state
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar, 3> fixedDischargeBoundary(const Scalar qlocal,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
std::array<Scalar, 3> fixedDischargeBoundary(const Scalar dischargeBoundary,
const Scalar waterDepthInside,
const Scalar velocityXInside,
const Scalar velocityYInside,
const Scalar gravity,
const GlobalPosition& nxy)
{
std::array<Scalar, 3> cellStateRight;
std::array<Scalar, 3> cellStateOutside;
using std::abs;
using std::sqrt;
using std::max;
// only impose if abs(q) > 0
if (abs(qlocal) > 1.0e-9)
if (abs(dischargeBoundary) > 1.0e-9)
{
const auto uboundIn = nxy[0]*velocityXLeft + nxy[1]*velocityYLeft;
const auto alphal = uboundIn + 2.0*sqrt(9.81 * waterDepthLeft);
const auto uboundIn = nxy[0]*velocityXInside + nxy[1]*velocityYInside;
const auto alphal = uboundIn + 2.0*sqrt(gravity * waterDepthInside);
//initial guess for hstar solved with newton
constexpr Scalar tol_hstar = 1.0E-12;
......@@ -84,8 +109,8 @@ std::array<Scalar, 3> fixedDischargeBoundary(const Scalar qlocal,
Scalar hstar = 0.1;
for (int i = 0; i < maxstep_hstar; ++i)
{
Scalar f_hstar = alphal - qlocal/hstar - 2 * sqrt(9.81 * hstar);
Scalar df_hstar = (f_hstar -(alphal - qlocal/(hstar + ink_hstar) - 2 * sqrt(9.81 * (hstar+ink_hstar))))/ink_hstar;
Scalar f_hstar = alphal - dischargeBoundary/hstar - 2 * sqrt(gravity * hstar);
Scalar df_hstar = (f_hstar -(alphal - dischargeBoundary/(hstar + ink_hstar) - 2 * sqrt(gravity * (hstar+ink_hstar))))/ink_hstar;
Scalar dx_hstar = -f_hstar/df_hstar;
hstar = max(hstar - dx_hstar,0.001);
......@@ -93,13 +118,13 @@ std::array<Scalar, 3> fixedDischargeBoundary(const Scalar qlocal,
break;
}
const auto qinner = (nxy[0] * waterDepthLeft * velocityYLeft) - (nxy[1] * waterDepthLeft * velocityXLeft);
cellStateRight[0] = hstar;
cellStateRight[1] = (nxy[0] * qlocal - nxy[1] * qinner)/hstar;
cellStateRight[2] = (nxy[1] * qlocal + nxy[0] * qinner)/hstar;
const auto qinner = (nxy[0] * waterDepthInside * velocityYInside) - (nxy[1] * waterDepthInside * velocityXInside);
cellStateOutside[0] = hstar;
cellStateOutside[1] = (nxy[0] * dischargeBoundary - nxy[1] * qinner)/hstar;
cellStateOutside[2] = (nxy[1] * dischargeBoundary + nxy[0] * qinner)/hstar;
}
return cellStateRight;
return cellStateOutside;
}
} // end namespace ShallowWater
......
......@@ -45,7 +45,7 @@ public:
* \return ustar_h friction used for the source term in shallow water models.
*/
const virtual Scalar computeUstarH(const Scalar waterDept, const Scalar frictionValue) = 0;
virtual Scalar computeUstarH(const Scalar waterDept, const Scalar frictionValue) const = 0;
/*!
* \brief Limit the friction for small water depth.
......@@ -76,7 +76,7 @@ public:
* \param roughnessHeight roughness height of the representive structure (e.g. largest grain size).
* \param waterDepth water depth.
*/
Scalar limitRoughH(const Scalar roughnessHeight, const Scalar waterDepth)
Scalar limitRoughH(const Scalar roughnessHeight, const Scalar waterDepth) const
{
using std::min;
using std::max;
......
......@@ -50,7 +50,7 @@ public:
*
* \return ustar_h friction used for the source term in shallow water models.
*/
const Scalar computeUstarH(const Scalar waterDepth, const Scalar frictionValue) final
Scalar computeUstarH(const Scalar waterDepth, const Scalar frictionValue) const final
{
using std::pow;
......
......@@ -40,6 +40,8 @@ template <typename Scalar>
class FrictionLawNikuradse : public FrictionLaw<Scalar>
{
public:
FrictionLawNikuradse(const Scalar gravity)
: gravity_(gravity) {}
/*!
* \brief Compute the friction ustar_h.
*
......@@ -49,7 +51,7 @@ public:
* \return ustar_h friction used for the source term in shallow water models.
*/
const Scalar computeUstarH(const Scalar waterDepth, const Scalar frictionValue) final
Scalar computeUstarH(const Scalar waterDepth, const Scalar frictionValue) const final
{
using std::pow;
using std::log;
......@@ -61,6 +63,8 @@ public:
ustar_h = pow(0.41,2.0)/pow(log((12*(waterDepth + rough_h))/frictionValue),2.0);
return ustar_h;
}
private:
Scalar gravity_;
};
} // end namespace Dumux
......
......@@ -3,7 +3,7 @@ Name = roughchannel
FrictionValue = 0.025 # [-]
BedSlope = 0.001 # [-]
Gravity = 9.81 # [m/s^2]
Discharge = -1.0 # [m/s] discharge at the inflow boundary
Discharge = -1.0 # [m^2/s] discharge per meter at the inflow boundary
FrictionLaw = Manning
[TimeLoop]
......
......@@ -144,15 +144,21 @@ public:
discharge_ = getParam<Scalar>("Problem.Discharge");
hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
const auto gravity = this->spatialParams().gravity();
if (getParam<std::string>("Problem.FrictionLaw") == ("Manning"))
{
const auto gravity = this->spatialParams().gravity();
auto ptrManning = std::make_shared<FrictionLawManning<Scalar>>(gravity);
frictionLaw_ = std::static_pointer_cast<FrictionLaw<Scalar>>(ptrManning);
frictionLaw_ = std::make_shared<FrictionLawManning<Scalar>>(gravity);
}
else if (getParam<std::string>("Problem.FrictionLaw") == ("Nikuradse"))
{
frictionLaw_ = std::make_shared<FrictionLawNikuradse<Scalar>>(gravity);
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<<"Change the FrictionLaw in params.input. This test is only valid for Manning!";
std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are 'Manning' and 'Nikuradse'!";
}
}
......@@ -175,7 +181,7 @@ public:
using std::abs;
using std::sqrt;
return std::pow(std::abs(discharge)*manningN/sqrt(bedSlope), 0.6);
return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6);
}
//! Udpate the analytical solution
......@@ -306,6 +312,7 @@ public:
insideVolVars.waterDepth(),
insideVolVars.velocity(0),
insideVolVars.velocity(1),
gravity,
nxy);
}
// impose water depth at the right side
......@@ -315,6 +322,7 @@ public:
insideVolVars.waterDepth(),
insideVolVars.velocity(0),
insideVolVars.velocity(1),
gravity,
nxy);
}
// no flow boundary
......
......@@ -96,7 +96,8 @@ public:
*
* \param element The current element
* \param scv The sub-control volume inside the element.
* \return the bed surface
*
* \return The bed surface
*/
Scalar bedSurface(const Element& element,
const SubControlVolume& scv) const
......
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