Commit 194c56fb authored by Leopold Stadler's avatar Leopold Stadler Committed by Timo Koch
Browse files

[shallowwater] added freeflow/shallowwater model

parent 87b6a923
......@@ -2,6 +2,7 @@ add_subdirectory(compositional)
add_subdirectory(navierstokes)
add_subdirectory(nonisothermal)
add_subdirectory(rans)
add_subdirectory(shallowwater)
install(FILES
properties.hh
......
#install headers
install(FILES
boundaryfluxes.hh
fluxvariables.hh
fluxvariablescache.hh
indices.hh
localresidual.hh
model.hh
properties.hh
volumevariables.hh
vtkoutputfields.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflox/shallowwater/)
// -*- 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 ShallowWater
* \brief Compute the boundary fluxes based on the Riemann invariants
*
* bdType can be:
* 0 = noFlow (do nothing)
* 1 = given discharge
* 2/3 = given h (if 3, h comes from a stage discharge curve)
*
*
*/
#ifndef DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#define DUMUX_SHALLOWWATER_BOUNDARYFLUXES_HH
#include <dumux/common/math.hh>
namespace Dumux {
namespace ShallowWater {
inline void hBoundary(const auto& nxy,
const auto& faceVolume,
const auto& minH,
const auto& cellStatesLeft,
auto& cellStatesRight,
const auto& bdType,
const auto& bdValue)
{
cellStatesRight[0] = bdValue - cellStatesRight[3];
auto uboundIn = nxy[0] * cellStatesLeft[1] + nxy[1] * cellStatesLeft[2] ;
auto uboundQut = uboundIn + 2.0 * std::sqrt(9.81 * cellStatesLeft[0]) - 2.0 * std::sqrt(9.81 * cellStatesRight[0]);
cellStatesRight[1] = (nxy[0] * uboundQut); // - nxy[1] * (-nxy[1] * cellStatesLeft[1] + nxy[0] * cellStatesLeft[1]));
cellStatesRight[2] = (nxy[1] * uboundQut); // + nxy[0] * (-nxy[1] * cellStatesLeft[1] + nxy[0] * cellStatesLeft[1]));
}
inline void qBoundary(const auto& nxy,
const auto& faceVolume,
const auto& minH,
const auto& cellStatesLeft,
auto& cellStatesRight,
const auto& bdType,
const auto& bdValue)
{
std::array<Scalar,3> cellStateRight;
using std::pow;
using std::abs;
using std::sqrt;
//olny impose if abs(q) > 0
if (std::abs(bdValue) > 1.0E-9){
auto qlocal = (bdValue) /faceVolume;
auto uboundIn = nxy[0] * cellStatesLeft[1] + nxy[1] * cellStatesLeft[2];
auto alphal = uboundIn + 2.0 * std::sqrt(9.81 * cellStatesLeft[0]);
//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;
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 (std::pow(dx_hstar,2.0) < tol_hstar){
break;
}
}
auto qinner = (nxy[0] * cellStatesLeft[0] * cellStatesLeft[2]) - (nxy[1] * cellStatesLeft[0] * cellStatesLeft[1]);
cellStatesRight[0] = hstar;
cellStatesRight[1] = (nxy[0] * qlocal - nxy[1] * qinner)/hstar;
cellStatesRight[2] = (nxy[1] * qlocal + nxy[0] * qinner)/hstar;
}
}
} // end namespcae ShallowWater
} // end namespace Dumux
#endif
// -*- 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 ShallowWaterModel
* \copydoc Dumux::ShallowWaterFluxVariables
*/
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_FLUXVARIABLES_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_FLUXVARIABLES_HH
#include <dumux/common/properties.hh>
#include <dumux/flux/fluxvariablesbase.hh>
namespace Dumux{
/*!
* \ingroup ShallowWaterModel
* \brief The flux variables class for the shallow water model.
*
*/
template<class TypeTag>
class ShallowWaterFluxVariables
: 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>;
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
//using DiffusionType = GetPropType<TypeTag, Properties::DiffusionType>;
static constexpr bool enableAdvection = ModelTraits::enableAdvection();
static constexpr bool enableDiffusion = ModelTraits::enableDiffusion();
public:
/*!
* \brief Returns the advective flux computed by the Riemann solver
*
*/
NumEqVector advectiveFlux() const
{
NumEqVector fluxVector(0.0);
if (enableAdvection)
{
return AdvectionType::flux(this->problem(),
this->element(),
this->fvGeometry(),
this->elemVolVars(),
this->scvFace(),
this->elemFluxVarsCache());
}
else
{
return fluxVector;
}
}
/*!
* \brief Returns the diffusive flux (e.g. diffusion of tracer)
*
*/
NumEqVector diffusiveFlux() const
{
NumEqVector fluxVector(0.0);
if (enableDiffusion)
{
/*return DiffusionType::flux(this->problem(),
this->element(),
this->fvGeometry(),
this->elemVolVars(),
this->scvFace(),
this->elemFluxVarsCache());
*/
return fluxVector;
}
else
{
return fluxVector;
}
}
};
} // end namespace Dumux
#endif
// -*- 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 ShallowWaterModel
* \brief Base class for the flux variables
*/
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_FLUXVARIABLESCACHE_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_FLUXVARIABLESCACHE_HH
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/box/fluxvariablescache.hh>
#include <dumux/flux/fluxvariablescaching.hh>
namespace Dumux {
// forward declaration
template<class TypeTag, DiscretizationMethod discMethod>
class ShallowWaterFluxVariablesCacheImplementation;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
//! The cache is dependent on the active physical processes (advection, diffusion, heat conduction)
//! For each type of process there is a base cache storing the data required to compute the respective fluxes
//! Specializations of the overall cache are provided for combinations of processes
//!
//! Caches are not used for the Shallow Water Model since we use a Riemann Solver for all kind of fluxes
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
/*!
* \ingroup ImplicitModel
* \brief The flux variables cache classes.
* Store data required for flux calculation. For each type of physical process (advection, diffusion, heat conduction)
* there is a base cache storing the data required to compute the respective fluxes. Specializations of the overall
* cache class are provided for different combinations of processes.
*/
template<class TypeTag>
using ShallowWaterFluxVariablesCache = ShallowWaterFluxVariablesCacheImplementation<TypeTag, GetPropType<TypeTag, Properties::FVGridGeometry>::discMethod>;
// the following classes choose the cache type: empty if the law disabled and the law's cache if it's enabled
// if advections is disabled the advection type is still instatiated if we use std::conditional_t and has to be a full type
// in order to prevent that instead of std::conditional_t we use this helper type is only dependent on the advection type
// if advection is enabled otherwise its an empty cache type
template<class TypeTag, bool EnableAdvection> class AdvectionCacheChooser : public FluxVariablesCaching::EmptyAdvectionCache {};
template<class TypeTag> class AdvectionCacheChooser<TypeTag, true> : public GetPropType<TypeTag, Properties::AdvectionType>::Cache {};
template<class TypeTag, bool EnableDiffusion> class DiffusionCacheChooser : public FluxVariablesCaching::EmptyDiffusionCache {};
template<class TypeTag> class DiffusionCacheChooser<TypeTag, true> : public GetPropType<TypeTag, Properties::AdvectionType>::Cache {};
// specialization for the cell centered tpfa method
template<class TypeTag>
class ShallowWaterFluxVariablesCacheImplementation<TypeTag, DiscretizationMethod::cctpfa>
{};
} // end namespace Dumux
#endif
// -*- 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 3 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 ShallowWaterModel
* \brief A helper class to have a dummy filler for the flux variable cache
*/
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_FLUXVARSCACHE_FILLER_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_FLUXVARSCACHE_FILLER_HH
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
namespace Dumux {
/*!
* \ingroup ShallowWaterModel
* \brief A helper class to have a dummy filler for the flux variables cache
*/
template<class TypeTag>
class ShallowWaterFluxVariablesCacheFiller
{
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using FluxVariablesCache = GetPropType<TypeTag, Properties::FluxVariablesCache>;
using Element = typename GridView::template Codim<0>::Entity;
public:
static constexpr bool isSolDependent = (true);
//! The constructor. Sets the problem pointer
ShallowWaterFluxVariablesCacheFiller(const Problem& problem) : problemPtr_(&problem) {}
/*!
* \brief function to fill the flux variables caches
*
* \param fluxVarsCacheContainer Either the element or global flux variables cache
* \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param elemVolVars The element volume variables
* \param scvf The corresponding sub-control volume face
* \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones)
*/
template<class FluxVariablesCacheContainer>
void fill(FluxVariablesCacheContainer& fluxVarsCacheContainer,
FluxVariablesCache& scvfFluxVarsCache,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
bool forceUpdateAll = false)
{
}
private:
const Problem& problem() const
{ return *problemPtr_; }
const Problem* problemPtr_;
};
} // end namespace Dumux
#endif
// -*- 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 ShallowWaterModel
* \copydoc Dumux::ShallowWaterIndices
*/
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_INDICES_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_INDICES_HH
#include <dumux/common/properties.hh>
namespace Dumux{
// \{
/*!
* \ingroup ShallowWaterModel
* \brief The common indices for the shallow water equations model.
*
* \tparam PVOffset The first index in a primary variable vector.
*/
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
static constexpr int massBalanceIdx = 0; //!< Index of the mass balance equation
static constexpr int momentumXBalanceIdx = 1; //!< Index of the x momentum balance equation
static constexpr int momentumYBalanceIdx = 2; //!< Index of the y momentum balance equation
static constexpr int waterdepthIdx = massBalanceIdx; //!< Index of the velocity in a solution vector
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
};
// \}
} // end namespace Dumux
#endif
// -*- 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 ShallowWaterModel
* \copydoc Add I/O fields specific to shallow water
*/
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_IO_FIELDS_HH
#define DUMUX_FREEFLOW_SHALLOW_WATER_IO_FIELDS_HH
#include <dumux/io/name.hh>
#include <string>
#include <dune/common/fvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
namespace Dumux{
/*!
* \ingroup ShallowWaterModel
* \brief Adds vtk output fields for the shallow water model
*/
class ShallowWaterIOFields
{
public:
template <class OutputModule>
static void initOutputModule(OutputModule& out)
{
using VolumeVariables = typename OutputModule::VolumeVariables;
out.addVolumeVariable([](const VolumeVariables& v){ return v.waterDepth(); }, "waterDepth");
out.addVolumeVariable([](const VolumeVariables& v){ return v.velocity(0); }, "velocityX");
out.addVolumeVariable([](const VolumeVariables& v){ return v.velocity(1); }, "velocityY");
out.addVolumeVariable([](const VolumeVariables& v){ return v.bedSurface(); }, "bedSurface");
out.addVolumeVariable([](const VolumeVariables& v){ return v.bedSurface() + v.waterDepth(); }, "freeSurface");
}
/* Restart is limited for shallow water models since only primary
* variables are regarded. Important parameters like the bedSurface,
* fricition,.. are missing so far.
*/
template <class ModelTraits>
static std::string primaryVariableName(int pvIdx)
{
std::string name;
switch(pvIdx){
case 0 : name = "waterDepth";
case 1 : name = "velocityX";
case 2 : name = "velocityY";
}
return name;
}
};
} // end namespace Dumux
#endif
// -*- 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 ShallowWaterModel
* \copydoc Dumux::ShallowWaterResidual
*/
#ifndef DUMUX_FREEFLOW_SHALLOW_WATER_LOCAL_RESIDUAL_HH
#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{
/*!
* \ingroup ShallowWaterModel
* \brief Element-wise calculation of the residual for the shallow water equations
*/
template<class TypeTag>
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>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
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 DiffusionType = GetPropType<TypeTag, Properties::DiffusionType>;
public:
using ParentType::ParentType;
/*!
* \brief Evaluate the rate of change of all conservation
* quantites (e.g. mass, momentum) within a sub-control
* volume of a finite volume element.
* \param problem The problem