Commit 0fbd12f7 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'master' into 110-add-scaling-field-input

parents b8964a2b 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.
......@@ -157,3 +157,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.
@}
**/
......@@ -9,21 +9,27 @@
namespace Dune{
namespace Dorie{
/**
/*-------------------------------------------------------------------------*//**
* @brief Base class for simulations.
* @details This class is used to handle simulations with a common interface.
* @author Santiago Ospina.
* @date 2018
* @todo Dune::Dorie::OutputPolicy can handle specific intermediate times
* of the stepping.
* @ingroup Models
*/
class SimulationBase
{
public:
/**
* @brief Constructs the SimulationBase.
*
* @param[in] output_policy The output policy.
* Defaults to OutputPolicy::EndOfStep.
* @param[in] adapt_policy The adapt policy.
* Defaults to AdaptiviyPolicy::None.
*/
/*-------------------------------------------------------------------------*//**
* @brief Constructs the SimulationBase.
*
* @param[in] output_policy The output policy. Defaults to
* OutputPolicy::EndOfStep.
* @param[in] adapt_policy The adapt policy. Defaults to
* AdaptiviyPolicy::None.
*/
SimulationBase(
OutputPolicy output_policy=OutputPolicy::EndOfStep,
AdaptivityPolicy adapt_policy=AdaptivityPolicy::None
......@@ -32,14 +38,14 @@ public:
, _adapt_policy(adapt_policy)
{}
/**
/*-----------------------------------------------------------------------*//**
* @brief Sets the output policy.
*
* @param[in] output_policy The output policy,
*/
void set_policy(OutputPolicy output_policy) {_output_policy = output_policy;}
/**
/*-----------------------------------------------------------------------*//**
* @brief Sets the adaptivity policy.
*
* @param[in] adapt_policy The adaptivity policy.
......@@ -49,7 +55,7 @@ public:
_adapt_policy = adapt_policy;
}
/**
/*-----------------------------------------------------------------------*//**
* @brief Returns the current output policy.
*
* @return The current output policy.
......@@ -57,7 +63,7 @@ public:
OutputPolicy output_policy() const {return _output_policy;}
/**
/*-----------------------------------------------------------------------*//**
* @brief Returns the current adaptivity policy.
*
* @return The current adaptivity policy.
......@@ -65,7 +71,7 @@ public:
AdaptivityPolicy adaptivity_policy() const {return _adapt_policy;}
/**
/*-----------------------------------------------------------------------*//**
* @brief Writes the data.
*/
void virtual write_data() const
......@@ -74,7 +80,7 @@ public:
DUNE_THROW(Dune::IOError,"Simulation can't write data!");
}
/**
/*-----------------------------------------------------------------------*//**
* @brief Mark the grid in order to improve the current model.
*/
virtual void mark_grid()
......@@ -83,14 +89,14 @@ public:
"Simulation can't mark the grid for adaptivity!");
}
/**
/*-----------------------------------------------------------------------*//**
* @brief Operations before adaptation of the grid
*/
virtual void pre_adapt_grid() {};
/**
* @brief Adapt the grid together it every dependency of the
* grid (e.g. solution vector and grid function spaces).
/*-----------------------------------------------------------------------*//**
* @brief Adapt the grid together it every dependency of the grid (e.g.
* solution vector and grid function spaces).
*/
virtual void adapt_grid()
{
......@@ -101,47 +107,48 @@ public:
<< "policy None. Method adapt() must not be called!");
}
/**
/*-----------------------------------------------------------------------*//**
* @brief Operations after adaptation of the grid
*/
virtual void post_adapt_grid() {};
/**
/*-----------------------------------------------------------------------*//**
* @brief Method that provides the begin time of the model.
*
* @return Begin time of the model.
*/
virtual double begin_time() const = 0;
/**
/*-----------------------------------------------------------------------*//**
* @brief Method that provides the end time of the model.
*
* @return End time of the model.
*/
virtual double end_time() const = 0;
/**
/*-----------------------------------------------------------------------*//**
* @brief Method that provides the current time of the model.
*
* @return Current time of the model.
*/
virtual double current_time() const = 0;
/**
/*-----------------------------------------------------------------------*//**
* @brief Suggest a time step to the model.
*
* @param[in] suggestion for the internal time step of the model.
* @param[in] dt { parameter_description }
* @param[in] suggestion for the internal time step of the model.
*/
virtual void suggest_timestep(double dt) = 0;
/**
/*-----------------------------------------------------------------------*//**
* @brief Performs one steps in direction to end_time(). The time-step
* should never result on a bigger step than the one suggested
* in suggest_timestep().
* should never result on a bigger step than the one suggested in
* suggest_timestep().
*/
virtual void step() = 0;
/**
/*-----------------------------------------------------------------------*//**
* @brief Runs the model performing steps until current_time() equals
* end_time()
*/
......@@ -170,7 +177,109 @@ private:
AdaptivityPolicy _adapt_policy;
};
/**
@{
@class SimulationBase
## Base class for simulations
### Stepping:
The method `step()` must be an abstract method that always has to be
implemented by the derived class. The step must succeed always, otherwise,