Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
dumux-repositories
dumux
Commits
33fe841b
Commit
33fe841b
authored
Mar 04, 2020
by
Leopold Stadler
Committed by
Timo Koch
Mar 18, 2020
Browse files
[test][shallowwater] Add bowl test for new upwind mobility
parent
0be25839
Changes
7
Hide whitespace changes
Inline
Side-by-side
test/freeflow/shallowwater/CMakeLists.txt
View file @
33fe841b
add_subdirectory
(
bowl
)
add_subdirectory
(
dambreak
)
add_subdirectory
(
roughchannel
)
test/freeflow/shallowwater/bowl/CMakeLists.txt
0 → 100644
View file @
33fe841b
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"
)
test/freeflow/shallowwater/bowl/main.cc
0 → 100644
View file @
33fe841b
// -*- 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 (bowl).
*/
#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 <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.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
::
Bowl
;
// 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
GridGeometry
=
GetPropType
<
TypeTag
,
Properties
::
GridGeometry
>
;
auto
gridGeometry
=
std
::
make_shared
<
GridGeometry
>
(
leafGridView
);
gridGeometry
->
update
();
// the problem (initial and boundary conditions)
using
Problem
=
GetPropType
<
TypeTag
,
Properties
::
Problem
>
;
auto
problem
=
std
::
make_shared
<
Problem
>
(
gridGeometry
);
// the solution vector
using
SolutionVector
=
GetPropType
<
TypeTag
,
Properties
::
SolutionVector
>
;
SolutionVector
x
(
gridGeometry
->
numDofs
());
problem
->
applyInitialSolution
(
x
);
auto
xOld
=
x
;
// the grid variables
using
GridVariables
=
GetPropType
<
TypeTag
,
Properties
::
GridVariables
>
;
auto
gridVariables
=
std
::
make_shared
<
GridVariables
>
(
problem
,
gridGeometry
);
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"
);
problem
->
updateAnalyticalSolution
(
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
,
gridGeometry
,
gridVariables
,
timeLoop
,
xOld
);
// the linear solver
using
LinearSolver
=
AMGBiCGSTABBackend
<
LinearSolverTraits
<
GridGeometry
>>
;
auto
linearSolver
=
std
::
make_shared
<
LinearSolver
>
(
leafGridView
,
gridGeometry
->
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
);
//! Compute L2 error for the initial solution (should be zero)
problem
->
computeL2error
(
timeLoop
->
time
(),
x
,
*
gridVariables
);
// time loop
timeLoop
->
start
();
do
{
nonLinearSolver
.
solve
(
x
,
*
timeLoop
);
// make the new solution the old solution
xOld
=
x
;
gridVariables
->
advanceTimeStep
();
// update the analytical solution
problem
->
updateAnalyticalSolution
(
timeLoop
->
time
());
problem
->
computeL2error
(
timeLoop
->
time
(),
x
,
*
gridVariables
);
// 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
;
}
test/freeflow/shallowwater/bowl/params.input
0 → 100644
View file @
33fe841b
[Problem]
Name = bowl
BowlDepthAtCenter = 0.03 # [m]
BowlParaboloidRadius = 1.6 # [m]
BowlInitialWaterElevationAtCenter = 0.005 # [m]
Gravity = 9.81 # [m/s^2]
[FluxLimiterLET]
UpperWaterDepth = 1.0E-5
LowerWaterDepth = 1.0E-7
UpwindFluxLimiting = true
[TimeLoop]
TEnd = 13.104 # [s] 6.552 is one period
MaxTimeStepSize = 0.1 # [s]
DtInitial = 0.01 # [s]
[Grid]
Positions0 = -2.0 2.0
Positions1 = -2.0 2.0
Cells0 = 50
Cells1 = 50
[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
0 → 100644
View file @
33fe841b
// -*- 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 (bowl).
*/
#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.
*
* The domain is 4 meters long and 4 meters wide. The center of the domain is loacted at
* x = 0 and y = 0. There is no flow over the boundaries and no friction is considered.
*
* This example is demanding for the implicit model if a high mesh resolution is applied
* (e.g. 150x150 cells) in combination with a large time step size. Using the new limiting
* (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.
*
* 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
* in SWASHES (Shallow Water Analytic Solutions for Hydraulic and Environmental Studies,
* https://www.idpoisson.fr/swashes/).
*
* This problem uses the \ref ShallowWaterModel
*/
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
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
;
using
SolutionVector
=
GetPropType
<
TypeTag
,
Properties
::
SolutionVector
>
;
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
;
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
;
}
//! 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
;
}
}
//! Udpate the analytical solution
void
updateAnalyticalSolution
(
const
Scalar
time
)
{
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
]);
}
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \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
;
}
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \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
{
BoundaryTypes
bcTypes
;
bcTypes
.
setAllNeumann
();
return
bcTypes
;
}
/*!
* \brief Specifies the neumann boundary
*
* 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
*/
NeumannFluxes
neumann
(
const
Element
&
element
,
const
FVElementGeometry
&
fvGeometry
,
const
ElementVolumeVariables
&
elemVolVars
,
const
ElementFluxVariablesCache
&
elemFluxVarsCache
,
const
SubControlVolumeFace
&
scvf
)
const
{
NeumannFluxes
values
(
0.0
);
const
auto
&
insideScv
=
fvGeometry
.
scv
(
scvf
.
insideScvIdx
());
const
auto
&
insideVolVars
=
elemVolVars
[
insideScv
];
const
auto
&
nxy
=
scvf
.
unitOuterNormal
();
const
auto
gravity
=
this
->
spatialParams
().
gravity
(
scvf
.
center
());
std
::
array
<
Scalar
,
3
>
boundaryStateVariables
;
//no flow with zero normal velocity and tangential velocity
const
auto
vNormalGhost
=
-
(
nxy
[
0
]
*
insideVolVars
.
velocity
(
0
)
+
nxy
[
1
]
*
insideVolVars
.
velocity
(
1
));
const
auto
vTangentialGhost
=
-
nxy
[
1
]
*
insideVolVars
.
velocity
(
0
)
+
nxy
[
0
]
*
insideVolVars
.
velocity
(
1
);
boundaryStateVariables
[
0
]
=
insideVolVars
.
waterDepth
();
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
],
insideVolVars
.
bedSurface
(),
insideVolVars
.
bedSurface
(),
gravity
,
nxy
);
values
[
Indices
::
massBalanceIdx
]
=
riemannFlux
[
0
];
values
[
Indices
::
velocityXIdx
]
=
riemannFlux
[
1
];
values
[
Indices
::
velocityYIdx
]
=
riemannFlux
[
2
];
return
values
;
}