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

[swe][frictionlaws] Make FrictionLaw a runtime parameter

Also some of the smaller suggestions of the review process are
integrated.
parent a5f86adb
......@@ -36,10 +36,10 @@ namespace ShallowWater {
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar, 3> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
{
std::array<Scalar, 3> cellStateRight;
......@@ -60,10 +60,10 @@ std::array<Scalar, 3> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar, 3> fixedDischargeBoundary(const Scalar qlocal,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
{
std::array<Scalar, 3> cellStateRight;
using std::abs;
......
......@@ -21,8 +21,8 @@
* \ingroup Fluidmatrixinteractions
* \copydoc Dumux::FrictionLaw
*/
#ifndef DUMUX_FRICTIONLAW_HH
#define DUMUX_FRICTIONLAW_HH
#ifndef DUMUX_MATERIAL_FLUIDMATRIX_FRICTIONLAW_HH
#define DUMUX_MATERIAL_FLUIDMATRIX_FRICTIONLAW_HH
#include <algorithm>
#include <cmath>
......@@ -45,7 +45,7 @@ public:
* \return ustar_h friction used for the source term in shallow water models.
*/
virtual Scalar computeUstarH(const Scalar waterDepth) =0;
const virtual Scalar computeUstarH(const Scalar waterDept, const Scalar frictionValue) = 0;
/*!
* \brief Limit the friction for small water depth.
......@@ -55,38 +55,39 @@ public:
* So the friciton term get's not extreme large for small water
* depths.
*
* ------------------------- minUpperh -----------
* ------------------------- minUpperH -----------
*
*
*
* ------------------------rough_h ---------------
* ------------------------roughnessHeight ---------------
* /\ /\ roughness /grain\
* -------------------------------bottom ------------------
* /////////////////////////////////////////////////
*
* For the limiting the LET model is used, which is usually applied in the
* porous media flow to limit the permeability due to the saturation. It employs
* the three empirical paramaters L, E and T, which describe the limiting curve.
* the three empirical paramaters L, E and T, which describe the limiting curve (mobility).
*
* \param rough_h roughness height of the representive structure (e.g. largest grain size).
* auto mobility = (mobility_max * pow(sw,L))/(pow(sw,L) + E * pow(1.0-sw,T));
*
* For the limitation of the roughness height L = 0.0, T = 2.0 and E = 1.0 are choosen.
* Therefore the calculation of the mobility is simplified significantly.
*
* \param roughnessHeight roughness height of the representive structure (e.g. largest grain size).
* \param waterDepth water depth.
*/
Scalar limitRoughH(Scalar rough_h, const Scalar waterDepth)
Scalar limitRoughH(const Scalar roughnessHeight, const Scalar waterDepth)
{
using std::pow;
using std::min;
using std::max;
const Scalar letL = 0.0; //!< empirical parameter of the LET model
const Scalar letT = 2.0; //!< empirical parameter of the LET model
const Scalar letE = 1.0; //!< empirical parameter of the LET model
Scalar mobility_max = 1.0; //!< maximal mobility
auto minUpperH = rough_h * 2.0;
auto sw = min(waterDepth * (1.0/minUpperH),1.0);
Scalar minUpperH = roughnessHeight * 2.0;
Scalar 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);
auto mobility = mobility_max /(1 + (1.0-sw)*(1.0-sw));
return roughnessHeight * (1.0 - mobility);
}
};
......
......@@ -21,8 +21,8 @@
* \ingroup Fluidmatrixinteractions
* \copydoc Dumux::FrictionLawManning
*/
#ifndef DUMUX_FRICTIONLAW_MANNING_HH
#define DUMUX_FRICTIONLAW_MANNING_HH
#ifndef DUMUX_MATERIAL_FLUIDMATRIX_FRICTIONLAW_MANNING_HH
#define DUMUX_MATERIAL_FLUIDMATRIX_FRICTIONLAW_MANNING_HH
#include <algorithm>
#include <cmath>
......@@ -37,33 +37,33 @@ namespace Dumux {
*/
template <typename Scalar>
class FrictionLawManning : FrictionLaw<Scalar>
class FrictionLawManning : public FrictionLaw<Scalar>
{
public:
FrictionLawManning(const Scalar manningN, const Scalar gravity)
: manningN_(manningN), gravity_(gravity) {}
FrictionLawManning(const Scalar gravity)
: gravity_(gravity) {}
/*!
* \brief Compute the friction ustar_h.
*
* \param waterDepth water depth.
* \param frictionValue The Manning friction value.
*
* \return ustar_h friction used for the source term in shallow water models.
*/
Scalar computeUstarH(const Scalar waterDepth)
const Scalar computeUstarH(const Scalar waterDepth, const Scalar frictionValue) final
{
using std::pow;
Scalar ustar_h = 0.0;
Scalar rough_h = pow(25.68/(1.0/manningN_),6.0);
Scalar rough_h = pow(25.68/(1.0/frictionValue),6.0);
rough_h = this->limitRoughH(rough_h, waterDepth);
auto cfric = pow((waterDepth + rough_h),1.0/6.0) * 1.0/(manningN_);
auto cfric = pow((waterDepth + rough_h),1.0/6.0) * 1.0/(frictionValue);
ustar_h = gravity_ / pow(cfric,2.0);
return ustar_h;
}
private:
Scalar manningN_;
Scalar gravity_;
};
......
......@@ -21,8 +21,8 @@
* \ingroup Fluidmatrixinteractions
* \copydoc Dumux::FrictionLawNikuradse
*/
#ifndef DUMUX_FRICTIONLAW_NIKURADSE_HH
#define DUMUX_FRICTIONLAW_NIKURADSE_HH
#ifndef DUMUX_MATERIAL_FLUIDMATRIX_FRICTIONLAW_NIKURADSE_HH
#define DUMUX_MATERIAL_FLUIDMATRIX_FRICTIONLAW_NIKURADSE_HH
#include <algorithm>
#include <cmath>
......@@ -37,33 +37,30 @@ namespace Dumux {
*/
template <typename Scalar>
class FrictionLawNikuradse : FrictionLaw<Scalar>
class FrictionLawNikuradse : public FrictionLaw<Scalar>
{
public:
FrictionLawNikuradse(const Scalar ks)
: ks_(ks) {}
/*!
* \brief Compute the friction ustar_h.
*
* \param waterDepth water depth.
* \param frictionValue The equivalent sand roughness.
*
* \return ustar_h friction used for the source term in shallow water models.
*/
Scalar computeUstarH(const Scalar waterDepth)
const Scalar computeUstarH(const Scalar waterDepth, const Scalar frictionValue) final
{
using std::pow;
using std::log;
Scalar ustar_h = 0.0;
Scalar rough_h = ks_;
Scalar rough_h = frictionValue;
rough_h = this->limitRoughH(rough_h, waterDepth);
ustar_h = pow(0.41,2.0)/pow(log((12*(waterDepth + rough_h))/ks_),2.0);
ustar_h = pow(0.41,2.0)/pow(log((12*(waterDepth + rough_h))/frictionValue),2.0);
return ustar_h;
}
private:
Scalar ks_;
};
} // end namespace Dumux
......
......@@ -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")
[Problem]
Name = damBreak
Name = dambreak
[TimeLoop]
TEnd = 1.0 # [s]
......
......@@ -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")
......@@ -23,7 +23,6 @@
*/
#include <config.h>
#include <ctime>
#include <iostream>
......@@ -32,7 +31,6 @@
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
......
[Problem]
Name = roughChannel
Name = roughchannel
FrictionValue = 0.025 # [-]
BedSlope = 0.001 # [-]
Gravity = 9.81 # [m/s^2]
Discharge = -1.0 # [m/s] discharge at the inflow boundary
FrictionLaw = Manning
[TimeLoop]
TEnd = 3600.0 # [s]
TEnd = 900.0 # [s]
MaxTimeStepSize = 60.0 # [s]
DtInitial = 1.0 # [s]
[Grid]
Positions0 = 0.0 1000.0
Positions1 = 0.0 10.0
Cells0 = 1000
Cells1 = 10
Positions0 = 0.0 500.0
Positions1 = 0.0 5.0
Cells0 = 500
Cells1 = 5
[Newton]
EnablePartialReassembly = true
......@@ -27,11 +27,14 @@
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include "spatialparams.hh"
#include <dumux/common/parameters.hh>
#include <dumux/freeflow/shallowwater/model.hh>
#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh>
namespace Dumux {
......@@ -106,7 +109,7 @@ struct EnableGridFluxVariablesCache<TypeTag, TTag::RoughChannel>
* h = \left(\frac{n*q}{\sqrt{I_s}} \right)^{3/5}
* \f]
*
* The formular of Gaukler Manning and Strickler is also used to calculate the analytic solution.
* The formula of Gaukler Manning and Strickler is also used to calculate the analytic solution.
*
* This problem uses the \ref ShallowWaterModel
*/
......@@ -137,7 +140,20 @@ public:
exactWaterDepth_.resize(fvGridGeometry->numDofs(), 0.0);
exactVelocityX_.resize(fvGridGeometry->numDofs(), 0.0);
constManningN_ = this->spatialParams().frictionValue();
bedSlope_ = getParam<Scalar>("Problem.BedSlope");
discharge_ = getParam<Scalar>("Problem.Discharge");
hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
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);
}
else
{
std::cout<<"Change the FrictionLaw in params.input. This test is only valid for Manning!";
}
}
//! Get the analytical water depth
......@@ -159,7 +175,7 @@ public:
using std::abs;
using std::sqrt;
return std::pow(std::abs(discharge)*manningN/sqrt(bedSlope), 0.6);
return std::pow(std::abs(discharge)*manningN/sqrt(bedSlope), 0.6);
}
//! Udpate the analytical solution
......@@ -167,7 +183,6 @@ public:
void updateAnalyticalSolution(const SolutionVector& curSol,
const GridVariables& gridVariables)
{
using std::pow;
using std::abs;
//compute solution for all elements
for (const auto& element : elements(this->fvGridGeometry().gridView()))
......@@ -219,22 +234,23 @@ public:
* 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,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
using std::hypot;
NumEqVector source(0.0);
const auto& globalPos = scv.center();
const auto& volVars = elemVolVars[scv];
const Scalar manningN = this->spatialParams().frictionValue(globalPos);
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>(manningN,gravity).computeUstarH(h);
auto uv = sqrt(pow(u,2.0) + pow(v,2.0));
const Scalar h = volVars.waterDepth();
const Scalar u = volVars.velocity(0);
const Scalar v = volVars.velocity(1);
const Scalar ustarH = frictionLaw_->computeUstarH(h, manningN);
const Scalar uv = hypot(u,v);
source[0] = 0.0;
source[1] = -ustarH * u * uv;
......@@ -307,7 +323,6 @@ public:
insideVolVars.velocity(1),
nxy);
}
// no flow boundary
else
{
......@@ -352,32 +367,25 @@ public:
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
values[0] = hBoundary_;
values[1] = 0.0;
values[1] = abs(discharge_)/hBoundary_;
values[2] = 0.0;
return values;
};
// \}
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 discharge_ = -1.0; // discharge at the inflow boundary
static constexpr Scalar channelLenght_ = 100.0;
static constexpr Scalar bedSurfaceLeft = 10.0;
Scalar bedSlope_;
Scalar discharge_; // discharge at the inflow boundary
std::shared_ptr<FrictionLaw<Scalar>> frictionLaw_;
static constexpr Scalar eps_ = 1.0e-6;
std::string name_;
};
......
......@@ -25,6 +25,7 @@
#define DUMUX_ROUGH_CHANNEL_SPATIAL_PARAMETERS_HH
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/common/parameters.hh>
namespace Dumux {
......@@ -49,15 +50,11 @@ class RoughChannelSpatialParams
public:
RoughChannelSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{}
/*! \brief Define the porosity in [-].
*
* \param globalPos The global position where we evaluate
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 1.0; }
{
gravity_ = getParam<Scalar>("Problem.Gravity");
bedSlope_ = getParam<Scalar>("Problem.BedSlope");
frictionValue_ = getParam<Scalar>("Problem.FrictionValue");
}
/*! \brief Define the gravitation.
*
......@@ -68,6 +65,15 @@ public:
return gravity_;
}
/*! \brief Define the gravitation.
*
* \return gravity constant
*/
Scalar gravity() const
{
return gravity_;
}
/*! \brief Define the friction value.
*
* \return friction value
......@@ -100,9 +106,9 @@ public:
}
private:
static constexpr Scalar gravity_ = 9.81;
static constexpr Scalar bedSlope_ = 0.001;
static constexpr Scalar frictionValue_ = 0.025;
Scalar gravity_;
Scalar bedSlope_;
Scalar frictionValue_;
};
} // 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