Commit 33f0c477 authored by Martin Utz's avatar Martin Utz
Browse files

[swe][frictionlaws] Include some small changes

- use std::array instead of std::vector.
- In general manningN is a spatial parameter, so put it there, even
  if it is a constant in the rough channel test.
- Rename the water depth from h to waterDepth.
- Fix the dambreak test, which was incidentially broken by the work on thei
  rouch channel test.
parent 0468997f
......@@ -25,7 +25,7 @@
#ifndef DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#define DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#include <vector>
#include <array>
#include <cmath>
namespace Dumux {
......@@ -35,14 +35,14 @@ namespace ShallowWater {
* \brief compute the cell state for fixed water depth boundary.
*/
template<class Scalar, class GlobalPosition>
std::vector<Scalar> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
std::array<Scalar, 3> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
{
std::vector<Scalar> cellStateRight(3);
std::array<Scalar, 3> cellStateRight;
cellStateRight[0] = waterDepthBoundary;
using std::sqrt;
......@@ -59,13 +59,13 @@ std::vector<Scalar> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
* \brief compute the cell state for a fixed discharge boundary.
*/
template<class Scalar, class GlobalPosition>
std::vector<Scalar> fixedDischargeBoundary(const Scalar qlocal,
std::array<Scalar, 3> fixedDischargeBoundary(const Scalar qlocal,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
{
std::vector<Scalar> cellStateRight(3);
std::array<Scalar, 3> cellStateRight;
using std::abs;
using std::sqrt;
using std::max;
......
......@@ -45,25 +45,23 @@ public:
/*!
* \brief Compute the friction ustar_h.
*
* \param h water depth.
* \param waterDepth water depth.
* \param manningN Mannings friction value.
* \return ustar_h friction used for the source term in shallow water models.
*/
Scalar computeUstarH(const Scalar h,const Scalar manningN, const Scalar gravity)
Scalar computeUstarH(const Scalar waterDepth, const Scalar manningN, const Scalar gravity)
{
using std::pow;
Scalar ustar_h = 0.0;
Scalar rough_h = pow(25.68/(1.0/manningN),6.0);
rough_h = this->limitRoughH(rough_h, h);
rough_h = this->limitRoughH(rough_h, waterDepth);
auto cfric = pow((h + rough_h),1.0/6.0) * 1.0/(manningN);
auto cfric = pow((waterDepth + rough_h),1.0/6.0) * 1.0/(manningN);
ustar_h = gravity / pow(cfric,2.0);
return ustar_h;
}
};
} // end namespace Dumux
......
......@@ -44,12 +44,12 @@ public:
/*!
* \brief Compute the friction ustar_h.
*
* \param h water depth.
* \param waterDepth water depth.
* \param ks the Strickler friction value.
* \return ustar_h friction used for the source term in shallow water models.
*/
Scalar computeUstarH(const Scalar h,const Scalar ks)
Scalar computeUstarH(const Scalar waterDepth,const Scalar ks)
{
using std::pow;
using std::log;
......@@ -57,8 +57,8 @@ public:
Scalar ustar_h = 0.0;
Scalar rough_h = ks;
rough_h = limitRoughH(rough_h, h);
ustar_h = pow(0.41,2.0)/pow(log((12*(h + rough_h))/ks),2.0);
rough_h = limitRoughH(rough_h, waterDepth);
ustar_h = pow(0.41,2.0)/pow(log((12*(waterDepth + rough_h))/ks),2.0);
return ustar_h;
}
......@@ -84,9 +84,9 @@ public:
* the three empirical paramaters L, E and T, which describe the limiting curve.
*
* \param rough_h roughness height of the representive structure (e.g. largest grain size).
* \param h water depth.
* \param waterDepth water depth.
*/
Scalar limitRoughH(Scalar rough_h, const Scalar h)
Scalar limitRoughH(Scalar rough_h, const Scalar waterDepth)
{
using std::pow;
using std::min;
......@@ -98,7 +98,7 @@ public:
Scalar mobility_max = 1.0; //!< maximal mobility
auto minUpperH = rough_h * 2.0;
auto sw = min(h * (1.0/minUpperH),1.0);
auto sw = min(waterDepth * (1.0/minUpperH),1.0);
sw = max(0.0,sw);
auto mobility = (mobility_max * pow(sw,letL))/(pow(sw,letL) + letE * pow(1.0-sw,letT));
return rough_h * (1.0 - mobility);
......
......@@ -5,7 +5,7 @@ dune_add_test(NAME test_shallowwater_dambreak
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_dambreak-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/dambreak-00001.vtu
${CMAKE_CURRENT_BINARY_DIR}/damBreak-00001.vtu
--zeroThreshold {"velocityY":1e-14}
--command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_dambreak params.input
-Problem.Name dambreak")
-Problem.Name damBreak")
......@@ -238,19 +238,6 @@ public:
const auto& nxy = scvf.unitOuterNormal();
const auto gravity = this->spatialParams().gravity(scvf.center());
//left side q_in
if (scfv.center()[0] < 0.0 + eps_)
{
}
//right side fixed h
else if (scfv.center()[0] > 100.0 - eps_)
// no flow boundary
auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
insideVolVars.waterDepth(),
insideVolVars.velocity(0),
......
......@@ -5,7 +5,7 @@ dune_add_test(NAME test_shallowwater_roughchannel
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_roughchannel-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/roughchannel-00001.vtu
${CMAKE_CURRENT_BINARY_DIR}/roughChannel-00001.vtu
--zeroThreshold {"velocityY":1e-14}
--command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_roughchannel params.input
-Problem.Name roughchannel")
-Problem.Name roughChannel")
......@@ -110,7 +110,7 @@ int main(int argc, char** argv) try
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
problem->updateAnalyticalSolution(x,*gridVariables,0.0);
problem->updateAnalyticalSolution(x,*gridVariables);
IOFields::initOutputModule(vtkWriter);
vtkWriter.write(0.0);
......@@ -141,7 +141,7 @@ int main(int argc, char** argv) try
nonLinearSolver.solve(x,*timeLoop);
// update the analytical solution
problem->updateAnalyticalSolution(x,*gridVariables,timeLoop->time()+timeLoop->timeStepSize());
problem->updateAnalyticalSolution(x,*gridVariables);
// make the new solution the old solution
xOld = x;
......
......@@ -137,7 +137,8 @@ public:
name_ = getParam<std::string>("Problem.Name");
exactWaterDepth_.resize(fvGridGeometry->numDofs(), 0.0);
exactVelocityX_.resize(fvGridGeometry->numDofs(), 0.0);
hBoundary_ = this->gauklerManningStrickler(discharge_,manningN_,bedSlope_);
constManningN_ = this->spatialParams().frictionValue();
hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
}
//! Get the analytical water depth
......@@ -165,8 +166,7 @@ public:
//! Udpate the analytical solution
template<class SolutionVector, class GridVariables>
void updateAnalyticalSolution(const SolutionVector& curSol,
const GridVariables& gridVariables,
const Scalar time)
const GridVariables& gridVariables)
{
using std::pow;
using std::abs;
......@@ -175,11 +175,9 @@ public:
{
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, curSol);
Scalar h = this->gauklerManningStrickler(discharge_,manningN_,bedSlope_);
Scalar h = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
Scalar u = abs(discharge_)/h;
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
......@@ -233,11 +231,11 @@ public:
const auto& elementIndex = scv.elementIndex();
const auto gravity = this->spatialParams().gravity(globalPos);
const auto manningN = this->spatialParams().frictionValue(globalPos);
auto h = volVars.waterDepth();
auto u = volVars.velocity(0);
auto v = volVars.velocity(1);
auto ustarH = FrictionLawManning<Scalar>().computeUstarH(h,manningN_, gravity);
auto ustarH = FrictionLawManning<Scalar>().computeUstarH(h,manningN,gravity);
auto uv = sqrt(pow(u,2.0) + pow(v,2.0));
source[0] = 0.0;
......@@ -291,7 +289,7 @@ public:
const auto& insideVolVars = elemVolVars[insideScv];
const auto& nxy = scvf.unitOuterNormal();
const auto gravity = this->spatialParams().gravity(scvf.center());
std::vector<Scalar> boundaryStateVariables;
std::array<Scalar, 3> boundaryStateVariables;
// impose discharge at the left side
if (scvf.center()[0] < 0.0 + eps_)
......@@ -374,9 +372,9 @@ private:
std::vector<Scalar> exactWaterDepth_;
std::vector<Scalar> exactVelocityX_;
Scalar hBoundary_;
Scalar constManningN_; // analytic solution is only available for const friction.
static constexpr Scalar bedSlope_ = 0.001;
static constexpr Scalar manningN_ = 0.025;
static constexpr Scalar discharge_ = -1.0; // discharge at the inflow boundary
static constexpr Scalar channelLenght_ = 100.0;
static constexpr Scalar bedSurfaceLeft = 10.0;
......
......@@ -68,6 +68,24 @@ public:
return gravity_;
}
/*! \brief Define the friction value.
*
* \return friction value
*/
Scalar frictionValue(const GlobalPosition& globalPos) const
{
return frictionValue_;
}
/*! \brief Define the friction value.
*
* \return friction value
*/
Scalar frictionValue() const
{
return frictionValue_;
}
/*! \brief Define the bed surface
*
* \param element The current element
......@@ -84,7 +102,7 @@ public:
private:
static constexpr Scalar gravity_ = 9.81;
static constexpr Scalar bedSlope_ = 0.001;
static constexpr Scalar frictionValue_ = 0.02;
static constexpr Scalar frictionValue_ = 0.025;
};
} // 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