Commit 04d8f17f authored by Leopold Stadler's avatar Leopold Stadler Committed by Timo Koch
Browse files

[test][dambreak] add analytical solution and test

parent b1e38a08
......@@ -30,6 +30,15 @@
namespace Dumux {
namespace ShallowWater{
template<typename Scalar>
struct RiemannSolution {
std::array<Scalar,3> flux;
Scalar waterDepth;
Scalar velocityX;
Scalar velocityY;
};
/*!
* \ingroup ShallowWater
* \brief Exact Riemann solver for Shallow water equations.
......@@ -46,15 +55,17 @@ namespace ShallowWater{
* \param vl velocityY on the left side
* \param vr velocityY on the right side
* \param grav gravity constant
* \param s sample point (default = 0 since x = 0 for flux computation)
*/
template<class Scalar>
std::array<Scalar,3> exactRiemann(const Scalar& dl,
const Scalar& dr,
const Scalar& ul,
const Scalar& ur,
const Scalar& vl,
const Scalar& vr,
const Scalar& grav)
RiemannSolution<Scalar> exactRiemann(const Scalar& dl,
const Scalar& dr,
const Scalar& ul,
const Scalar& ur,
const Scalar& vl,
const Scalar& vr,
const Scalar& grav,
const Scalar& s = 0.0)
{
Scalar fl = 0.0; //function fl
......@@ -66,7 +77,6 @@ std::array<Scalar,3> exactRiemann(const Scalar& dl,
Scalar cr = 0.0; //celerity sqrt(g*h) right side
Scalar c = 0.0; //celerity
Scalar cs = 0.0; //celerity
Scalar s = 0.0; //sample point
Scalar dcrit = 0.0; //critical water depth
Scalar dmin = 0.0; //minimum depth
Scalar gel = 0.0; //gravity component
......@@ -102,9 +112,6 @@ std::array<Scalar,3> exactRiemann(const Scalar& dl,
cr = sqrt(grav * max(dr,0.0));
dcrit = (ur - ul) - 2.0*(cl+cr);
//we use s as 0.0 since we are at the x = 0
s = 0.0;
// dry case
if ((dl <= 0)||(dr <= 0)||(dcrit >= 0))
{
......@@ -448,13 +455,16 @@ std::array<Scalar,3> exactRiemann(const Scalar& dl,
}
}
//============================== Flux computation ===================
std::array<Scalar,3> flux;
RiemannSolution<Scalar> result;
flux[0] = d*u;
flux[1] = d * u * u + 0.5 * grav * d * d;
flux[2] = d*u*v;
result.flux[0] = d * u;
result.flux[1] = d * u * u + 0.5 * grav * d * d;
result.flux[2] = d * u * v;
result.waterDepth = d;
result.velocityX = u;
result.velocityY = v;
return flux;
return result;
}
} // end namespace ShallowWater
......
......@@ -73,10 +73,8 @@ std::array<Scalar,3> riemannProblem(Scalar waterDepthLeft,
Scalar gravity,
GlobalPosition nxy)
{
std::array<Scalar,3> riemannFlux;
using std::max;
//Hydrostatic reconstrucion after Audusse
Scalar dzl = max(0.0,bedSurfaceRight - bedSurfaceLeft);
Scalar waterDepthLeftReconstructed = max(0.0, waterDepthLeft - dzl);
......@@ -98,18 +96,18 @@ std::array<Scalar,3> riemannProblem(Scalar waterDepthLeft,
velocityXRight = nxy[0] * tempFlux + nxy[1] * velocityYRight;
velocityYRight = -nxy[1] * tempFlux + nxy[0] * velocityYRight;
riemannFlux = ShallowWater::exactRiemann(waterDepthLeftReconstructed,
waterDepthRightReconstructed,
velocityXLeft,
velocityXRight,
velocityYLeft,
velocityYRight,
gravity);
ShallowWater::RiemannSolution<Scalar> riemannResult = ShallowWater::exactRiemann(waterDepthLeftReconstructed,
waterDepthRightReconstructed,
velocityXLeft,
velocityXRight,
velocityYLeft,
velocityYRight,
gravity);
//redo rotation
tempFlux = riemannFlux[1];
riemannFlux[1] = nxy[0] * tempFlux - nxy[1] * riemannFlux[2];
riemannFlux[2] = nxy[1] * tempFlux + nxy[0] * riemannFlux[2];
tempFlux = riemannResult.flux[1];
riemannResult.flux[1] = nxy[0] * tempFlux - nxy[1] * riemannResult.flux[2];
riemannResult.flux[2] = nxy[1] * tempFlux + nxy[0] * riemannResult.flux[2];
//Add reconstruction flux from Audusse reconstruction
Scalar hgzl = 0.5 * (waterDepthLeftReconstructed + waterDepthLeft) * (waterDepthLeftReconstructed - waterDepthLeft);
......@@ -125,9 +123,9 @@ std::array<Scalar,3> riemannProblem(Scalar waterDepthLeft,
*/
std::array<Scalar,3> localFlux{0.0,0.0,0.0};
localFlux[0] = riemannFlux[0] * mobility;
localFlux[1] = riemannFlux[1] - hdxzl;
localFlux[2] = riemannFlux[2] - hdyzl;
localFlux[0] = riemannResult.flux[0] * mobility;
localFlux[1] = riemannResult.flux[1] - hdxzl;
localFlux[2] = riemannResult.flux[2] - hdyzl;
return localFlux;
}
......
dune_symlink_to_source_files(FILES "params.input" references)
dune_symlink_to_source_files(FILES "params.input")
dune_add_test(NAME test_shallowwater_dambreak
SOURCES main.cc
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/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
--command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_dambreak params.input
-Problem.Name dambreak")
......@@ -136,6 +136,9 @@ int main(int argc, char** argv) try
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);
......@@ -165,6 +168,9 @@ int main(int argc, char** argv) try
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();
......
......@@ -2,9 +2,9 @@
Name = damBreak
[TimeLoop]
TEnd = 1.5 # [s]
MaxTimeStepSize = 0.01 # [s]
DtInitial = 0.01 # [s]
TEnd = 1.0 # [s]
MaxTimeStepSize = 0.005 # [s]
DtInitial = 0.005 # [s]
[Grid]
Positions0 = 0.0 20.0
......
......@@ -31,6 +31,7 @@
#include <dumux/freeflow/shallowwater/model.hh>
#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/flux/shallowwater/riemannproblem.hh>
#include <dumux/flux/shallowwater/exactriemann.hh>
namespace Dumux
......@@ -120,7 +121,8 @@ class DamBreakProblem : public ShallowWaterProblem<TypeTag>
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
enum {
// copy some indices for convenience
......@@ -140,6 +142,58 @@ public:
: ParentType(fvGridGeometry)
{
name_ = getParam<std::string>("Problem.Name");
exactWaterDepth_.resize(fvGridGeometry->numDofs(), 0.0);
exactVelocityX_.resize(fvGridGeometry->numDofs(), 0.0);
}
//! Get the analytical water depth
const std::vector<Scalar>& getExactWaterDepth()
{
return exactWaterDepth_;
}
//! Get the analytical velocity
const std::vector<Scalar>& getExactVelocityX()
{
return exactVelocityX_;
}
//! Udpate the analytical solution
void updateAnalyticalSolution(const SolutionVector& curSol,
const GridVariables& gridVariables,
const Scalar time)
{
//compute solution for all elements
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, curSol);
auto elemId = this->fvGridGeometry().elementMapper().index(element);
const auto& globalPos = element.geometry().center();
//compute the position s and the initial water depth at the gate, velocities are zero
Scalar s = (globalPos[0] - gatePosition_)/time;
Scalar waterDepthLeft = initialWaterDepthLeft_;
Scalar waterDepthRight = initialWaterDepthRight_;
ShallowWater::RiemannSolution<Scalar> riemannResult = ShallowWater::exactRiemann(waterDepthLeft,
waterDepthRight,
0.0,
0.0,
0.0,
0.0,
9.81,
s);
exactWaterDepth_[elemId] = riemannResult.waterDepth;
exactVelocityX_[elemId] = riemannResult.velocityX;
}
}
/*!
......@@ -183,7 +237,6 @@ public:
* \param elemVolVars
* \param scvf
*/
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
......@@ -256,20 +309,20 @@ public:
PrimaryVariables values(0.0);
values[0] = 1.0;
values[0] = initialWaterDepthRight_;
values[1] = 0.0;
values[2] = 0.0;
// water level on the left side of the gate
if (globalPos[0] < 10.0 + eps_)
{
values[0] = 4.0;
values[0] = initialWaterDepthLeft_;
}
//water level on the right side of the gate
else
{
values[0] = 1.0;
values[0] = initialWaterDepthRight_;
}
return values;
......@@ -280,6 +333,13 @@ public:
private:
std::vector<Scalar> exactWaterDepth_;
std::vector<Scalar> exactVelocityX_;
static constexpr Scalar initialWaterDepthLeft_ = 4.0;
static constexpr Scalar initialWaterDepthRight_ = 1.0;
static constexpr Scalar channelLenght_ = 20.0;
static constexpr Scalar gatePosition_ = 10.0;
static constexpr Scalar eps_ = 1.0e-6;
std::string name_;
......
This diff is collapsed.
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