Commit 0468997f authored by Martin Utz's avatar Martin Utz
Browse files

[swe][frictionlaws] Add a friction law test

The test provides a simple channel flow based on the friction law after
Manning. Normal flow is assumed, therefore the formular after Gaukler,
Manning and Strickler can be used to get an analytic solution.
parent 3a6fefe8
......@@ -25,7 +25,7 @@
#ifndef DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#define DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#include <array>
#include <vector>
#include <cmath>
namespace Dumux {
......@@ -35,18 +35,14 @@ namespace ShallowWater {
* \brief compute the cell state for fixed water depth boundary.
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar, 3> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
const Scalar waterDepthLeft,
const Scalar waterDepthRight,
const Scalar velocityXLeft,
const Scalar velocityXRight,
const Scalar velocityYLeft,
const Scalar velocityYRight,
const Scalar gravity,
const GlobalPosition& nxy)
std::vector<Scalar> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
{
std::array<Scalar, 3> cellStateRight;
std::vector<Scalar> cellStateRight(3);
cellStateRight[0] = waterDepthBoundary;
using std::sqrt;
......@@ -63,25 +59,20 @@ std::array<Scalar, 3> fixedWaterDepthBoundary(const Scalar waterDepthBoundary,
* \brief compute the cell state for a fixed discharge boundary.
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar, 3> fixedDischargeBoundary(const Scalar dischargeBoundary,
const Scalar waterDepthLeft,
const Scalar waterDepthRight,
const Scalar velocityXLeft,
const Scalar velocityXRight,
const Scalar velocityYLeft,
const Scalar velocityYRight,
const Scalar gravity,
const GlobalPosition& nxy,
const Scalar faceVolume)
std::vector<Scalar> fixedDischargeBoundary(const Scalar qlocal,
const Scalar waterDepthLeft,
const Scalar velocityXLeft,
const Scalar velocityYLeft,
const GlobalPosition& nxy)
{
std::array<Scalar, 3> cellStateRight;
std::vector<Scalar> cellStateRight(3);
using std::abs;
using std::sqrt;
using std::max;
// only impose if abs(q) > 0
if (abs(dischargeBoundary) > 1.0e-9)
if (abs(qlocal) > 1.0e-9)
{
const auto qlocal = dischargeBoundary/faceVolume;
const auto uboundIn = nxy[0]*velocityXLeft + nxy[1]*velocityYLeft;
const auto alphal = uboundIn + 2.0*sqrt(9.81 * waterDepthLeft);
......
......@@ -28,7 +28,7 @@
#include <algorithm>
#include <cmath>
#include"nikuradse.hh"
#include "nikuradse.hh"
namespace Dumux {
/*!
......@@ -39,27 +39,27 @@ namespace Dumux {
*/
template <typename Scalar>
class FrictionLawManning : FrictionLawNikuradse
class FrictionLawManning : FrictionLawNikuradse<Scalar>
{
public:
/*!
* \brief Compute the friction ustar_h.
*
* \param h water depth.
* \param ks the Strickler friction value.
* \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 ks, const Scalar gravity)
Scalar computeUstarH(const Scalar h,const Scalar manningN, const Scalar gravity)
{
using std::pow;
using std::log;
Scalar ustar_h = 0.0;
Scalar rough_h = pow(25.68/(1.0/ks),6.0);
Scalar rough_h = pow(25.68/(1.0/manningN),6.0);
rough_h = limitRoughH(rough_h, h);
rough_h = this->limitRoughH(rough_h, h);
auto cfric = pow((h + rough_h),1.0/6.0) * 1.0/(ks);
auto cfric = pow((h + rough_h),1.0/6.0) * 1.0/(manningN);
ustar_h = gravity / pow(cfric,2.0);
return ustar_h;
}
......
......@@ -46,6 +46,7 @@ public:
*
* \param h 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)
......@@ -67,7 +68,7 @@ public:
* We define a water depth minUpperH. If the water depth is
* smaller, we start to limit the friciton.
* So the friciton term get's not extreme large for small water
* depths
* depths.
*
* ------------------------- minUpperh -----------
*
......@@ -79,11 +80,11 @@ public:
* /////////////////////////////////////////////////
*
* For the limiting the LET model is used, which is usually applied in the
* porouse media flow to limit the permeability due to the saturation. It employs
* 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.
*
* \param rough_h roughness height of the representive structure (e.g. largest grain size).
* \params h water depth.
* \param h water depth.
*/
Scalar limitRoughH(Scalar rough_h, const Scalar h)
{
......@@ -91,18 +92,17 @@ public:
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
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(h * (1.0/minUpperH),1.0);
sw = max(0.0,sw);
auto mobility = (krw * pow(sw,letL))/(pow(sw,letL) + letE * pow(1.0-sw,letT));
auto mobility = (mobility_max * pow(sw,letL))/(pow(sw,letL) + letE * pow(1.0-sw,letT));
return rough_h * (1.0 - mobility);
}
};
} // end namespace Dumux
......
add_subdirectory(dambreak)
add_subdirectory(roughchannel)
......@@ -238,6 +238,19 @@ 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),
......
dune_symlink_to_source_files(FILES "params.input")
dune_add_test(NAME test_shallowwater_roughchannel
SOURCES main.cc
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
--zeroThreshold {"velocityY":1e-14}
--command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_roughchannel params.input
-Problem.Name roughchannel")
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup ShallowWaterTests
* \brief A test for the shallow water model (rough channel).
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#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>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include "problem.hh"
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::RoughChannel;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(fvGridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x);
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
problem->updateAnalyticalSolution(x,*gridVariables,0.0);
IOFields::initOutputModule(vtkWriter);
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
// the linear solver
using LinearSolver = Dumux::AMGBackend<TypeTag>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
//! set some check point at the end of the time loop
timeLoop->setCheckPoint(tEnd);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
nonLinearSolver.solve(x,*timeLoop);
// update the analytical solution
problem->updateAnalyticalSolution(x,*gridVariables,timeLoop->time()+timeLoop->timeStepSize());
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
if (timeLoop->isCheckPoint())
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by newton controller
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
}
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (const Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
}
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
[Problem]
Name = roughChannel
[TimeLoop]
TEnd = 3600.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
[Newton]
EnablePartialReassembly = true
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup ShallowWaterTests
* \brief A test for the Shallow water model (rough channel).
*/
#ifndef DUMUX_ROUGH_CHANNEL_TEST_PROBLEM_HH
#define DUMUX_ROUGH_CHANNEL_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include "spatialparams.hh"
#include <dumux/freeflow/shallowwater/model.hh>
#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
namespace Dumux {
template <class TypeTag>
class RoughChannelProblem;
// Specify the properties for the problem
namespace Properties {
// Create new type tags
namespace TTag {
struct RoughChannel { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
} // end namespace TTag
template<class TypeTag>
struct Grid<TypeTag, TTag::RoughChannel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::RoughChannel>
{ using type = Dumux::RoughChannelProblem<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::RoughChannel>
{
private:
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = RoughChannelSpatialParams<FVGridGeometry, Scalar>;
};
template<class TypeTag>
struct EnableFVGridGeometryCache<TypeTag, TTag::RoughChannel>
{ static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::RoughChannel>
{ static constexpr bool value = false; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::RoughChannel>
{ static constexpr bool value = false; };
} // end namespace Properties
/*!
* \ingroup ShallowWaterTests
* \brief A simple flow in a rough channel with friction law after Manning.
*
* The domain is 1000 meters long and 10 meters wide. At the left border a discharge
* boundary condition is applied and at the right border a water depth boundary condition.
* All other boundaries are set to no-flow. Normal flow is assumed, therefor the water depth
* at the right border can be calculated with the formular of Gaukler-Manning-Strickler.
*
* \f[
* v_m = 1/n * R_{hy}^{2/3} * I_s^{1/2}
* \f]
*
* With the mean velocity
* \f[
* v_m = \frac{q}/{h}
* \f]
* the friction value n after Manning
* the hydraulic radius R_{hy} equal to the water depth h (because normal flow is assumed)
* the bed slope I_s and the unity inflow discharge q.
*
* Therefore h can be calculated with
*
* \f[
* 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.
*
* This problem uses the \ref ShallowWaterModel
*/
template <class TypeTag>
class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
{
using ParentType = ShallowWaterProblem<TypeTag>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
public:
RoughChannelProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
name_ = getParam<std::string>("Problem.Name");
exactWaterDepth_.resize(fvGridGeometry->numDofs(), 0.0);
exactVelocityX_.resize(fvGridGeometry->numDofs(), 0.0);
hBoundary_ = this->gauklerManningStrickler(discharge_,manningN_,bedSlope_);
}
//! Get the analytical water depth
const std::vector<Scalar>& getExactWaterDepth()
{
return exactWaterDepth_;
}
//! Get the analytical velocity
const std::vector<Scalar>& getExactVelocityX()
{
return exactVelocityX_;
}
//! Calculate the water depth with Gaukler-Manning-Strickler
Scalar gauklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope)
{
using std::pow;
using std::abs;
using std::sqrt;
return std::pow(std::abs(discharge)*manningN/sqrt(bedSlope), 0.6);
}
//! Udpate the analytical solution
template<class SolutionVector, class GridVariables>
void updateAnalyticalSolution(const SolutionVector& curSol,