Commit 9839303e authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'master' into 92-add-a-logging-framework-for-command-line-output

parents a1ad48f8 113155f9
......@@ -20,6 +20,9 @@
* Code coverage reports for unit tests. Detailed reports can be retrieved from
the artifacts of the `test:unit-tests` test in the folder
`build-cmake/test/coverage`. The total coverage is reported to GitLab.
* New class [`SimulationTransport`](dune/dorie/interface/transport_simulation.hh)
which solves the transport equation with a finite volume scheme for give
grid functions of water content and velocity.
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......@@ -111,6 +114,10 @@
It dumps its data into a YAML file which is then loaded for writing
the default configuration files and the cheat sheet.
### Fixed
* Solver in `RichardsSimulation` was using the wrong time variable.
[!116](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/116)
### Deprecated
* The configuration file key `[parameters.interpolation]` is deprecated due to
the new scheme for storing parameterization data. DORiE now only supports
......@@ -135,6 +142,8 @@
CI variable `REBUILD_BASE_IMAGE` to a non-empty value. Previously, the
builds in stage `setup` where triggered by any manual pipeline run.
* The test system folder was renamed from `testing` to [`test`](test).
* File `operator_DG.hh` is renamed
[`richards_operator_DG.hh`](dune/dorie/solver/richards_operator_DG.hh)
### Fixed
* Removed VTK output from `dorie_mass_conservation` executable, which would
......@@ -198,4 +207,4 @@
## 1.0.0 (2018-03-28)
First stable version.
\ No newline at end of file
First stable version.
......@@ -158,3 +158,20 @@ create_default_config(
${PROJECT_SOURCE_DIR}/dune/dorie
${CMAKE_CURRENT_SOURCE_DIR}/parameters.css
)
# Model: Transport
scrape_parameters(
SOURCE_DIR ${PROJECT_SOURCE_DIR}/dune/dorie/model/transport
XML_FILE ${CMAKE_CURRENT_SOURCE_DIR}/transport-parameters.xml
OUTPUT config.yml
MODEL "transport"
)
# Transport config (temporally hidden to the user)
# TODO: Include in API once transport is ready
create_default_config(
${CMAKE_CURRENT_BINARY_DIR}/config.yml
"transport-config.ini;transport-parameters.html;transport-parameters.rst"
${PROJECT_SOURCE_DIR}/dune/dorie
${CMAKE_CURRENT_SOURCE_DIR}/parameters.css
)
<?xml version="1.0" encoding="UTF-8"?>
<!--
If you want to use any special characters, you will need to define them here.
A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
-->
<!DOCTYPE naughtyxml [
<!ENTITY alpha "&#945;">
<!ENTITY beta "&#946;">
<!ENTITY eta "&#951;">
<!ENTITY tau "&#964;">
<!ENTITY times "&#215;">
]>
<!--
XML file hierarchy:
<dorie> -> <category> -> <parameter> -> (parameter attributes)
Possible parameter attributes:
definition: meaning of the parameter, will only show up in html output
suggestion: standard value in created parameter files
values: possible values, will only show up in html output
comment: extra comment, will only show up in parameter files
All attributes are optional.
The parser supports rudimentary markdown / styling. You can add a paragraph by
adding an empty line, make text **bold** or ``monospaced``.
-->
<dorie>
<category name="output">
<parameter name="verbose">
<definition> Overall verbosity of the program </definition>
<suggestion> 0 </suggestion>
<values> 0, 1, 2, 3 </values>
</parameter>
<parameter name="outputPath">
<definition> Path to the directory where most of the outputs are stored.
DORiE will attempt to create it if it does not exist. </definition>
<suggestion> ./ </suggestion>
<values> path </values>
</parameter>
<parameter name="fileName">
<definition> Base file name for VTK output. </definition>
<values> string </values>
</parameter>
<parameter name="vertexData">
<definition> Plot vertex based (``true``) or cell-centered (``false``)
data into VTK files. Vertex based data might render sharp
parameterization boundaries inappropriately.
System tests and plotting functions (``dorie plot``) require
cell-centered data.
</definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="subsamplingLevel">
<definition> Plot VTK files with virtually refined grids. VTK only
supports bilinear triangulations and displays higher-order solutions
inappropriately. Use level 0 for order 1, and level N for order N.
For level &gt; 0, the printed grid does not resemble the actual grid.
This parameter defaults to 0 if not given in the config file. Notice
that subsampling significantly increases the output file size!
</definition>
<values> int </values>
<suggestion> 0 </suggestion>
</parameter>
<parameter name="asciiVtk">
<definition> Defines whether VTK files should be written as ASCII (``true``)
or binary (``false``). ASCII is easier to parse in case you want to write
your own post-processing, but takes a lot more space on your hard drive.
</definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
</category>
<category name="boundary">
<parameter name="file">
<definition> Path to the boundary condition file. </definition>
<values> path </values>
</parameter>
<parameter name="fileType">
<definition> Type of spatial segmentation of the boundaries specified in the BC file </definition>
<values> rectangularGrid </values>
<suggestion> rectangularGrid </suggestion>
<comment> Choose type of boundary segmentation: rectangularGrid </comment>
</parameter>
<parameter name="interpolateBCvalues">
<definition> Whether to interpolate between the boundary conditions
at different times linearly (``true``) or not at all (``false``). May require
different boundary condition files. </definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
</category>
<category name="initial">
<parameter name="value">
<definition> Initial constant value over the whole domain. </definition>
<values> float </values>
<suggestion> 0 </suggestion>
</parameter>
</category>
<category name="time">
<parameter name="start">
<definition> Starting time in seconds. </definition>
<values> float </values>
<suggestion> 0 </suggestion>
</parameter>
<parameter name="end">
<definition> Ending time in seconds. </definition>
<values> float </values>
<suggestion> 1E6 </suggestion>
</parameter>
<parameter name="minTimestep">
<definition> Minimum time step that is allowed before DORiE stops running.
with an error, in seconds. </definition>
<values> float </values>
<suggestion> 0.1 </suggestion>
</parameter>
<parameter name="startTimestep">
<definition> Value of the first time step in seconds. </definition>
<values> float </values>
<suggestion> 10 </suggestion>
</parameter>
<parameter name="maxTimestep">
<definition> Largest allowed time step in seconds. Use this to control the
density of your output. </definition>
<values> float </values>
<suggestion> 1E5 </suggestion>
</parameter>
<parameter name="minIterations">
<definition> Minimum number of Newton iterations of the solver per
time step. At maxTimestep, the Newton solver is not allowed to calculate more
than this number of iterations. </definition>
<values> int </values>
<suggestion> 1 </suggestion>
</parameter>
<parameter name="maxIterations">
<definition> Maximum number of Newton iterations of the solver per
time step. At minTimestep, the Newton solver is not allowed to calculate more
than this number of iterations. </definition>
<values> int </values>
<suggestion> 12 </suggestion>
</parameter>
<parameter name="timestepIncreaseFactor">
<definition> Factor the current time step is multiplied with when increasing
the time step. </definition>
<values> float &gt; 1 </values>
<suggestion> 1.5 </suggestion>
</parameter>
<parameter name="timestepDecreaseFactor">
<definition> Factor the current time step is multiplied with when decreasing
the time step. </definition>
<values> float &lt; 1 </values>
<suggestion> 0.5 </suggestion>
</parameter>
</category>
<category name="parameters">
<parameter name="molecularDiffusion">
<definition> Molecular diffusion of the fluid phase </definition>
<values> float &gt; 0 </values>
<suggestion> 2E-9 </suggestion>
</parameter>
</category>
<category name="numerics">
<parameter name="timestepMethod">
<definition> Numerical scheme to perform time steps in the simulation
</definition>
<values> string </values>
<suggestion> explicit_euler, implicit_euler, heun, shu3, rk4, alex2, fractional, alex3
</suggestion>
</parameter>
<parameter name="courant">
<definition> Courant number for explicit methods. Explicit methods have a
maximum time step according to the CFL condition.
</definition>
<values> float &gt; 0 </values>
<suggestion> 0.8 </suggestion>
</parameter>
</category>
</dorie>
......@@ -18,5 +18,7 @@ USE_MATHJAX = YES
# subdirectory from a directory tree whose root is specified with the INPUT tag.
# exclude versions of the operators that are not used
# EXCLUDE += @top_srcdir@/src/operator_DG_old.hh \
# @top_srcdir@/src/operator_DG_no_upwind.hh
# EXCLUDE += @top_srcdir@/...
# Include documents with other extesions
FILE_PATTERNS += *.dox
\ No newline at end of file
......@@ -45,9 +45,9 @@ namespace Dorie{
virtual BoundaryCondition::Type getBCtype (const Domain& xGlobal, const RF& time) const = 0;
virtual RF getHeadValue (const Domain& xGlobal, const RF& time) const = 0;
virtual RF getDirichletValue (const Domain& xGlobal, const RF& time) const = 0;
virtual RF getFluxValue (const Domain& xGlobal, const RF& time) const = 0;
virtual RF getNeumannValue (const Domain& xGlobal, const RF& time) const = 0;
virtual RF getNextTimeStamp (const RF& time) const = 0;
......@@ -117,16 +117,22 @@ namespace Dorie{
public:
/// Read necessary parameters. Then read out BC data file with rectangular boundary decomposition.
/** \throw Dune::IOError BC File path not specified
/*-----------------------------------------------------------------------*//**
* @brief Read necessary parameters. Then read out BC data file
* with rectangular boundary decomposition.
*
* @throws Dune::IOError BC File path not specified
*
* @param[in] config The configuration
*/
RectangularGrid (const ParameterTree& config)
: BCReadoutInterface<Traits>()
, extensions(config.get<Domain>("grid.extensions"))
, bcFilePath(config.get<std::string>("boundary.file"))
, tStart(config.get<RF>("time.start"))
, tEnd(config.get<RF>("time.end"))
, interpolate(config.get<bool>("boundary.interpolateBCvalues"))
, interpolate(config.get<bool>("boundary.interpolateBCvalues"))
, eps(1e-8)
{
// Check if bcFilePath is specified
......@@ -141,7 +147,7 @@ namespace Dorie{
}
/// Returns the time value of the next BC change. Returns -1 if there is none.
RF getNextTimeStamp (const RF& time) const
RF getNextTimeStamp (const RF& time) const override
{
for(auto it : bcTimeStamps){
if(it>time)
......@@ -167,7 +173,7 @@ namespace Dorie{
* \param time Current time
* \return BoundaryCondition::Type
*/
BoundaryCondition::Type getBCtype (const Domain& xGlobal, const RF& time) const
BoundaryCondition::Type getBCtype (const Domain& xGlobal, const RF& time) const override
{
const typename BC::Side bcSide = getBCside(xGlobal);
......@@ -183,7 +189,7 @@ namespace Dorie{
return bcData[bcSide][bcIndex].type[bcTimeIndex];
}
/// Return Matric Head Value at given global position and time
/// Return Dirihclet Value at given global position and time
/** * Determine boundary side
* * Determine boundary part index
* * Determine time index
......@@ -191,11 +197,11 @@ namespace Dorie{
* * Read out BC data \see Dune::Dorie::RectangularGrid::bcSectionData
* \param xGlobal Global coordinates
* \param time Current time
* \return Matric Head Value
* \return Dirihclet Value
* \throw Dune::Exception Boundary side not determined successfully
* \throw Dune::Exception Dirichlet BC is queried at no-flow boundary
*/
RF getHeadValue (const Domain& xGlobal, const RF& time) const
RF getDirichletValue (const Domain& xGlobal, const RF& time) const override
{
const typename BC::Side bcSide = getBCside(xGlobal);
......@@ -218,7 +224,7 @@ namespace Dorie{
return val0;
}
/// Return Flux Value at given global position and time
/// Return Neumann Value at given global position and time
/** * Determine boundary side
* * Determine boundary part index
* * Determine time index
......@@ -226,10 +232,10 @@ namespace Dorie{
* * Read out BC data \see Dune::Dorie::RectangularGrid::bcSectionData
* \param xGlobal Global coordinates
* \param time Current time
* \return Flux Value
* \return Neumann Value
* \throw Dune::Exception Boundary side not determined successfully
*/
RF getFluxValue (const Domain& xGlobal, const RF& time) const
RF getNeumannValue (const Domain& xGlobal, const RF& time) const override
{
const typename BC::Side bcSide = getBCside(xGlobal);
......
#ifndef DUNE_DORIE_CFL_CONDITION_HH
#define DUNE_DORIE_CFL_CONDITION_HH
#include <dune/common/float_cmp.hh>
#include <limits>
namespace Dune{
namespace Dorie{
/**
* @brief Computes the CFL-condition (\f$\mathcal{CFL}\f$) for a grid
* function representing a velocity field. In particular, the
* return value is the suggested timestep for the lowest CFL
* ratio for each cell:
*
* @f{eqnarray*}{
* \mathcal{CFL} = \min_{T\in\mathcal{T}_h} ||\beta||^{-1}_{[L^\infty(T)]^d}h_T
* @f}
* where \f$\mathcal{T}_h\f$ is a triangulation, \f$\beta\f$ the
* velocity field, and \f$h_T\f$ the diameter of the cell
* \f$T\in \mathcal{T}_h\f$.
* As is usual with the CFL-condition, one can restrict the time
* step (\f$\Delta t\f$) for explicit schemes using a given courant
* number \f$\varrho \le 1\f$:
*
* @f{eqnarray*}{
* \Delta t \le \varrho \cdot \mathcal{CFL}
* @f}
*
* @param[in] gf Grid function representing the veolcity field
* \f$\beta\f$.
* @param[in] intorder The integration order. This value determines the
* quadrature points to evaluate
* \f$||\beta||_{[L^\infty(T)]^d}h_T\f$.
*
* @tparam GF The grid function type.
* @tparam RF The range field. Type to represent CFL values.
*
* @return The CFL-condition (\f$\mathcal{CFL}\f$).
*/
template<class GF, class RF = double>
RF cfl_condition(const GF& gf, unsigned int intorder = 1)
{
using Traits = typename GF::Traits;
using DomainField = typename Traits::DomainFieldType;
using RangeField = typename Traits::RangeFieldType;
const auto& gv = gf.getGridView();
RF suggested_dt = std::numeric_limits<RF>::max();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
auto geo = cell.geometry();
// Check that cell is not a vertex
assert(geo.corners()>1);
DomainField h_T = 0.;
// Cell diameter: maximum distance between all vertices
for (int i = 1; i < geo.corners(); ++i)
for (int j = 0; j < i-1; ++j)
h_T = std::max(h_T,(geo.corner(i)-geo.corner(j)).two_norm());
// Check that h_T is not zero in debug mode, who knows.
assert(Dune::FloatCmp::gt(h_T,0.));
RangeField max_water_flux = 0.;
// Find the maximum velocity in the cell
for (const auto& point : quadratureRule(geo, intorder))
{
typename Traits::RangeType water_flux;
const auto& x = point.position();
gf.evaluate(cell,x,water_flux);
max_water_flux = std::max(max_water_flux,water_flux.two_norm());
}
// Calculate the cfl for this cell
if (Dune::FloatCmp::gt(max_water_flux,0.))
suggested_dt = std::min(suggested_dt,h_T/max_water_flux);
}
// Communicate the lowest cfl number to all procesors
suggested_dt = gv.comm().min(suggested_dt);
return suggested_dt;
}
}
}
#endif
\ No newline at end of file
// These groups should gather information strongly related portions of the code.
// In them we aim to explain to other developers:
// - The rationale of the design
// - How careful someone has to be when modifying it
// - Specific naming or other conventions upheld in the code
/**
@defgroup Common Common
@{
@todo document models
@}
@defgroup Models Models
@defgroup LocalOperators Local operators
@{
Local operators are in many senses the heart of dorie; in them, we transform
finite element basis into a vector of residuals and a mass matrix taking into
account the specific equation and method we want to solve. These
transformations follow the structure of dune-pdelab. In particular,
some care has to be taken regarding the different frames of reference when
working with coordinates and gradients. It can get particularly messy
when working with intersections because one has to express the same
position/gradient with respect to the inside, outside, and face entities,
and sometimes with respect to the macro-grid (e.g. global coordinates).
Therefore, in dorie we try to use the following suffixes convention :
| Convention | Meaning |
|------------|---------|
| `<description>` + `_i` | `<description>` with respect to the **inside** entity. |
| `<description>` + `_o` | `<description>` with respect to the **outside** entity. |
| `<description>` + `_f` | `<description>` with respect to the **face** entity (or intersection). |
| `<description>` + `_g` | `<description>` with respect to the **global** grid. |
| `<description>` + `_n` | `<description>` in **normal** direction. |
Notice that when working with gradients, the suffix for normal direction can
be combined with the first three suffixes (`_n_i`, `_n_o`, and `_n_f`). When
it's not specified, is understood that normal directions are referenced with
respect to the face (`_n` \f$\equiv\f$. `_n_f`).
@todo Update Dune::Dorie::Operator::RichardsDGSpatialOperator to the
local operator convention.
@}
**/
This diff is collapsed.
......@@ -74,7 +74,7 @@ public:
using CON = Dune::PDELab::P0ParallelConstraints;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed, 1>>;
Dune::PDELab::ISTL::VectorBackend<>>;
/// GFS Constraints Container type
using CC = typename Type::template ConstraintsContainer<RF>::Type;
......@@ -105,7 +105,7 @@ public:
using CON = Dune::PDELab::P0ParallelConstraints;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed, 1>>;
Dune::PDELab::ISTL::VectorBackend<>>;
/// GFS Constraints Container type
using CC = typename Type::template ConstraintsContainer<RF>::Type;
......@@ -194,6 +194,16 @@ public:
}
};
/// A wrapper for begin and end time stamps of a time interval
/** \tparam TF Arithmetic data type for time stamps
*/
template<typename TF>
struct TimeInterval {
TF begin;
TF end;
};
/// Traits struct defining basic types based on grid an geometry
template<class GridType, Dune::GeometryType::BasicType GeometryType>
struct BaseTraits
......@@ -201,6 +211,10 @@ struct BaseTraits
static constexpr int dim = GridType::dimension;
static constexpr Dune::GeometryType::BasicType GridGeometryType = GeometryType;
using TF = double;
using TimeField = TF;
using TimeInterval = Dune::Dorie::TimeInterval<TimeField>;
using RF = double;
using RangeField = RF;
using Array = std::vector<RangeField>;
......
......@@ -37,9 +37,8 @@
* \n
* Easy access links for the most important documentation pages:
* \see main Main Program Function
* \see Dune::Dorie::RichardsEquationParameter Class handling Richards Equation Parameter query functions
* \see Dune::Dorie::FlowBoundary Class handling Boundary Condition query functions
* \see Dune::Dorie::FlowSource Class handling Source Term query functions
* \see Dune::Dorie::FlowSource Class handling Source Term query functions
*/
template<typename Traits>
......@@ -183,7 +182,7 @@ int main(int argc, char** argv)
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
else if (gtype == "rectangular"){
......@@ -210,7 +209,7 @@ int main(int argc, char** argv)
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
else{ // no adaptivity
......@@ -233,7 +232,7 @@ int main(int argc, char** argv)
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
}
......@@ -266,7 +265,7 @@ int main(int argc, char** argv)
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
else if (gtype == "rectangular"){
......@@ -293,7 +292,7 @@ int main(int argc, char** argv)
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
else{ // no adaptivity
......@@ -316,7 +315,7 @@ int main(int argc, char** argv)
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
}
......
......@@ -93,7 +93,7 @@ namespace Dune{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getHeadValue(xGlobal, time);
return bcDataHandler->getDirichletValue(xGlobal, time);
}
return 0.0;
}
......@@ -111,7 +111,7 @@ namespace Dune{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getFluxValue(xGlobal, time);
return bcDataHandler->getNeumannValue(xGlobal, time);
}
return 0.0;
}
......
// These groups should gather information strongly related portions of the code.
// In them we aim to explain to other developers:
// - The rationale of the design
// - How careful someone has to be when modifying it
// - Specific naming or other conventions upheld in the code
/**
@defgroup RichardsModel Richards
@{
@ingroup Models
@todo document richards model
@}
**/
......@@ -59,7 +59,7 @@ struct RichardsDGWeights
//! Operator internal BC handling
struct BCType
{
enum Type { vonNeumann, //!< Fixed flux at boundary
enum Type { Neumann, //!< Fixed flux at boundary
dirichlet, //!< Fixed matric head at boundary
none //!< No BC specified
};
......@@ -79,7 +79,7 @@ inline auto read_experimental_operator_settings(
RichardsDGUpwinding::Type upwinding;
RichardsDGWeights::Type weights;
const auto method_str = inifile.get<std::string>("numerics.experimental.method");
const auto method_str = inifile.get<std::string>("dg.experimental.method");
if (method_str == "SIPG")
method = RichardsDGMethod::SIPG;
else if (method_str == "NIPG")
......@@ -92,7 +92,7 @@ inline auto read_experimental_operator_settings(
DUNE_THROW(Dune::IOError, "Unknown DG method '" + method_str + "'!");
const auto upwinding_str
= inifile.get<std::string>("numerics.experimental.upwinding");
= inifile.get<std::string>("dg.experimental.upwinding");
if (upwinding_str == "none")
upwinding = RichardsDGUpwinding::none;
else if (upwinding_str == "semiUpwind")
......@@ -102,7 +102,7 @@ inline auto read_experimental_operator_settings(
else