Commit 2a4bd16d authored by Santiago Ospina's avatar Santiago Ospina
Browse files

Merge branch 'master' into 62-gradientfluxadapter-problem



 * Created a new folder for finite elements and finite element maps

# Conflicts:
#	doc/doxygen/Doxylocal
#	dune/dorie/model/richards/richards.hh
Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parents b83f4444 113155f9
......@@ -4,20 +4,9 @@ build-cmake/
# Exclude generated files
python/dorie/wrapper/pf_from_file.py
python/dorie/wrapper/test_dorie.py
python/dorie/cli/cmds.py
python/dorie/dorie/cli/cmds.py
test/maps/cell_ids.h5
# Ignore temporary and auto-generated files #
*~
*.pyc
*.log
*.DS_Store
*/__pycache__/
*./DS_Store
*/.DS_Store?
*/._*
*/.Spotlight-V100
*/.Trashes
*/ehthumbs.db
*/Thumbs.db
# Ignore auxiliary Python files
__pycache__/
*.egg-info/
......@@ -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>
add_subdirectory("flux_reconstruction")
add_subdirectory("finite_element")
#install headers
install(FILES
......
......@@ -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 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 common
@}
@defgroup FluxReconstruction Flux reconstruction
@{
@ingroup Common
@todo document flux reconstruction
@}
**/
add_subdirectory("dg")
add_subdirectory("raviart_thomas_dg")
add_subdirectory("raviart_thomas_volume")
#install headers
install(FILES
raviart_thomas_dg_fem.hh
raviart_thomas_fem.hh
raviart_thomas_volume_fem.hh
skeleton_fem.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/dorie/common/finite_element)
\ No newline at end of file
#install headers
install(FILES
local_coefficients.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/dorie/common/finite_element/dg)
\ No newline at end of file
#ifndef DUNE_DORIE_RAVIART_THOMAS_FINITE_ELEMENT_HH
#define DUNE_DORIE_RAVIART_THOMAS_FINITE_ELEMENT_HH
#ifndef DUNE_DORIE_DG_LOCAL_COEFFICIENTS_HH
#define DUNE_DORIE_DG_LOCAL_COEFFICIENTS_HH
#include <vector>
#include <dune/localfunctions/utility/localfiniteelement.hh>
#include <dune/localfunctions/raviartthomas/raviartthomassimplex.hh>
#include <dune/localfunctions/common/localkey.hh>
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief General local keys for dg finite elements.
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @ingroup FiniteElement
*/
class DGLocalCoefficients
{
public:
/*-----------------------------------------------------------------------*//**
* @brief Constructs the coeffcients.
* @brief Constructs the coeffcients.
* @details All locall keys are associated to codimension 0.
*
* @param[in] s The size of the local basis.
*/
DGLocalCoefficients(unsigned int s)
: _s(s), li(_s)
......@@ -53,44 +61,8 @@ private:
std::vector<Dune::LocalKey> li;
};
template<unsigned int dim, class D, class R,
class SF=R, class CF=SF>
class RaviartThomasSimplexLocalFiniteElement
: public Dune::GenericLocalFiniteElement<Dune::RaviartThomasBasisFactory<dim, SF, CF>,
Dune::RaviartThomasCoefficientsFactory<dim>,
Dune::RaviartThomasL2InterpolationFactory<dim, SF> >
{
using LocalCoefficients = DGLocalCoefficients;
typedef Dune::GenericLocalFiniteElement<
Dune::RaviartThomasBasisFactory<dim, SF, CF>,
Dune::RaviartThomasCoefficientsFactory<dim>,
Dune::RaviartThomasL2InterpolationFactory<dim, SF> > Base;
public:
using Traits = Dune::LocalFiniteElementTraits<
typename Base::Traits::LocalBasisType,
LocalCoefficients,
typename Base::Traits::LocalInterpolationType>;
/** \todo Please doc me */
RaviartThomasSimplexLocalFiniteElement(const Dune::GeometryType &gt, unsigned int k)
: Base(gt, k), local_coefficients((dim == 2) ? (k+1)*(k+3) : (k+1)*(k+2)*(k+4)/2)
{
// TODO assert dim 2 and 3 and gt ==simplex
}
const typename Traits::LocalCoefficientsType& localCoefficients () const
{
return local_coefficients;
}
private:
LocalCoefficients local_coefficients;
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_RAVIART_THOMAS_FINITE_ELEMENT_HH
#endif // DUNE_DORIE_DG_LOCAL_COEFFICIENTS_HH
// 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