Commit f9efc5cd authored by Timo Koch's avatar Timo Koch
Browse files

[cleanup][shallowwater] Cleanup test and model

parent 46a34d11
......@@ -54,6 +54,7 @@ public:
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const auto& nxy = scvf.unitOuterNormal();
const auto gravity = problem.spatialParams().gravity(scvf.center());
auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
outsideVolVars.waterDepth(),
......@@ -63,7 +64,7 @@ public:
outsideVolVars.velocity(1),
insideVolVars.bedSurface(),
outsideVolVars.bedSurface(),
insideVolVars.gravity(),
gravity,
nxy);
NumEqVector localFlux(0.0);
......
......@@ -26,8 +26,7 @@
#define DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#include <array>
#include <algorithm>
#include <dumux/common/math.hh>
#include <cmath>
namespace Dumux {
namespace ShallowWater {
......@@ -36,23 +35,24 @@ namespace ShallowWater {
* \brief compute the cell state for fixed water depth boundary.
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar,3> fixedWaterDepthBoundary(Scalar waterDepthBoundary,
Scalar waterDepthLeft,
Scalar waterDepthRight,
Scalar velocityXLeft,
Scalar velocityXRight,
Scalar velocityYLeft,
Scalar velocityYRight,
Scalar gravity,
GlobalPosition nxy)
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::array<Scalar,3> cellStateRight;
std::array<Scalar, 3> cellStateRight;
cellStateRight[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]);
auto uboundIn = nxy[0] * velocityXLeft + nxy[1] * velocityYLeft ;
auto uboundQut = uboundIn + 2.0 * sqrt(9.81 * waterDepthLeft) - 2.0 * sqrt(9.81 * cellStateRight[0]);
cellStateRight[1] = (nxy[0] * uboundQut); // we only use the normal part
cellStateRight[2] = (nxy[1] * uboundQut); // we only use the normal part
......@@ -63,45 +63,46 @@ std::array<Scalar,3> fixedWaterDepthBoundary(Scalar waterDepthBoundary,
* \brief compute the cell state for a fixed discharge boundary.
*/
template<class Scalar, class GlobalPosition>
std::array<Scalar,3> fixedDischargeBoundary(Scalar dischargeBoundary,
Scalar waterDepthLeft,
Scalar waterDepthRight,
Scalar velocityXLeft,
Scalar velocityXRight,
Scalar velocityYLeft,
Scalar velocityYRight,
Scalar gravity,
GlobalPosition nxy,
Scalar faceVolume)
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::array<Scalar,3> cellStateRight;
using std::pow;
std::array<Scalar, 3> cellStateRight;
using std::abs;
using std::sqrt;
//olny impose if abs(q) > 0
if (abs(dischargeBoundary) > 1.0E-9){
auto qlocal = (dischargeBoundary) /faceVolume;
auto uboundIn = nxy[0] * velocityXLeft + nxy[1] * velocityYLeft;
auto alphal = uboundIn + 2.0 * sqrt(9.81 * waterDepthLeft);
// only impose if abs(q) > 0
if (abs(dischargeBoundary) > 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);
//initial guess for hstar solved with newton
Scalar hstar = 0.1;
Scalar tol_hstar = 1.0E-12;
Scalar ink_hstar = 1.0E-9;
int maxstep_hstar = 30;
constexpr Scalar tol_hstar = 1.0E-12;
constexpr Scalar ink_hstar = 1.0E-9;
constexpr int maxstep_hstar = 30;
for(int i = 0; i < maxstep_hstar; ++i){
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 dx_hstar = -f_hstar/df_hstar;
hstar = max(hstar - dx_hstar,0.001);
if (pow(dx_hstar,2.0) < tol_hstar){
if (dx_hstar*dx_hstar < tol_hstar)
break;
}
}
auto qinner = (nxy[0] * waterDepthLeft * velocityYLeft) - (nxy[1] * waterDepthLeft * velocityXLeft);
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;
......
......@@ -27,7 +27,7 @@
#include <dumux/common/properties.hh>
#include <dumux/flux/fluxvariablesbase.hh>
namespace Dumux{
namespace Dumux {
/*!
* \ingroup ShallowWaterModel
......@@ -36,13 +36,11 @@ namespace Dumux{
*/
template<class TypeTag>
class ShallowWaterFluxVariables
: public FluxVariablesBase<GetPropType<TypeTag,
Properties::Problem>,
: public FluxVariablesBase<GetPropType<TypeTag, Properties::Problem>,
typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView,
typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView,
typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView>
{
//using ParentType = FluxVariablesBase<TypeTag>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
......@@ -62,7 +60,6 @@ class ShallowWaterFluxVariables
static constexpr bool enableAdvection = ModelTraits::enableAdvection();
static constexpr bool enableDiffusion = ModelTraits::enableDiffusion();
public:
/*!
......@@ -75,15 +72,10 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) const
{
NumEqVector fluxVector(0.0);
if (enableAdvection)
{
return AdvectionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
}
else
{
return fluxVector;
}
return NumEqVector(0.0);
}
/*!
......@@ -96,17 +88,11 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) const
{
NumEqVector fluxVector(0.0);
// TODO: add diffusive flux (e.g. tracer and viscosity)
if (enableDiffusion)
{
// TODO: add diffusive flux (e.g. tracer and viscosity)
return fluxVector;
}
return NumEqVector(0.0);
else
{
return fluxVector;
}
return NumEqVector(0.0);
}
};
......
......@@ -24,9 +24,7 @@
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_INDICES_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_INDICES_HH
#include <dumux/common/properties.hh>
namespace Dumux{
namespace Dumux {
// \{
/*!
......@@ -37,7 +35,6 @@ namespace Dumux{
*/
struct ShallowWaterIndices
{
static constexpr int dimXIdx = 0; //!< Index of the x-component of a vector of size dim
static constexpr int dimYIdx = 1; //!< Index of the y-component of a vector of size dim
......@@ -49,7 +46,6 @@ struct ShallowWaterIndices
static constexpr int velocityXIdx = momentumXBalanceIdx; //!< Index of the velocity in a solution vector
static constexpr int velocityYIdx = momentumYBalanceIdx; //!< Index of the velocity in a solution vector
static constexpr int velocityOffset = velocityXIdx; //!< Offset vor the velocity index
};
// \}
......
......@@ -27,7 +27,7 @@
#include <dumux/io/name.hh>
#include <string>
namespace Dumux{
namespace Dumux {
/*!
* \ingroup ShallowWaterModel
......
......@@ -25,8 +25,6 @@
#define DUMUX_FREEFLOW_SHALLOW_WATER_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
#include <dune/common/hybridutilities.hh>
namespace Dumux{
......@@ -35,10 +33,10 @@ namespace Dumux{
* \brief Element-wise calculation of the residual for the shallow water equations
*/
template<class TypeTag>
class ShallowWaterResidual : public GetPropType<TypeTag, Properties::BaseLocalResidual>
class ShallowWaterResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
......@@ -50,7 +48,6 @@ class ShallowWaterResidual : public GetPropType<TypeTag, Properties::BaseLocalRe
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
public:
......@@ -77,7 +74,6 @@ public:
storage[Indices::massBalanceIdx] = volVars.waterDepth();
storage[Indices::momentumXBalanceIdx] = volVars.waterDepth() * volVars.velocity(0);
storage[Indices::momentumYBalanceIdx] = volVars.waterDepth() * volVars.velocity(1);
return storage;
}
......@@ -99,10 +95,8 @@ public:
{
NumEqVector flux(0.0);
FluxVariables fluxVars;
auto advectiveFlux = fluxVars.advectiveFlux(problem, element, fvGeometry, elemVolVars, scvf);
auto diffusiveFlux = fluxVars.diffusiveFlux(problem, element, fvGeometry, elemVolVars, scvf);
flux = advectiveFlux + diffusiveFlux;
flux += fluxVars.advectiveFlux(problem, element, fvGeometry, elemVolVars, scvf);
flux += fluxVars.diffusiveFlux(problem, element, fvGeometry, elemVolVars, scvf);
return flux;
}
};
......
......@@ -93,7 +93,6 @@ struct ShallowWaterModelTraits
//! Enable diffusion
static constexpr bool enableDiffusion() { return false; }
};
/*!
......@@ -111,9 +110,7 @@ struct ShallowWaterVolumeVariablesTraits
using ModelTraits = MT;
};
///////////////////////////////////////////////////////////////////////////
// properties for the shallow water model
///////////////////////////////////////////////////////////////////////////
namespace Properties {
//! Type tag for shallow water equation model inherits from model properties
......@@ -121,45 +118,21 @@ namespace TTag {
struct ShallowWater { using InheritsFrom = std::tuple<ModelProperties>; };
}// end namespace TTag
//! The grid flux variables cache with an empty ShallowWaterFluxVariablesCacheFiller
template<class TypeTag>
struct GridFluxVariablesCache<TypeTag, TTag::ShallowWater>
{
private:
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FluxVariablesCache = GetPropType<TypeTag, Properties::FluxVariablesCache>;
using FluxVariablesCacheFiller = GetPropType<TypeTag, Properties::FluxVariablesCacheFiller>;
public:
using type = CCTpfaGridFluxVariablesCache<Problem, FluxVariablesCache, FluxVariablesCacheFiller, enableCache>;
};
//////////////////////////////////////////////////////////////////
// Make new type tags
// Define properties
//////////////////////////////////////////////////////////////////
template<class TypeTag>
struct SolutionDependentAdvection<TypeTag, TTag::ShallowWater> { static constexpr bool value = true; };
template<class TypeTag>
struct SolutionDependentMolecularDiffusion<TypeTag, TTag::ShallowWater> { static constexpr bool value = false; };
template<class TypeTag>
struct SolutionDependentHeatConduction<TypeTag, TTag::ShallowWater> { static constexpr bool value = false; };
template<class TypeTag>
struct EnableThermalNonEquilibrium<TypeTag, TTag::ShallowWater> { static constexpr bool value = false; };
//////////////////////////////////////////////////////////////////
// Type tags
//////////////////////////////////////////////////////////////////
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::ShallowWater> {using type = ShallowWaterModelTraits;};
struct ModelTraits<TypeTag, TTag::ShallowWater>
{ using type = ShallowWaterModelTraits; };
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::ShallowWater> {using type = ShallowWaterResidual<TypeTag>;};
struct LocalResidual<TypeTag, TTag::ShallowWater>
{ using type = ShallowWaterResidual<TypeTag>; };
template<class TypeTag>
struct FluxVariables<TypeTag, TTag::ShallowWater> {using type = ShallowWaterFluxVariables<TypeTag>;};
struct FluxVariables<TypeTag, TTag::ShallowWater>
{ using type = ShallowWaterFluxVariables<TypeTag>; };
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::ShallowWater>
......@@ -167,24 +140,22 @@ struct VolumeVariables<TypeTag, TTag::ShallowWater>
private:
using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
using MT = GetPropType<TypeTag, Properties::ModelTraits>;
using Traits = ShallowWaterVolumeVariablesTraits<PV, MT>;
public:
using type = ShallowWaterVolumeVariables<Traits>;
};
template<class TypeTag>
struct FluxVariablesCache<TypeTag, TTag::ShallowWater>
{
using type = FluxVariablesCaching::EmptyCache< GetPropType<TypeTag, Properties::Scalar> >;
};
{ using type = FluxVariablesCaching::EmptyCache< GetPropType<TypeTag, Properties::Scalar> >; };
template<class TypeTag>
struct FluxVariablesCacheFiller<TypeTag, TTag::ShallowWater> { using type = FluxVariablesCaching::EmptyCacheFiller; };
struct FluxVariablesCacheFiller<TypeTag, TTag::ShallowWater>
{ using type = FluxVariablesCaching::EmptyCacheFiller; };
template<class TypeTag>
struct IOFields<TypeTag, TTag::ShallowWater> {using type = ShallowWaterIOFields;};
struct IOFields<TypeTag, TTag::ShallowWater>
{ using type = ShallowWaterIOFields; };
template<class TypeTag>
struct AdvectionType<TypeTag, TTag::ShallowWater>
......
......@@ -28,7 +28,7 @@
#include <dumux/common/properties.hh>
#include "model.hh"
namespace Dumux{
namespace Dumux {
/*!
* \ingroup ShallowWaterModel
......@@ -39,11 +39,6 @@ class ShallowWaterProblem : public FVProblem<TypeTag>
{
using ParentType = FVProblem<TypeTag>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
......@@ -53,10 +48,12 @@ public:
*
* \param fvGridGeometry The finite volume grid geometry
* \param spatialParams The spatial parameter class
* \param paramGroup The parameter group in which to look for runtime parameters first (default is "")
*/
ShallowWaterProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<SpatialParams> spatialParams)
: ParentType(fvGridGeometry)
std::shared_ptr<SpatialParams> spatialParams,
const std::string& paramGroup = "")
: ParentType(fvGridGeometry, paramGroup)
, spatialParams_(spatialParams)
{}
......@@ -66,9 +63,11 @@ public:
* \param fvGridGeometry The finite volume grid geometry
* \param paramGroup The parameter group in which to look for runtime parameters first (default is "")
*/
ShallowWaterProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
ShallowWaterProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
const std::string& paramGroup = "")
: ShallowWaterProblem(fvGridGeometry,
std::make_shared<SpatialParams>(fvGridGeometry))
std::make_shared<SpatialParams>(fvGridGeometry),
paramGroup)
{}
......@@ -82,11 +81,8 @@ public:
// \}
protected:
GlobalPosition gravity_;
// material properties
std::shared_ptr<SpatialParams> spatialParams_;
private:
std::shared_ptr<SpatialParams> spatialParams_; //!< the spatial parameters
};
} // end namespace Dumux
......
......@@ -24,8 +24,6 @@
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_VOLUME_VARIABLES_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_VOLUME_VARIABLES_HH
#include <dumux/common/properties.hh>
namespace Dumux {
/*!
......@@ -38,17 +36,6 @@ class ShallowWaterVolumeVariables
using Indices = typename Traits::ModelTraits::Indices;
using Scalar = typename Traits::PrimaryVariables::value_type;
enum {
// indices for primary variables
massBalanceIdx = Indices::massBalanceIdx,
momentumXBalanceIdx = Indices::momentumXBalanceIdx,
momentumYBalanceIdx = Indices::momentumYBalanceIdx,
waterdepthIdx = Indices::waterdepthIdx,
velocityXIdx = Indices::velocityXIdx,
velocityYIdx = Indices::velocityYIdx,
velocityOffset = Indices::velocityOffset
};
public:
using PrimaryVariables = typename Traits::PrimaryVariables;
......@@ -60,12 +47,11 @@ public:
{
priVars_ = elemSol[scv.localDofIndex()];
gravity_ = problem.spatialParams().gravityAtPos(element.geometry().center());
bedSurface_ = problem.spatialParams().bedSurface(element,scv);
}
/*!
* \brief Return the intrusion factor (dummy variable).
* \brief Return the extrusion factor (dummy variable).
*
*/
Scalar extrusionFactor() const
......@@ -77,7 +63,7 @@ public:
*/
Scalar waterDepth() const
{
return priVars_[waterdepthIdx];
return priVars_[Indices::waterdepthIdx];
}
/*!
......@@ -88,7 +74,7 @@ public:
Scalar velocity(int directionIndex) const
{
return priVars_[velocityOffset + directionIndex];
return priVars_[Indices::velocityOffset + directionIndex];
}
/*!
......@@ -100,20 +86,9 @@ public:
return bedSurface_;
}
/*!
* \brief Return the gravity constant.
*
*/
Scalar gravity() const
{
return gravity_;
}
private:
PrimaryVariables priVars_;
Scalar bedSurface_;
Scalar gravity_;
};
} // end namespace Dumux
......
......@@ -23,7 +23,6 @@
*/
#include <config.h>
#include "problem.hh"
#include <ctime>
#include <iostream>
......@@ -37,7 +36,6 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
......@@ -47,32 +45,7 @@
#include <dumux/assembly/fvassembler.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n";
std::cout << errorMessageOut
<< "\n";
}
}
#include "problem.hh"
////////////////////////
// the main function
......@@ -92,7 +65,7 @@ int main(int argc, char** argv) try
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv, usage);
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
......@@ -105,7 +78,6 @@ int main(int argc, char** argv) try
// 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);
......@@ -207,12 +179,12 @@ int main(int argc, char** argv) try
return 0;
}
catch (Dumux::ParameterException &e)
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;