From f563b6ae0f2c5d04a3decb859b74ed55d9605882 Mon Sep 17 00:00:00 2001
From: Martin Utz <martin.utz@baw.de>
Date: Thu, 13 Jun 2019 12:42:24 +0200
Subject: [PATCH] [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.
---
dumux/freeflow/shallowwater/boundaryfluxes.hh | 81 ++++++++++++-------
.../frictionlaws/frictionlaw.hh | 4 +-
.../frictionlaws/manning.hh | 2 +-
.../frictionlaws/nikuradse.hh | 6 +-
.../shallowwater/roughchannel/params.input | 2 +-
.../shallowwater/roughchannel/problem.hh | 18 +++--
.../roughchannel/spatialparams.hh | 3 +-
7 files changed, 77 insertions(+), 39 deletions(-)
diff --git a/dumux/freeflow/shallowwater/boundaryfluxes.hh b/dumux/freeflow/shallowwater/boundaryfluxes.hh
index 178db2c80b..97833bf74f 100644
--- a/dumux/freeflow/shallowwater/boundaryfluxes.hh
+++ b/dumux/freeflow/shallowwater/boundaryfluxes.hh
@@ -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
diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh
index c66da30da8..fb57e2af46 100644
--- a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh
+++ b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh
@@ -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;
diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh
index 53062431cf..db06d92b83 100644
--- a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh
+++ b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh
@@ -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;
diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh
index 87767d975a..ec3a9d91b3 100644
--- a/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh
+++ b/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh
@@ -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
diff --git a/test/freeflow/shallowwater/roughchannel/params.input b/test/freeflow/shallowwater/roughchannel/params.input
index 1afacd7100..3ab64984f0 100644
--- a/test/freeflow/shallowwater/roughchannel/params.input
+++ b/test/freeflow/shallowwater/roughchannel/params.input
@@ -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]
diff --git a/test/freeflow/shallowwater/roughchannel/problem.hh b/test/freeflow/shallowwater/roughchannel/problem.hh
index 4198c27290..501e0ee6d7 100644
--- a/test/freeflow/shallowwater/roughchannel/problem.hh
+++ b/test/freeflow/shallowwater/roughchannel/problem.hh
@@ -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
diff --git a/test/freeflow/shallowwater/roughchannel/spatialparams.hh b/test/freeflow/shallowwater/roughchannel/spatialparams.hh
index 888d168765..e7b6d17f46 100644
--- a/test/freeflow/shallowwater/roughchannel/spatialparams.hh
+++ b/test/freeflow/shallowwater/roughchannel/spatialparams.hh
@@ -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
--
GitLab