Skip to content

GitLab

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

[test][shallowwater] Cleanup bowl test, run in parallel too

parent 33fe841b
Expand all Hide whitespace changes
Inline Side-by-side
test/freeflow/shallowwater/bowl/CMakeLists.txt
View file @ a3f39818
dune_symlink_to_source_files(FILES "params.input")
dumux_add_test(NAME test_shallowwater_bowl
SOURCES main.cc
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_bowl-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/bowl-00001.vtu
--zeroThreshold {"velocityY":1e-14}
--command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_bowl params.input
-Problem.Name bowl")
SOURCES main.cc
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_bowl-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/bowl-00013.vtu
--zeroThreshold {"velocityY":1e-14}
--command "${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_bowl")
dumux_add_test(NAME test_shallowwater_bowl_parallel
TARGET test_shallowwater_bowl
CMAKE_GUARD MPI_FOUND
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_shallowwater_bowl-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/s0002-bowl-parallel-00013.pvtu
--zeroThreshold {"velocityY":1e-14,"process rank":100}
--command "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater_bowl -Problem.Name bowl-parallel")
test/freeflow/shallowwater/bowl/main.cc
View file @ a3f39818
......@@ -42,7 +42,83 @@
#include <dumux/assembly/fvassembler.hh>
#include "problem.hh"
#include "properties.hh"
template<class SolutionVector, int i>
class SolutionComponent
{
const SolutionVector& sol_;
public:
SolutionComponent(const SolutionVector& sol) : sol_(sol) {}
auto operator[] (std::size_t idx) const { return sol_[idx][i]; }
auto size() const { return sol_.size(); }
};
//! Compute the analytical solution at time t
template<class SolutionVector, class Problem>
SolutionVector computeAnalyticalSolution(const double t, const Problem& problem)
{
const auto& gg = problem.gridGeometry();
SolutionVector exactSolution(gg.numDofs());
for (const auto& element : elements(gg.gridView()))
{
const auto eIdx = gg.elementMapper().index(element);
const auto& globalPos = element.geometry().center();
exactSolution[eIdx] = problem.analyticalSolution(t, globalPos);
}
return exactSolution;
}
//! Compute L2 error wrt the analytical solution at time t
template<class SolutionVector, class GridGeometry>
typename SolutionVector::block_type
computeL2Error(const double t,
const SolutionVector& exactSolution,
const SolutionVector& curSol,
const GridGeometry& gridGeometry)
{
typename SolutionVector::block_type l2Error(0.0);
for (const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
{
auto fvGeometry = localView(gridGeometry);
fvGeometry.bindElement(element);
using std::pow;
for (auto&& scv : scvs(fvGeometry))
{
auto localDiffSq = exactSolution[scv.dofIndex()] - curSol[scv.dofIndex()];
for (int i = 0; i < localDiffSq.size(); ++i)
localDiffSq[i] *= localDiffSq[i];
l2Error += localDiffSq;
}
}
// sum over processes if we are running with multiple processes in parallel
if (gridGeometry.gridView().comm().size() > 0)
l2Error = gridGeometry.gridView().comm().sum(l2Error);
// square root of sum of squared errors is our absolute discrete l2 error
for (int i = 0; i < l2Error.size(); ++i)
l2Error[i] = std::sqrt(l2Error[i]);
return l2Error;
}
//! Compute and print L2 error wrt the analytical solution at time t
template<class SolutionVector, class GridGeometry>
void computeAndPrintL2Error(const double t,
const SolutionVector& exactSolution,
const SolutionVector& curSol,
const GridGeometry& gridGeometry)
{
const auto l2Error = computeL2Error(t, exactSolution, curSol, gridGeometry);
auto numElements = gridGeometry.gridView().size(0);
// sum over processes if we are running with multiple processes in parallel
if (gridGeometry.gridView().comm().size() > 0)
numElements = gridGeometry.gridView().comm().sum(numElements);
if (gridGeometry.gridView().comm().rank() == 0)
std::cout << "L2 error in (h, u, v) at t = " << t << " seconds for " << numElements << " elements: "
<< std::scientific << l2Error << std::endl;
}
////////////////////////
// the main function
......@@ -86,7 +162,7 @@ int main(int argc, char** argv) try
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
SolutionVector x;
problem->applyInitialSolution(x);
auto xOld = x;
......@@ -97,16 +173,20 @@ int main(int argc, char** argv) try
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto tEnd = 2.0*problem->oscillationPeriodInSeconds();
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
// intialize the vtk output module and analytical solution
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
problem->updateAnalyticalSolution(0.0);
auto exactSolution = computeAnalyticalSolution<SolutionVector>(0.0, *problem);
const auto exactWaterDepth = SolutionComponent<SolutionVector, 0>(exactSolution);
const auto exactVelocityX = SolutionComponent<SolutionVector, 1>(exactSolution);
const auto exactVelocityY = SolutionComponent<SolutionVector, 2>(exactSolution);
vtkWriter.addField(exactWaterDepth, "exactWaterDepth");
vtkWriter.addField(exactVelocityX, "exactVelocityX");
vtkWriter.addField(exactVelocityY, "exactVelocityY");
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter);
vtkWriter.write(0.0);
......@@ -126,11 +206,8 @@ int main(int argc, char** argv) try
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
//! set some check point at the end of the time loop
timeLoop->setCheckPoint(tEnd);
//! Compute L2 error for the initial solution (should be zero)
problem->computeL2error(timeLoop->time(), x, *gridVariables);
//! Compute L2 error for the initial solution (should be zero because initial solution is exact)
computeAndPrintL2Error(0.0, exactSolution, x, *gridGeometry);
// time loop
timeLoop->start(); do
......@@ -141,15 +218,15 @@ int main(int argc, char** argv) try
xOld = x;
gridVariables->advanceTimeStep();
// update the analytical solution
problem->updateAnalyticalSolution(timeLoop->time());
problem->computeL2error(timeLoop->time(), x, *gridVariables);
// update the analytical solution and print current l2 error
exactSolution = computeAnalyticalSolution<SolutionVector>(timeLoop->time(), *problem);
computeAndPrintL2Error(timeLoop->time(), exactSolution, x, *gridGeometry);
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
if (timeLoop->isCheckPoint())
// write vtk output every 10 time steps
if (!(timeLoop->timeStepIndex() % 10))
vtkWriter.write(timeLoop->time());
// report statistics of this time step
......@@ -176,28 +253,14 @@ int main(int argc, char** argv) try
return 0;
}
// Error handling
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;
return 2;
}
test/freeflow/shallowwater/bowl/params.input
View file @ a3f39818
......@@ -11,29 +11,14 @@ LowerWaterDepth = 1.0E-7
UpwindFluxLimiting = true
[TimeLoop]
TEnd = 13.104 # [s] 6.552 is one period
MaxTimeStepSize = 0.1 # [s]
DtInitial = 0.01 # [s]
DtInitial = 0.1 # [s]
[Grid]
Positions0 = -2.0 2.0
Positions1 = -2.0 2.0
Cells0 = 50
Cells1 = 50
Cells0 = 40
Cells1 = 40
[Newton]
EnableAbsoluteResidualCriterion = true
MaxAbsoluteResidual = 1e-6
RetryTimeStepReductionFactor = 0.5
SatisfyResidualAndShiftCriterion = "false"
EnableDynamicOutput = false
EnablePartialReassembly = false
TargetSteps = 6
MaxSteps = 10
UseLineSearch = true
[LinearSolver]
MaxIterations = 200
ResidualReduction = 1e-5
test/freeflow/shallowwater/bowl/problem.hh
View file @ a3f39818
......@@ -24,59 +24,11 @@
#ifndef DUMUX_BOWL_TEST_PROBLEM_HH
#define DUMUX_BOWL_TEST_PROBLEM_HH
#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>
namespace Dumux {
template <class TypeTag>
class BowlProblem;
// Specify the properties for the problem
namespace Properties {
// Create new type tags
namespace TTag {
struct Bowl { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
} // end namespace TTag
template<class TypeTag>
struct Grid<TypeTag, TTag::Bowl>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::Bowl>
{ using type = Dumux::BowlProblem<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::Bowl>
{
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
public:
using type = BowlSpatialParams<GridGeometry, Scalar, VolumeVariables>;
};
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::Bowl>
{ static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::Bowl>
{ static constexpr bool value = false; };
} // end namespace Properties
/*!
* \ingroup ShallowWaterTests
* \brief A wetting and drying test with sloshing water in a bowl.
......@@ -89,11 +41,14 @@ struct EnableGridVolumeVariablesCache<TypeTag, TTag::Bowl>
* (UpwindFluxLimiting = true) will help to improve the convergence for such cases.
*
* This test uses a low mesh resoultion and only ensures that UpwindFluxLimiting for the mobility
* works.
* works. For low mesh resolution the solution is very diffusive so that the oscillation is dampened.
* This gets better with grid refinement (not tested here).
*
* The results are checked against a analytical solution which is based on the "Thacker-Solution"
* (William Thacker, "Some exact solutions to the nonlinear shallow-water wave equations", Journal
* of Fluid Mechanics, 107:499–508, 1981). Further examples and details on the solution are given
* of Fluid Mechanics, 107:499–508, 1981, doi: https://doi.org/10.1017/S0022112081001882).
* This implements the oscillating solution in a circular bowl (Section 4 in the paper).
* Further examples and details on the solution are given
* in SWASHES (Shallow Water Analytic Solutions for Hydraulic and Environmental Studies,
* https://www.idpoisson.fr/swashes/).
*
......@@ -103,125 +58,93 @@ template <class TypeTag>
class BowlProblem : 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 GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
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 Scalar = GetPropType<TypeTag, Properties::Scalar>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
public:
BowlProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
using std::sqrt;
using std::pow;
name_ = getParam<std::string>("Problem.Name");
gravity_ = getParam<Scalar>("Problem.Gravity");
bowlDepthAtCenter_ = getParam<Scalar>("Problem.BowlDepthAtCenter");
bowlParaboloidRadius_ = getParam<Scalar>("Problem.BowlParaboloidRadius");
bowlInitialWaterElevationAtCenter_ = getParam<Scalar>("Problem.BowlInitialWaterElevationAtCenter");
bowlAnalyticParameterOmega_ = sqrt(8.0 * gravity_ * bowlDepthAtCenter_) / bowlParaboloidRadius_;
bowlAnalyticParameterA_ = (pow(bowlDepthAtCenter_ + bowlInitialWaterElevationAtCenter_, 2.0) -
pow(bowlDepthAtCenter_, 2.0)) /
(pow(bowlDepthAtCenter_ + bowlInitialWaterElevationAtCenter_, 2.0)+
pow(bowlDepthAtCenter_, 2.0));
exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0);
updateAnalyticalSolution(0.0);
}
//! Get the analytical water depth
const std::vector<Scalar>& getExactWaterDepth()
{
return exactWaterDepth_;
}
//! Get the analctic water depth at
Scalar calculateAnalyticWaterDepth(Scalar time, Scalar x, Scalar y)
{
using std::max;
// Thacker (1981) Eq. (43)
using std::sqrt;
using std::pow;
using std::cos;
auto waterDepth = max(bowlDepthAtCenter_* ((sqrt(1.0 - pow(bowlAnalyticParameterA_, 2.0)))/
(1.0 - bowlAnalyticParameterA_ * cos(bowlAnalyticParameterOmega_ * time))-
1.0 - (pow(x, 2.0) + pow(y, 2.0)) / pow(bowlParaboloidRadius_, 2.0) *
((1.0- pow(bowlAnalyticParameterA_, 2.0)) /
(pow(1.0 - bowlAnalyticParameterA_ *
cos(bowlAnalyticParameterOmega_ * time), 2.0)) - 1.0)) -
bowlDepthAtCenter_ *((pow(x, 2.0) + pow(y, 2.0))/
pow(bowlParaboloidRadius_, 2.0) - 1.0), 1.0E-5);
return waterDepth;
bowlAnalyticParameterOmega_ = sqrt(8.0 * getParam<Scalar>("Problem.Gravity") * bowlDepthAtCenter_) / bowlParaboloidRadius_;
std::cout << "One full oscillation period of the water table is: "
<< oscillationPeriodInSeconds() << " seconds." << std::endl;
// Thacker (1981) Eq. (50)
const auto D0PlusEta = bowlDepthAtCenter_ + getParam<Scalar>("Problem.BowlInitialWaterElevationAtCenter");
const auto D0PlusEtaSquared = D0PlusEta*D0PlusEta;
const auto D0Squared = bowlDepthAtCenter_*bowlDepthAtCenter_;
bowlAnalyticParameterA_ = (D0PlusEtaSquared - D0Squared)/(D0PlusEtaSquared + D0Squared);
// check constraint Thacker (1981) text after Eq. (49)
if (bowlAnalyticParameterA_ >= 1.0)
DUNE_THROW(Dune::InvalidStateException, "Parameter A has to be smaller than unity!");
}
//! Compute L2 error
void computeL2error(const Scalar time,
const SolutionVector& curSol,
const GridVariables& gridVariables)
{
Scalar l2error = 0.0;
//first ensure that the Analytical solution is up-to-date
updateAnalyticalSolution(time);
for (const auto& element : elements(this->gridGeometry().gridView()))
{
const auto eIdx = this->gridGeometry().elementMapper().index(element);
const auto& globalPos = element.geometry().center();
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, curSol);
exactWaterDepth_[eIdx] = calculateAnalyticWaterDepth(time, globalPos[0], globalPos[1]);
for (auto&& scv : scvs(fvGeometry))
{
using std::pow;
l2error += pow(exactWaterDepth_[eIdx] - elemVolVars[scv].waterDepth(), 2.0);
}
}
using std::sqrt;
l2error = sqrt(l2error);
l2error = this->gridGeometry().gridView().comm().sum(l2error);
if (this->gridGeometry().gridView().comm().rank() == 0)
{
std::cout << "L2 error at t = "
<< time << " seconds "
<< " for "
<< std::setw(6) << this->gridGeometry().gridView().size(0)
<< " elements: "
<< std::scientific
<< l2error
<< std::endl;
}
}
//! One oscillation period of the water table (analytically this goes on forever)
Scalar oscillationPeriodInSeconds() const
{ return 2*M_PI/bowlAnalyticParameterOmega_; }
//! Udpate the analytical solution
void updateAnalyticalSolution(const Scalar time)
//! Get the analytical water depth at time t and position pos
PrimaryVariables analyticalSolution(const Scalar t, const GlobalPosition& pos) const
{
for (const auto& element : elements(this->gridGeometry().gridView()))
{
const auto eIdx = this->gridGeometry().elementMapper().index(element);
const auto& globalPos = element.geometry().center();
exactWaterDepth_[eIdx] = calculateAnalyticWaterDepth(time, globalPos[0], globalPos[1]);
}
using std::sqrt;
using std::cos;
using std::sin;
// see Thacker (1981) Eq. (51) for formula
const auto radiusSquared = pos[0]*pos[0] + pos[1]*pos[1];
const auto LSquared = bowlParaboloidRadius_*bowlParaboloidRadius_;
const auto A = bowlAnalyticParameterA_;
const auto omega = bowlAnalyticParameterOmega_;
const auto D0 = bowlDepthAtCenter_;
const auto oneMinusASq = 1.0 - A*A;
const auto oneMinusACosOmegaT = 1.0 - A*cos(omega*t);
const auto ratioSq = oneMinusASq / (oneMinusACosOmegaT*oneMinusACosOmegaT);
const auto localRadiusSq = radiusSquared / LSquared;
// bowl depth function (cf. D in Thacker (1981))
const auto D = D0*(1.0 - localRadiusSq);
// height above equilibrium water level (cf. h in Thacker (1981))
const auto h = D0*(sqrt(ratioSq) - 1.0 - localRadiusSq*(ratioSq - 1.0));
// see remark about total water depth in Thacker (1981) beginning section 2
const auto analyticalWaterDepth = h + D;
const auto halfOmegaASinOmegaT = 0.5*omega*A*sin(omega*t);
const auto analyticalVelocityX = pos[0]*halfOmegaASinOmegaT/oneMinusACosOmegaT;
const auto analyticalVelocityY = pos[1]*halfOmegaASinOmegaT/oneMinusACosOmegaT;
// The radius of the shoreline (where h + D = 0), Eq. (48)
const auto h0 = D0*(sqrt(ratioSq) - 1.0); // h in the middle of the bowl (r=0)
const auto shoreLineRadiusSquared = LSquared*(D0/(D0 + h0));
// outside shoreline the water height and velocity is zero
if (radiusSquared > shoreLineRadiusSquared)
return { 0.0, 0.0, 0.0 };
else
return { analyticalWaterDepth, analyticalVelocityX, analyticalVelocityY };
}
/*!
......@@ -231,41 +154,10 @@ public:
/*!
* \brief The problem name
*
* This is used as a prefix for files generated by the simulation.
*/
const std::string& name() const
{
return name_;
}
/*!
* \brief Evaluate the source term for all balance equations within a given
* sub-control-volume.
*
* This is the method for the case where the source term is
* potentially solution dependent and requires some quantities that
* are specific to the fully-implicit method.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param scv The sub control volume
*
* For this method, the \a values parameter stores the conserved quantity rate
* generated or annihilate per volume unit. Positive values mean
* that the conserved quantity is created, negative ones mean that it vanishes.
* 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
{
NumEqVector source (0.0);
return source;
}
{ return name_; }
// \}
......@@ -277,8 +169,6 @@ public:
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param globalPos The position for which the boundary type is set
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
......@@ -293,12 +183,8 @@ public:
* We need the Riemann invariants to compute the values depending of the boundary type.
* Since we use a weak imposition we do not have a dirichlet value. We impose fluxes
* based on q, h, etc. computed with the Riemann invariants
*
* \param element
* \param fvGeometry
* \param elemVolVars
* \param scvf
*/
template<class ElementFluxVariablesCache>
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
......@@ -321,16 +207,12 @@ public:
boundaryStateVariables[1] = nxy[0] * vNormalGhost - nxy[1] * vTangentialGhost;
boundaryStateVariables[2] = nxy[1] * vNormalGhost + nxy[0] * vTangentialGhost;
auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
boundaryStateVariables[0],
insideVolVars.velocity(0),
boundaryStateVariables[1],
insideVolVars.velocity(1),
boundaryStateVariables[2],