Commit 00e61325 authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos

Merge branch '111-add-dg-local-operator-to-the-simulationtransport-class' into 'master'

Resolve "Add DG local operator to the SimulationTransport class"

Closes #111

See merge request !112
parents 52fa9ee5 8235700d
......@@ -124,6 +124,9 @@ prep:update-dune-clang:
# --- Build jobs ---
build:system-tests: &build-tests
stage: build
# Reduce concurrent jobs due to RAM limitations
variables:
MAKE_FLAGS: -j 2
script:
- CMAKE_FLAGS="$CMAKE_FLAGS -DCOVERAGE_REPORT=On"
$DUNECONTROL --only=dorie configure
......
......@@ -35,6 +35,7 @@
* Steady state initial condition in Richards model !176
* Changes to config file parameters listed per version in user docs !175
* Control negative transport solution by a check policy !181
* DG solver for solute transport model !112
### Changed
* Data structures for storing and accessing parameter information !55
......
......@@ -9,16 +9,16 @@ set(DORIE_MAX_RORDER_2 6)
set(DORIE_MAX_RORDER_3 6)
# Maximum polynomial orders of transport model for available targets
set(DORIE_MAX_TORDER_2 0)
set(DORIE_MAX_TORDER_3 0)
set(DORIE_MAX_TORDER_2 3)
set(DORIE_MAX_TORDER_3 3)
# Maximum polynomial orders of Richards model for default targets
set(DORIE_MAX_DEFAULT_RORDER_2 3)
set(DORIE_MAX_DEFAULT_RORDER_3 1)
# Maximum polynomial orders of transport model for default targets
set(DORIE_MAX_DEFAULT_TORDER_2 0)
set(DORIE_MAX_DEFAULT_TORDER_3 0)
set(DORIE_MAX_DEFAULT_TORDER_2 3)
set(DORIE_MAX_DEFAULT_TORDER_3 1)
#
# .. cmake_function:: dorie_compile_instance
......@@ -107,11 +107,13 @@ function(dorie_compile_instance)
endif()
# Add the executable to the default targets
if ((ARGS_RORDER LESS_EQUAL DORIE_MAX_DEFAULT_RORDER_${ARGS_DIMENSION})
AND ((NOT ARGS_TORDER)
OR ARGS_TORDER LESS_EQUAL DORIE_MAX_DEFAULT_TORDER_${ARGS_DIMENSION})
)
add_dependencies(${ARGS_MODEL} ${exe_name})
if (ARGS_RORDER LESS_EQUAL DORIE_MAX_DEFAULT_RORDER_${ARGS_DIMENSION})
if (ARGS_MODEL STREQUAL "richards")
add_dependencies(${ARGS_MODEL} ${exe_name})
elseif ((ARGS_TORDER LESS_EQUAL DORIE_MAX_DEFAULT_TORDER_${ARGS_DIMENSION})
AND (ARGS_RORDER EQUAL ARGS_TORDER))
add_dependencies(${ARGS_MODEL} ${exe_name})
endif()
endif()
# set compile definitions
......
......@@ -232,7 +232,7 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
* **semiUpwind:** Apply upwinding to conductivity factor (only).
* **fullUpwind:** Apply upwinding to conductivity.
**Not recommended for DG**.
**Not recommended for DG, but for FV**.
</definition>
<values> none, semiUpwind, fullUpwind </values>
<suggestion> none </suggestion>
......
......@@ -172,22 +172,6 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
<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>
......@@ -215,7 +199,7 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
<category name="numerics">
<parameter name="timestepMethod">
<definition> Numerical scheme to perform time steps in the simulation.
``alex2`` and ``implicit_euler`` are implicit methods.
``alex2`` and ``implicit_euler`` are implicit methods.
``explicit_euler`` is a explicit method.
</definition>
<values> explicit_euler, implicit_euler, alex2 </values>
......@@ -234,15 +218,59 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
</parameter>
<parameter name="FEorder">
<definition> Order of the finite element method used. </definition>
<suggestion> 0 </suggestion>
<values> 0 </values>
<definition> Polynomial order of the DG method used. Setting this value
to 0 (zero) selects the finite volume (FV) solver (only compatible to
structured rectangular grids).
</definition>
<suggestion> 1 </suggestion>
<values> 0, 1, 2, 3 </values>
<comment> Select '0' for the finite volume (FV) solver </comment>
</parameter>
<parameter name="upwinding">
<definition> Enable upwinding for the convective flux term.
Upwinding leads to quasi-optimal error estimates of convective flux
problems and is therefore **recommended for DG schemes**. For FV,
upwinding is always enabled.
</definition>
<values> true, false </values>
<suggestion> true </suggestion>
</parameter>
<parameter name="DGMethod">
<definition> DG discretization method for skeleton terms.
* **SIPG:** Symmetric Interior Penalty
* **NIPG:** Non-Symmetric Interior Penalty
* **OBB:** Oden, Babuska, Baumann: no penalty term
* **IIP:** Incomplete Interior Penalty: no symmetry term
</definition>
<values> SIPG, NIPG, OBB, IIP </values>
<suggestion> SIPG </suggestion>
</parameter>
<parameter name="DGWeights">
<definition> Apply harmonic weighting to skeleton term contributions
in DG.
</definition>
<values> true, false </values>
<suggestion> true </suggestion>
</parameter>
<parameter name="penaltyFactor">
<definition> Penalty factor to be used in the Discontinuous Galerkin scheme
</definition>
<values> float </values>
<suggestion> 10 </suggestion>
</parameter>
</category>
<category name="fluxReconstruction">
<parameter name="enable">
<definition> Apply the flux reconstruction method to the solved solute
<definition> Apply the flux reconstruction method to the solved solute
concentration and obtain conservative fluxes.
</definition>
<suggestion> true </suggestion>
......@@ -251,8 +279,8 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
<parameter name="checkJumps">
<definition> Check that flux reconstruction engine is creating conforming
normal fluxes up to ``jumpTol``. ProTip: Setting warnings together
with a very low tolerance will let you track the changes on the
normal fluxes up to ``jumpTol``. ProTip: Setting warnings together
with a very low tolerance will let you track the changes on the
quality of the flux reconstruction.
</definition>
<suggestion> none </suggestion>
......@@ -260,8 +288,8 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
<definition> Whenever ``checkJumps`` is activated, it check that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
......@@ -301,15 +329,15 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
</category>
<category name="solverParameters">
<parameter name="reduction">
<definition> Required relative defect reduction.
<parameter name="reduction" hidden="true">
<definition> Required relative defect reduction in the linear solver.
Reduce this value to increase precision.
</definition>
<suggestion> 1E-6 </suggestion>
<suggestion> 1E-8 </suggestion>
<values> float </values>
</parameter>
<parameter name="minDefect">
<parameter name="min_defect" hidden="true">
<definition> Minimum absolute defect at which linear solver stops.
Reduce this value to increase precision.
</definition>
......
......@@ -100,13 +100,17 @@ are the available combinations of options:
| +-------------------+-----------------------------+------------------------------+----------------+
| | 3 | <= 1 | ✗ | ✓* |
+----------------------+-------------------+-----------------------------+------------------------------+----------------+
|``richards+transport``| 2 | <= 3 | 0 | ✗ |
|``richards+transport``| 2 | <= 3 | same as ``richards`` | ✓*\ :sup:`†` |
| +-------------------+-----------------------------+------------------------------+----------------+
| | 3 | <= 1 | 0 | ✗ |
| | 3 | <= 1 | same as ``richards`` | ✓*\ :sup:`†` |
+----------------------+-------------------+-----------------------------+------------------------------+----------------+
* *: Only for ``richards.numerics.FEorder`` > 0. Finite volume solvers do not
support unstructured grids.
* :sup:`†`: Grid refinement on *rectangular* grids produces hanging nodes which
are currently not considered in the
:ref:`flux reconstruction <man-flux_reconstruction>`. This can result in
erroneous results computed by the transport solver.
.. _DUNE: https://www.dune-project.org/
.. _DUNE-PDELab: https://www.dune-project.org/modules/dune-pdelab/
......@@ -15,6 +15,8 @@ because the basis function gradients are always zero. In this case, the only
option for evaluating flux data is flux reconstruction which has to be enabled
in the :doc:`config file <config-file>`.
.. _man-flux_reconstruction:
Flux reconstruction
-------------------
The flux reconstruction is a projection of the fluxes used in the Discontinuous Galerkin method into a vector field function. Using correct elements, this procedure can ensure that fluxes in normal direction to the element are *equivalent* to those computed by the Discontinuous Galerkin method, and most importantly, it can also ensure the continuity of them. Hence, the resulting vector field is useful to compute other problems that rely on the fluxes of the water (i.e. solute transport).
......@@ -50,4 +52,8 @@ Legend:
Usage
-----
To activate/deactivate flux reconstruction use the keyword ``richards.numeric.fluxReconstruction = true/false`` in the :doc:`config file <config-file>`.
To activate/deactivate flux reconstruction use the keyword
``numeric.fluxReconstruction = true/false`` in the
:doc:`config file <config-file>`. Flux reconstruction of the water flow is
always enabled in a coupled Richards and transport simulation because it is
required by the discretization of the convection-diffusion equation.
#ifndef DUNE_DORIE_RAVIART_THOMAS_PROJECTION_HH
#define DUNE_DORIE_RAVIART_THOMAS_PROJECTION_HH
// hard-coded bool helpers
/*-------------------------------------------------------------------------*//**
* @brief Helper for the simplex raviart thomas flux reconstruction.
* @details This helper returns a bool that says whether the flux
* reconstruction method in simplices is available for a given
* dimension and order.
*
* @tparam dim Grid dimension (2 or 3)
* @tparam order Order of the flux reconstruction.
*/
#define DORIE_SUPPORT_RT_SIMPLEX(dim, order) (order <= 2)
/*-------------------------------------------------------------------------*//**
* @brief Helper for the cube raviart thomas flux reconstruction.
* @details This helper returns a bool that says whether the flux
* reconstruction method in cubes is available for a given
* dimension and order.
*
* @tparam dim Grid dimension (2 or 3)
* @tparam order Order of the flux reconstruction.
*/
#define DORIE_SUPPORT_RT_CUBE(dim, order) (dim==2) ? (order <= 2) : (order <= 1)
#include <vector>
......
......@@ -240,6 +240,9 @@
@defgroup LocalOperators Local Operators
@{
@brief
Operators for assembling residuals and jacobians at element-local level
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
......
......@@ -46,7 +46,7 @@ using Geo = Dune::GeometryType::BasicType;
#warning Coupling for cubes is not available for this order and dimension. Case (YaspGrid) will be ignored
#endif
#if DORIE_TORDER > 0
#if (DORIE_TORDER > 0) && (DORIE_RORDER > 0)
#if DORIE_SUPPORT_RT_CUBE(DORIE_DIM, DORIE_RT_ORDER)
template class ModelRichardsTransportCoupling<ModelRichardsTransportCouplingTraits<BaseTraits<UGGrid<DORIE_DIM>, Geo::cube>,
DORIE_RORDER,
......
......@@ -9,8 +9,11 @@
// Checks for custom setup
#ifndef DORIE_TORDER
#error DORIE_TORDER macro must be defined for Transport factory builds
#elif DORIE_TORDER != 0
#error "Transport FV solver only supports DORIE_TORDER=0"
#else
static_assert(DORIE_TORDER >= 0,
"Negative polynomial orders do not make sense");
static_assert(DORIE_TORDER <= 6,
"DORiE only supports polynomial orders up to 6");
#endif
#ifdef DORIE_RT_ORDER
......@@ -129,26 +132,55 @@ public:
if (grid_mode == GridMode::gmsh)
{
raise_error_fv_grid();
}
else // GridMode::rectangular
{
if (grid_adaptive) {
#if (DORIE_RORDER == 0 || DORIE_TORDER == 0)
raise_error_fv_grid();
}
// grid not adaptive
else {
using GridType = StructuredGrid<dim>;
#elif DORIE_SUPPORT_RT_SIMPLEX(DORIE_DIM, DORIE_RT_ORDER)
using GridType = UnstructuredGrid<dim>;
GridCreator<GridType> grid_creator(config, helper);
#if DORIE_SUPPORT_RT_CUBE(DORIE_DIM, DORIE_RT_ORDER)
mod = std::make_shared<ModTransportCube<GridType, DORIE_RORDER, DORIE_TORDER>>(
mod = std::make_shared<ModTransportSimplex<GridType,
DORIE_RORDER,
DORIE_TORDER>>(
config_richards, config_transport, grid_creator, helper
);
#else
#else
raise_error_flux_reconstruction();
#endif
}
else // GridMode::rectangular
{
#if DORIE_SUPPORT_RT_CUBE(DORIE_DIM, DORIE_RT_ORDER)
if (grid_adaptive) {
#if (DORIE_RORDER == 0 || DORIE_TORDER == 0)
raise_error_fv_grid();
#else
using GridType = UnstructuredGrid<dim>;
GridCreator<GridType> grid_creator(config, helper);
mod = std::make_shared<ModTransportSimplex<
GridType, DORIE_RORDER, DORIE_TORDER>>
(
config_richards,
config_transport,
grid_creator,
helper
);
#endif
}
// grid not adaptive
else {
using GridType = StructuredGrid<dim>;
GridCreator<GridType> grid_creator(config, helper);
mod = std::make_shared<ModTransportCube<GridType,
DORIE_RORDER,
DORIE_TORDER>>(
config_richards, config_transport, grid_creator, helper
);
}
#else
raise_error_flux_reconstruction();
#endif
}
#endif
}
// set policy and return
......
#ifndef DUNE_DORIE_LOP_CFG_HH
#define DUNE_DORIE_LOP_CFG_HH
#include <tuple>
#include <string>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune {
namespace Dorie {
namespace Operator {
/**
* \addtogroup Models
* \{
*/
//! Interior penalty type
enum class DGMethod
{
NIPG, //!< Non-symmetric interior penalty
SIPG, //!< Symmetric interior penalty (default)
OBB, //!< Oden, Babuska, Baumann (no penalty term)
IIP //!< Incomplete interior penalty (no symmetry term)
};
//! Gradient weighting
enum class DGWeights
{
weightsOn, //!< harmonic weighting of flux contributions
weightsOff //!< no weighting of flux contributions
};
//! Upwinding type
enum class Upwinding
{
none, //!< No upwinding
fullUpwind //!< Upwind
};
/// Read upwinding operator settings from a parameter file
/* \return Upwinding option
*/
inline Operator::Upwinding read_operator_upwinding_settings(
const Dune::ParameterTree& inifile,
const std::shared_ptr<spdlog::logger>& log)
{
const auto upwinding_str
= inifile.get<std::string>("numerics.upwinding");
log->debug(" Upwinding: {}", upwinding_str);
if (upwinding_str == "false")
return Operator::Upwinding::none;
else if (upwinding_str == "true")
return Operator::Upwinding::fullUpwind;
else {
log->error("Unknown upwinding setting: {}", upwinding_str);
DUNE_THROW(Dune::IOError, "Unknown upwinding setting");
}
}
/// Read DG weights operator settings from a parameter file
/* \return DG weights option
*/
inline Operator::DGWeights read_operator_dg_weights_settings(
const Dune::ParameterTree& inifile,
const std::shared_ptr<spdlog::logger>& log)
{
const auto weights_str = inifile.get<bool>("numerics.DGWeights");
log->debug(" Harmonic weights for interface fluxes (DG): {}",
weights_str);
if (weights_str)
return Operator::DGWeights::weightsOn;
else
return Operator::DGWeights::weightsOff;
}
/// Read DG method operator settings from a parameter file
/* \return DG method option
*/
inline Operator::DGMethod read_operator_dg_method_settings(
const Dune::ParameterTree& inifile,
const std::shared_ptr<spdlog::logger>& log)
{
const auto method_str = inifile.get<std::string>("numerics.DGMethod");
log->debug(" DG method: {}", method_str);
if (method_str == "SIPG")
return Operator::DGMethod::SIPG;
else if (method_str == "NIPG")
return Operator::DGMethod::NIPG;
else if (method_str == "OBB")
return Operator::DGMethod::OBB;
else if (method_str == "IIP")
return Operator::DGMethod::IIP;
else {
log->error("Unknown DG method: {}", method_str);
DUNE_THROW(Dune::IOError, "Unknown DG method");
}
}
/**
* \}
* // group Models
*/
} // namespace Operator
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_LOP_CFG_HH
......@@ -15,7 +15,6 @@
@defgroup BuildProcess Model Build Process
@{
@ingroup Models
@brief How DORiE libraries and executables are built and linked
## Overview
......@@ -59,8 +58,8 @@
| ---------- | ----- | -------- | -------- |
| `richards` | 2 | <= 3 | none |
| ^ | 3 | <= 1 | none |
| `transport` | 2 | <= 3 | 0 |
| ^ | 3 | <= 1 | 0 |
| `transport` | 2 | <= 3 | =`RORDER` |
| ^ | 3 | <= 1 | =`RORDER` |
Available pre-defined targets:
......@@ -68,8 +67,8 @@
| ---------- | ----- | -------- | -------- |
| `richards` | 2 | <= 6 | none |
| ^ | 3 | ^ | none |
| `transport` | 2 | <= 4 | 0 |
| ^ | 3 | ^ | 0 |
| `transport` | 2 | <= 6 | <= 3 |
| ^ | 3 | ^ | <= 3 |
If the provided instance are not sufficient for your needs, use the CMake
function `dorie_compile_instance` to register additional model instance
......
#ifndef DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#define DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#ifndef DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
#define DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
#include <iostream>
#include <string>
......@@ -22,24 +22,25 @@
namespace Dune{
namespace Dorie{
/// Adaptivity base class. Does nothing.
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U, int order>
class AdaptivityHandlerBase
class RichardsAdaptivityHandler
{
private:
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
public:
AdaptivityHandlerBase () = default;
virtual ~AdaptivityHandlerBase () = default;
RichardsAdaptivityHandler () = default;
virtual ~RichardsAdaptivityHandler () = default;
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
/// Do nothing.
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
};
/// Specialized AdaptivityHandlerBase class providing functions and members for adaptive grid refinement
......@@ -49,7 +50,7 @@ template<typename Traits,
typename Boundary,
typename U,
int order>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>
class WaterFluxAdaptivityHandler : public RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>
{
private:
using RF = typename Traits::RF;
......@@ -61,7 +62,7 @@ private:
/// Adaptivity GFS
using AGFS = typename AGFSHelper::Type;
/// Error estimator local operator
using ESTLOP = Dune::Dorie::Operator::FluxErrorEstimator<Traits,Param,Boundary>;
using ESTLOP = Dune::Dorie::Operator::WaterFluxErrorEstimator<Traits,Param,Boundary>;
/// Empty constraints container
using NoTrafo = Dune::PDELab::EmptyTransformation;
/// Solution vector type
......@@ -90,8 +91,8 @@ public:
* \param _inifile Parameter file parser
* \param grid Grid to adapt (reference is not saved)
*/
AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>(),
WaterFluxAdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>(),
inifile(_inifile),
maxLevel(inifile.get<int>("adaptivity.maxLevel")),
minLevel(inifile.get<int>("adaptivity.minLevel")),
......@@ -220,7 +221,7 @@ public:
* define which create() function is enabled at compile time.
*/
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U, int order>
struct AdaptivityHandlerFactory
struct RichardsAdaptivityHandlerFactory
{
private:
......@@ -246,7 +247,7 @@ private:
using Grid = typename Traits::Grid;
using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>;
using AHB = RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>;
Dune::ParameterTree& inifile;
Grid& grid;
......@@ -257,17 +258,17 @@ public:
/*---------------------------------------------------------------------*//**
* @brief Create the factory, taking the parameters for building an
* AdaptivityHandler
* RichardsAdaptivityHandler
*
* @param _inifile Parameter file parser
* @param _grid Grid to adapt (reference is not saved)
*/
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
RichardsAdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
inifile(_inifile),
grid(_grid)
{ }
/// Create a dummy AdaptivityHandler
/// Create a dummy RichardsAdaptivityHandler
/** Assert that adaptive grid (UG) is used.
*/
template<bool active = enabled>
......@@ -276,15 +277,15 @@ public:
return std::make_unique<AHB>();
}
/// Create a true AdaptivityHandler
/// Create a true RichardsAdaptivityHandler
template<bool active = enabled>
std::enable_if_t<active,std::unique_ptr<AHB>> create ()
{
std::unique_ptr<AHB> p;
if (inifile.hasSub("adaptivity") &&