add documentation for transport models

parent 2f01b72f
add_subdirectory("groups")
# shortcut for creating the Doxyfile.in and Doxyfile
add_doxygen_target()
......@@ -3,7 +3,6 @@
# Location of source files
INPUT += @top_srcdir@/dune/dorie/
INPUT += @top_srcdir@/doc/doxygen/groups
# extract private class members
EXTRACT_PRIVATE = YES
......@@ -20,3 +19,6 @@ USE_MATHJAX = YES
# exclude versions of the operators that are not used
# EXCLUDE += @top_srcdir@/...
# Include documents with other extesions
FILE_PATTERNS += *.dox
\ No newline at end of file
/**
@defgroup Common Common
@{
@todo document models
@}
**/
\ No newline at end of file
/**
@defgroup Models Models
@{
@todo document models
@}
@defgroup RichardsModel Richards
@{
@ingroup Models
@todo document richards model
@}
@defgroup TransportModel Transport
@{
@ingroup Models
@todo document transport model
@}
**/
\ No newline at end of file
/** @defgroup LocalOperators Local operators
// These groups should gather information of code that belong to fairly similar
// categories. In them we aim to explain the code: why is done in that specific
// way, how carful someone has to be to modify it, specific name conventions, etc.
// These descriptions are not for regular users but for future developers (maybe
// the same developer in some months/years later).
/**
@defgroup Common Common
@{
@todo document models
@}
@defgroup Models Models
@defgroup RichardsModel Richards
@{
@ingroup Models
@todo document richards model
@}
@defgroup TransportModel Transport
@{
@ingroup Models
@todo document transport model
@}
@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. In particular,
some care has been taken regarding the different frames of reference when
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
particularly 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 grid (e.g. global coordinates). Therefore,
in dorie we try to use the following suffixes convention:
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. |
| `<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 it is with respect to the face
(`_n` \f$\equiv\f$. `_n_f`).
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 convention.
@todo Update Dune::Dorie::Operator::RichardsDGSpatialOperator to the
local operator convention.
@}
**/
**/
\ No newline at end of file
......@@ -10,20 +10,130 @@ namespace Dune{
namespace Dorie{
/**
@addtogroup Models
@{
## 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,
it should throw a `SimulationStepException`. It means that time adaptivity
must be done internally so that the time step succeeds.
### Writing data:
The method `step()` must be able to decide whether to write data or not.
(Say that in one case one wants that transport writes that every step and in
another case one wants that it writes at the end of the Richards step. The
first case can't be covered writing data in the `run()` method of the
coupled simulation.) Then, enum class called `OutputPolicy` must decide
whether to call the method `write_data()`. This policy can be set and be
extracted by the methods `set_policy(OutputPolicy)` and `output_policy()`
respectively.
The method `write_data()` should be in charge of writing the data.
In principle, this method should be only used following the `OutputPolicy`.
Usually called at the end of each step and at the beginning of the
simulation. Since the `step()` method is in charge of it, it is assumed that
this method writes the *current* state of the model (e.g. intermediate
results of the coupled simulation for the transport state).
### Adaptivity:
For the purpose of having a generic `run()` method in this class, the sketch
of the adaptivity must be provided here. First, adaptivity must be never
done inside the method `step()`. Additionally, if the model can perform
adaptivity, it must provide the method `mark_grid()`, otherwise, it must
throw a `NotImplemented` exception. For completion, there must be a method
called `adapt_grid()` that pack all the known degrees of freedom and, the
grid function spaces and perform the adaptivity.
Similar to the idea of `OutputPolicy` there must be an `AdaptivityPolicy`
set by the methods `set_policy(AdaptivityPolicy)` and `adaptivity_policy()`
respectively which will be set by default to `AdaptivityPolicy::None`.
Having this separation between marking the grid and adaptivity of solutions
leave coupled models to decide which model will mark the grid inside its
method of `mark_grid()` and still handle the adaptivity implementing
`adapt_grid()` separately.
### Time control:
A minimal interface is only providing `begin_time()`, `end_time()` and
`current_time()`.
In order to keep solution vectors minimal at the time of adaptivity, DORiE
expect the several models synchronized in time. Therefore, a good and still
minimal interface is only suggesting the timestep to the model via
`suggest_timestep(double dt)`. Such suggestion has to be taken into account
within the model, and the model must ensure that the timestep is done in
`step()` is never bigger than the suggested one so that coupled models can
achieve synchronization. Naturally, models have to ensure that
`begin_time()` \f$=\f$ `current_time()` at when they are constructed.
### Local Operator Setup:
Since adaptivity can modify the complete setup of the simulation, and since
the method `adapt_grid()` is not necessarily called by the model.
There is a need for the method `post_adapt_grid()` which is in charge of
setup everything again to the new grid. Same argument holds for
`pre_adapt_grid()`.
### Run:
With all the definitions given before, a general algorithm to run the
program can be given by following code:
```
virtual void run()
{
if (output_policy() != OutputPolicy::None)
write_data();
while( current_time() < end_time() )
{
step();
if (adaptivity_policy() != AdaptivityPolicy::None)
if ( current_time() < end_time() )
{
mark_grid();
pre_adapt_grid();
adapt_grid();
post_adapt_grid();
}
}
}
```
@}
*/
/*-------------------------------------------------------------------------*//**
* @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.
* @copyright MIT License.
* @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 +142,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 +159,7 @@ public:
_adapt_policy = adapt_policy;
}
/**
/*-----------------------------------------------------------------------*//**
* @brief Returns the current output policy.
*
* @return The current output policy.
......@@ -57,7 +167,7 @@ public:
OutputPolicy output_policy() const {return _output_policy;}
/**
/*-----------------------------------------------------------------------*//**
* @brief Returns the current adaptivity policy.
*
* @return The current adaptivity policy.
......@@ -65,7 +175,7 @@ public:
AdaptivityPolicy adaptivity_policy() const {return _adapt_policy;}
/**
/*-----------------------------------------------------------------------*//**
* @brief Writes the data.
*/
void virtual write_data() const
......@@ -74,7 +184,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 +193,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 +211,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()
*/
......
......@@ -35,9 +35,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>
......
......@@ -21,7 +21,7 @@ namespace Operator {
/*-------------------------------------------------------------------------*//**
* @brief Spatial local operator for the transport equation in unsaturated
* media in a finite volume scheme.
*
* @details It solves the spatial part of the transport equation:
* @f{eqnarray*}{
* \partial_t[\theta C_w] +
* \nabla\cdot [\textbf{j}_w C_w] -
......@@ -33,12 +33,14 @@ namespace Operator {
* C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
* \Gamma_D
* @f}
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup LocalOperators
*
* @todo Use diffusion coefficient grid function.
* @todo Implement outflow boundary condition.
*
* @ingroup LocalOperators
*
* @tparam Boundary Type of the class providing boundary
* conditions
......@@ -416,8 +418,25 @@ private:
/**
* @brief Temporal local operator for the transport equation in unsaturated
* media in a finite volume scheme.
*
* @details It solves the temporal part of the transport equation:
* @f{eqnarray*}{
* \partial_t[\theta C_w] +
* \nabla\cdot [\textbf{j}_w C_w] -
* \nabla [\theta \mathsf{D}_{eff}\nabla C_w]=0 &\qquad \text{in }
* \Omega\\
* C_w = g &\qquad \text{on } \Gamma_D
* \subseteq\partial\Omega\\
* \nabla C_w \cdot \textbf{n} = \textbf{j}_{\scriptscriptstyle
* C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
* \Gamma_D
* @f}
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup LocalOperators
*
* @bug The water content is not being updated by the set time to the
* right state of the grid function.
*
* @tparam GridFunctionWaterContent Type of a grid function which provides
* the water content
......@@ -501,6 +520,7 @@ public:
void setTime (RF t)
{
_time = t;
_gf_water_content->setTime(t);
}
private:
......
......@@ -11,103 +11,124 @@
namespace Dune{
namespace Dorie{
/// Boundary type and condition value queries
/** This class containts functions that return the type of boundary conditions
* and the values of the boundary condition parameters. Both types of queries
* can be time dependent.
*/
template<typename Traits>
class SoluteBoundary
{
private:
/*-------------------------------------------------------------------------*//**
* @brief Boundary type and condition value queries for solute.
* @details This class containts functions that return the type of boundary
* conditions and the values of the boundary condition parameters.
* Both types of queries can be time dependent.
* @author Dion Häfner
* @date 2016-2018
* @copyright MIT License.
* @ingroup TransportModel
* @see Same as Dune::Dorie::FlowBoundary.
*
* @tparam Traits The Dune::Dorie::BaseTraits that defines basic data
* types.
*/
template<typename Traits>
class SoluteBoundary
{
private:
typedef typename Traits::RangeField RF;
typedef typename Traits::Domain Domain;
typedef typename Traits::IntersectionDomain ID;
typedef typename Traits::Intersection Intersection;
typedef typename Traits::RangeField RF;
typedef typename Traits::Domain Domain;
typedef typename Traits::IntersectionDomain ID;
typedef typename Traits::Intersection Intersection;
enum {dim = Traits::dim};
enum {dim = Traits::dim};
//! Object for handling the BC data file and function queries
std::unique_ptr<BCReadoutInterface<Traits>> bcDataHandler;
//! Object for handling the BC data file and function queries
std::unique_ptr<BCReadoutInterface<Traits>> bcDataHandler;
public:
public:
/// Read BC file type and create data handling object
/** \see Dune::Dorie::BCReadoutInterface
* \throw Dune::IOError BC File Type not supported
*/
SoluteBoundary(const Dune::ParameterTree& config)
{
std::string bcFileType = config.get<std::string>("boundary.fileType");
if(bcFileType == "rectangularGrid") {
bcDataHandler = std::make_unique<RectangularGrid<Traits>>(config);
}
else
DUNE_THROW(IOError,"unknown bcFileType specified!");
}
/*-----------------------------------------------------------------------*//**
* @brief Read BC file type and create data handling object
* @see Dune::Dorie::BCReadoutInterface
*
* @param[in] config The configuration parameter three.
* @throws Dune::IOError BC File Type not supported.
*/
SoluteBoundary(const Dune::ParameterTree& config)
{
std::string bcFileType = config.get<std::string>("boundary.fileType");
if(bcFileType == "rectangularGrid") {
bcDataHandler = std::make_unique<RectangularGrid<Traits>>(config);
}
else
DUNE_THROW(IOError,"unknown bcFileType specified!");
}
/// Return next time at which any boundary condition changes. Return -1.0 by default
/** \param time Current time
*/
RF getNextTimeStamp (const RF& time) const {
if(bcDataHandler)
return bcDataHandler->getNextTimeStamp(time);
return -1.0;
}
/*-----------------------------------------------------------------------*//**
* Return next time at which any boundary condition changes. Return -1.0 by
* default
*
* @param time Current time
*
* @return The next time stamp.
*/
RF getNextTimeStamp (const RF& time) const {
if(bcDataHandler)
return bcDataHandler->getNextTimeStamp(time);
return -1.0;
}
/// Return Boundary Condition Type at certain position and time
/** \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xGlobal Global coordinate
* \return Boundary condition type enumerator
*/
BoundaryCondition::Type bc (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
return bcDataHandler->getBCtype(xGlobal, time);
}
return BoundaryCondition::Other; // unknown
}
/*-----------------------------------------------------------------------*//**
* @brief Return Boundary Condition Type at certain position and time
*
* @param is Intersection enitity.
* @param x Position in local entity coordinates.
* @param time Time value.
*
* @return Boundary Condition type enumerator (BoundaryCondition::Type).
*/
BoundaryCondition::Type bc (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
return bcDataHandler->getBCtype(xGlobal, time);
}
return BoundaryCondition::Other; // unknown
}
/**
* @brief Dirichlet boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of matric head
*/
RF g (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getDirichletValue(xGlobal, time);
}
return 0.0;
}
/*-----------------------------------------------------------------------*//**
* @brief Dirichlet boundary condition at certain position and time
*
* @param is Intersection enitity
* @param x Position in local entity coordinates
* @param time Time value
*
* @return Value of matric head
*/
RF g (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getDirichletValue(xGlobal, time);
}
return 0.0;
}
/**
* @brief Neumann boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of flux
*/
RF j (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getNeumannValue(xGlobal, time);
}
return 0.0;
}
};
/*-----------------------------------------------------------------------*//**
* @brief Neumann boundary condition at certain position and time
*
* @param is Intersection enitity
* @param x Position in local entity coordinates
* @param time Time value
*
* @return Value of flux
*/
RF j (const Intersection& is, const ID& x, const RF& time) const
{
const Domain xGlobal = is.geometry().global(x);
if(is.boundary()) {
if(bcDataHandler)
return bcDataHandler->getNeumannValue(xGlobal, time);
}
return 0.0;
}
};
}
}
......
......@@ -6,9 +6,41 @@
namespace Dune {
namespace Dorie {
template<typename T>
using SoluteInitial = Dune::PDELab::ConstGridFunction<typename T::GridView,typename T::RangeField>;
#ifndef DOXYGEN /* Actual code */
template<typename T>
using SoluteInitial = Dune::PDELab::ConstGridFunction<
typename T::GridView,
typename T::RangeField>;
#else /* Documentation of the code */
/**
* @brief Class for initial condition of the solute concentration.
* @remark This class is equivalent to Dune::PDELab::ConstGridFunction.
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup TransportModel
* @see Same as Dune::Dorie::FlowInitial.
*
* @todo Implement useful shape functions (e.g.
* <a href="https://en.wikipedia.org/wiki/Kernel_(statistics)">kernels</a>
* ) that allow more realistic boundary conditions for the transport
* simulation.
*
* @tparam Traits The Dune::Dorie::BaseTraits that defines basic data
* types.
*/
template<typename T>
class SoluteInitial
: public Dune::PDELab::ConstGridFunction<
typename T::GridView,
typename T::RangeField>
{};
#endif
}
}
......
......@@ -31,11 +31,15 @@ namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Traits for the class TransportSimulation. In particular, this
* class extends BaseTraits to be used for the DG implementation of
* the Transport simulation
*
* @brief Traits for the class TransportSimulation.
* @details This class extends BaseTraits to be used for the
* implementation of the Transport simulation.
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup TransportModel
*
* @todo Allow higher order methods.
*