Commit 89376315 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'master' into 121-add-linear-interpolator-and-use-it-for-scaling-fields

parents 41e93e7e 358286cc
......@@ -41,6 +41,8 @@
* Initial conditions generated from H5 input data !130
* Generic Python VTK file reader !143
* Linear interpolator for initial conditions and scaling fields !145
* Parameterizations for hydrodynamic dispersion in solute transport !141
* Generic Python VTK file reader !143, !150
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......
......@@ -197,7 +197,8 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
to the `CMAKE_FLAGS` in the above command, replacing `<hdf5-path>` with the
path chosen as installation prefix when configuring HDF5.
#### Experimental Features
### Experimental Features
The local operator implementing Richards equation's discretization supports
multiple scheme settings. Setting these via the config file is disabled by
default. You can enable this feature by reconfiguring DORiE with the CMake flag
......
......@@ -7,7 +7,8 @@ add_custom_target(unit_tests
# Add the system test build target and test command
add_custom_target(build_system_tests)
add_custom_target(system_tests
COMMAND ctest --output-on-failure --exclude-regex ^ut.+$
COMMAND ctest
--output-on-failure --exclude-regex "^\\(ut.+\\|python_tests\\)$"
)
# create the mapping datafile as preparation for all tests
......
......@@ -119,7 +119,7 @@ endfunction()
# copy BC & parameter files
file(COPY 2d_infiltr.bcdat 3d_infiltr.bcdat
2d_solute.bcdat 3d_solute.bcdat
param.yml
richards_param.yml transport_param.yml
DESTINATION .)
# Random field generator
......
......@@ -32,10 +32,15 @@ adding an empty line, make text **bold** or ``monospaced``.
<dorie>
<category name="simulation">
<parameter name="mode">
<definition> Mode of simulation in DORiE. </definition>
<definition> Sets the simulation mode of DORiE.
(``richards``) mode calculates the matric head with a
DG scheme and produce water fluxes with Raviart Thomas reconstruction.
(``richards+transport``) mode calculates (``richards``) mode and use
the generated water flux and saturation at each step to calculate the
solute transport model for unsaturated media with a DG scheme.
</definition>
<suggestion> richards </suggestion>
<values> richards, richards+transport </values>
<comment> richards, richards+transport </comment>
</parameter>
</category>
......
volumes:
sand:
index: 0
type: MvG
parameters:
alpha: -2.3
n: 4.17
k0: 2.2E-5
theta_r: 0.03
theta_s: 0.31
tau: -1.1
silt:
index: 1
type: MvG
parameters:
alpha: -0.7
n: 1.3
k0: 1.0E-5
theta_r: 0.01
theta_s: 0.41
tau: 0.0
scaling:
type: None
......@@ -227,6 +227,7 @@ adding an empty line, make text **bold** or ``monospaced``.
<definition> YAML file containing the parameter definitions.
</definition>
<values> path </values>
<suggestion> richards_param.yml </suggestion>
</parameter>
</category>
......
volumes:
sand:
index: 0
richards:
type: MvG
parameters:
alpha: -2.3
n: 4.17
k0: 2.2E-5
theta_r: 0.03
theta_s: 0.31
tau: -1.1
silt:
index: 1
richards:
type: MvG
parameters:
alpha: -0.7
n: 1.3
k0: 1.0E-5
theta_r: 0.01
theta_s: 0.41
tau: 0.0
scaling:
type: None
......@@ -67,6 +67,16 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> endOfTransportStep, endOfRichardsStep, none </values>
</parameter>
<parameter name="writeDispersionTensor">
<definition> Defines whether VTK files should include the hydrodynamic
dispersion tensor. Tensors are written in 3D and have 9 componentents
independently of the world dimension. This can be easily be visualizated
in Paraview with the ``Tensor Glyph`` filter.
</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
......@@ -109,6 +119,13 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="dirichletMode">
<definition> Type of the input value for the dirichlet condition.
</definition>
<values> soluteConcentration, totalSolute </values>
<suggestion> soluteConcentration </suggestion>
</parameter>
</category>
<category name="initial">
......@@ -222,11 +239,12 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="parameters">
<parameter name="effHydrDips">
<definition> Effective hydrodynamic dispersion </definition>
<values> float &gt; 0 </values>
<suggestion> 2E-9 </suggestion>
<category name="parameters">
<parameter name="file">
<definition> YAML file containing the parameter definitions.
</definition>
<values> path </values>
<suggestion> transport_param.yml </suggestion>
</parameter>
</category>
......@@ -253,8 +271,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="FEorder">
<definition> Order of the finite element method used. Values &gt; 1 are not
thoroughly tested. </definition>
<definition> Order of the finite element method used. </definition>
<suggestion> 0 </suggestion>
<values> 0 </values>
</parameter>
......
volumes:
sand:
index: 0
transport:
type: Deff_MQ1 + Dhm_iso
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
silt:
index: 1
transport:
type: Deff_MQ1 + Dhm_iso
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
\ No newline at end of file
......@@ -64,11 +64,16 @@ Output Files
.. object:: Transport
- ``solute``: Solute concentration in water phase
:math:`c_w \, [\mathrm{kg}/\mathrm{m}^3]`
:math:`c_w \, [\mathrm{kg}\mathrm{m}^{-3}]`
- ``solute_total``: Total solute concentration
:math:`c_t \, [\mathrm{kg}/\mathrm{m}^3]`
:math:`c_t \, [\mathrm{kg}\mathrm{m}^{-3}]`
- ``micro_peclet``: Microscopic peclet number
:math:`\mathsf{pe_\mu} \, [-]`
- ``flux_RT<k-1>``: Reconstructed solute flux
:math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}/\mathrm{m}^{-2}\mathrm{s}^{-1}]`
:math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}\,\mathrm{m}^{-2}\mathrm{s}^{-1}]`
(if enabled!)
- ``eff_hd_dispersion``: Hydrodynamic dispersion tensor
:math:`\mathsf{D}_\mathsf{hd} \, [\mathrm{m}^2\mathrm{s}^{-1}]`
(if enabled!)
* VTK Parallel Collection Files (``.pvtu``): Merging multiple VTK files in case
......
......@@ -7,7 +7,7 @@ Cook Book
Introduction
============
This part of the documentation is intended for first-time DORiE users. It explaines the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement.
This part of the documentation is intended for first-time DORiE users. It explains the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement.
Prerequisites
-------------
......@@ -57,7 +57,7 @@ Let's begin with the ``output``. We want the program to give us at least some ou
.. code-block:: none
[output]
verbose = 1
logLevel = info
outputPath = ./sand
fileName = sand
......
This diff is collapsed.
......@@ -53,12 +53,17 @@ Generic VTK Reader
------------------
This is a flexible VTK XML file reader which can evaluate all data sets written
by DORiE.
by DORiE. It is derived from the abstract base class ``Mapping`` and thus
supports accessing the data arrays stored inside like a Python ``dict``. Use
the data array names as keys.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKReader
:members:
.. automethod:: __init__
The ``VTKReader`` returns instances of ``VTKDataArray``, which is a wrapper
around the underlying VTK data array.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKDataArray
:members:
dorie.utilities.check_path
--------------------------
......
Sphinx<2.0
Sphinx
sphinx_rtd_theme
sphinx-argparse
sphinx-markdown-tables
......
......@@ -816,6 +816,8 @@ namespace Dorie{
std::cout << "outflow : ";
else if(bcData[i][j].type[k]==BoundaryCondition::Dirichlet)
std::cout << "dirichlet : ";
else if(bcData[i][j].type[k]==BoundaryCondition::Outflow)
std::cout << "outflow : ";
}
index++;
}
......
......@@ -16,6 +16,72 @@
namespace Dune{
namespace Dorie{
/**
* @brief Converts PDELab grid function to match a dune-function.
* @note This adapter differs from the one provided in PDELab in that this
* one has the current interface of dune-grid to write vtk files.
* The main difference using this adapter is that it allows to write
* more than 3 components.
*
* @tparam GF PDELab grid function type.
*/
template<class GF>
class VTKGridFunctionAdapter
{
using GridFunction = typename std::decay<GF>::type;
using GridView = typename GF::Traits::GridViewType;
using Entity = typename GridView::template Codim<0>::Entity;
using Range = typename GF::Traits::RangeType;
using Domain = typename GF::Traits::DomainType;
public:
/**
* @brief Constructs the object
*
* @param[in] gf pointer to a constant grid function
*/
VTKGridFunctionAdapter(std::shared_ptr<const GF> gf)
: _gf(gf)
, _entity(nullptr)
{}
/**
* @brief Binds an entity to the grid function
*
* @param[in] e Entity
*/
void bind(const Entity& e)
{
_entity = &e;
}
/**
* @brief Unbinds the previously binded entity
*/
void unbind()
{
_entity = nullptr;
}
/**
* @brief Evaluates the grid function at a given coordinate
*
* @param[in] x Local coordinate with respect to the last binded entity
*/
Range operator()(const Domain& x) const
{
Range y;
assert(_entity);
_gf->evaluate(*_entity,x,y);
return y;
}
private:
const Entity* _entity;
std::shared_ptr<const GF> _gf;
};
/*-------------------------------------------------------------------------*//**
* @brief Class for grid fucntion vtk sequence writer. This class works
* exacly the same way as a VTKSequence writer but it receives
......@@ -39,13 +105,20 @@ public:
* @tparam GF Type of the grid function
*/
template<class GF>
void addVertexData(const std::shared_ptr<const GF> &gf_ptr, std::string name)
void addVertexData(std::shared_ptr<const GF> gf_ptr, std::string name)
{
static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value,
"GridFunctionVTKSequenceWriter only works with only one GridView type.");
auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
Dune::VTKSequenceWriter<GV>::addVertexData(vtk_gf_ptr);
Dune::Dorie::VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
const int vtk_ncomps = GF::Traits::dimRange;
const int dim = GF::Traits::dimDomain;
auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector
: VTK::FieldInfo::Type::scalar;
auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps);
this->vtkWriter()->addVertexData(vtk_gf,vtk_info);
}
/*-----------------------------------------------------------------------*//**
......@@ -62,14 +135,21 @@ public:
static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value,
"GridFunctionVTKSequenceWriter only works with only one GridView type.");
auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
Dune::VTKSequenceWriter<GV>::addCellData(vtk_gf_ptr);
Dune::Dorie::VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
const int vtk_ncomps = GF::Traits::dimRange;
const int dim = GF::Traits::dimDomain;
auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector
: VTK::FieldInfo::Type::scalar;
auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps);
this->vtkWriter()->addCellData(vtk_gf,vtk_info);
}
/*-----------------------------------------------------------------------*//**
* @brief Clear all the attached VTKFunctions. (Necessary for each
* iteration if VTKFunctions are managed internaly with pointers
* instead of references)
* @brief Clear all the attached VTKFunctions.
* @details This is necessary for each iteration if VTKFunctions are
* managed internaly with pointers instead of references
*/
void clear()
{
......
......@@ -20,7 +20,7 @@ namespace Dorie {
* @brief Class for analytic initial conditions.
*
* @ingroup InitialConditions
* @author Santiago Ospina
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @tparam T BaseTraits
......
......@@ -17,7 +17,7 @@ namespace Dorie{
* @brief Abstract Base Class for initial conditions.
*
* @ingroup InitialConditions
* @author Santiago Ospina
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @tparam T BaseTraits
......
......@@ -18,7 +18,7 @@ namespace Dorie{
* @brief Class for initial conditions from hdf5 data files.
*
* @ingroup InitialConditions
* @author Santiago Ospina
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @tparam T BaseTraits
......
......@@ -13,7 +13,7 @@ namespace Dorie{
* @brief Abstract Base Class for initial conditions.
*
* @ingroup InitialConditions
* @author Santiago Ospina
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @tparam T BaseTraits
......
#ifndef DUNE_DORIE_PARAMETERIZATION_READER_HH
#define DUNE_DORIE_PARAMETERIZATION_READER_HH
#include <unordered_map>
#include <memory>
#include <string>
#include <set>
#include <yaml-cpp/yaml.h>
#include <dune/common/exceptions.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Interface and default reader for parameterization factories
*
* @tparam Param The parameter type to build
*/
template<class Param>
struct ParameterizationFactory
{
virtual ~ParameterizationFactory() = default;
/**
* @brief Parameterization selector
* @details Given a type_node in yaml, it selects a -polymorphic-
* parameter.
*
* @param[in] type_node YAML node describing the type of parameterization to
* choose from.
* @param[in] name Name that describes the parametrized material.
*
* @return Shared pointer to the selected parameter.
*/
virtual std::shared_ptr<Param> selector(
const YAML::Node& type_node,
const std::string& name) const = 0;
/**
* @brief Parameterization reader
* @details It reads parameterizations and its parameters. For each volume,
* it assings a index to a parameterization (choosen by the
* parameterization selector) and later filled with the paramaters
* in the yaml file.
*
* @param[in] param_file YAML node with information about type and
* parameters
* @param[in] param_name Name of the parameterization.
* @param[in] log Logger
*
* @return (Unsorted) Map with parameterizations for each material.
*/
std::unordered_map<int, std::shared_ptr<Param>>
reader( const YAML::Node& param_file,
const std::string& param_name,
const std::shared_ptr<spdlog::logger> log) const
{
// mapping: index on grid to parameterization object (will be returned)
std::unordered_map<int, std::shared_ptr<Param>> ret;
// set to check for duplicate indices
std::set<int> param_index_set;
// iterate over all volume nodes
for (auto&& node : param_file["volumes"])
{
const auto volume = node.first.as<std::string>();
const YAML::Node param_info = node.second;
// read relevant data
const auto index = param_info["index"].as<int>();
if (param_index_set.count(index) != 0) {
log->error("Index '{}' in volume '{}' of parameterization '{}'"
"already exists.",
index, volume, param_name);
DUNE_THROW(IOError, "Parameter file indices must be unique");
}
param_index_set.emplace(index);
// get parameterization object
const YAML::Node& type_node = param_info[param_name];
auto parameterization = this->selector(type_node, volume);
auto parameter_values = parameterization->parameters();
// fill parameterization with parameter entries
log->trace(" Reading parameters. Volume: {}, Index: {}",
volume, index);
auto yaml_parameters = yaml_sequence_to_map(type_node["parameters"]);
for (auto&& [name, value] : parameter_values)
{
try {
value = yaml_parameters[name].template as<double>();
}
catch (...) {
log->error(" Parameter '{}' is missing in '{}' parameterization,"
" volume {}, index {}.",
name, param_name, volume, index);
DUNE_THROW(IOError, "Missing parameters!");
}
}
ret.emplace(index, parameterization);
}
return ret;
}
};
}
}
#endif // DUNE_DORIE_PARAMETERIZATION_READER_HH
\ No newline at end of file
......@@ -63,7 +63,7 @@ namespace Setup
const std::string inifilename = argv[1];
// Read ini file
log->info("Reading parameter file: {}", inifilename);
log->info("Reading configuration file: {}", inifilename);
Dune::ParameterTree inifile;
Dune::ParameterTreeParser ptreeparser;
ptreeparser.readINITree(inifilename, inifile);
......
......@@ -13,7 +13,7 @@ namespace Dune{
enum Type {
Neumann, //!< Fixed flux at boundary
Dirichlet, //!< Fixed matric head at boundary
Outflow, //!< Fixed flux at boundary
Outflow, //!< Fixed outflow boundary condition
Other //!< UNUSED
};
......
......@@ -195,16 +195,6 @@ 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
......@@ -214,7 +204,6 @@ struct BaseTraits
using TF = double;
using TimeField = TF;
using TimeInterval = Dune::Dorie::TimeInterval<TimeField>;
using RF = double;
using RangeField = RF;
......
......@@ -8,6 +8,9 @@
#include <memory>
#include <utility>
#include <yaml-cpp/yaml.h>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/gridview.hh>
......@@ -113,6 +116,61 @@ inline bool is_none (const std::string& in)
return false;
}
/*-------------------------------------------------------------------------*//**
* @brief Transform a Yaml node with sequences to map
* @details This function check all the nodes within 'in' and transform every
* sequence into a map. In particular it transform vectors and
* tensors into single key-value map pair by adding a suffix
* indicating the component.
*
* @param[in] in Map Yaml node.
*
* @return Map Yaml node with only key-double pairs.
*/
inline YAML::Node yaml_sequence_to_map(const YAML::Node& in)
{
YAML::Node out;
for (auto&& node : in)
{
const auto key = node.first.as<std::string>();
const YAML::Node& value = node.second;
if (value.IsScalar()) {
out[key] = value;
} else if (value.IsSequence()) {
if (value.size() == 1) { // scalar
out[key] = value[0].as<double>();
} else if (value.size() == 2) { // assume vector in 2D
out[key+"_x"] = value[0].as<double>();
out[key+"_y"] = value[1].as<double>();
} else if (value.size() == 3) { // assume vector in 3D
out[key+"_x"] = value[0].as<double>();
out[key+"_y"] = value[1].as<double>();
out[key+"_z"] = value[2].as<double>();
} else if (value.size() == 4) { // assume tensor in 2D
out[key+"_xx"] = value[0].as<double>();
out[key+"_xy"] = value[1].as<double>();
out[key+"_yx"] = value[2].as<double>();
out[key+"_yy"] = value[3].as<double>();
} else if (value.size() == 9) { // assume tensor in 3D
out[key+"_xx"] = value[0].as<double>();
out[key+"_xy"] = value[1].as<double>();
out[key+"_xz"] = value[2].as<double>();
out[key+"_yx"] = value[3].as<double>();
out[key+"_yy"] = value[4].as<double>();
out[key+"_yz"] = value[5].as<double>();
out[key+"_zx"] = value[6].as<double>();
out[key+"_zy"] = value[7].as<double>();
out[key+"_zz"] = value[8].as<double>();
}
} else {
DUNE_THROW(Dune::IOError,
"Yaml sequence '" << key << "' is not convertible to map pair!");
}
}
return out;
}
} // namespace Dorie
} // namespace Dune
......
......@@ -9,7 +9,8 @@ namespace Dune{
/*---------------------------------------------------------------------*//**
* @brief Converts a parametrization P that contains conductivity
* information into a grid function (in PDELab sense)
*
* @ingroup RichardsParam
*
* @tparam T BasicTraits with information about the function domain
* @tparam P Parametrization class with saturated conductivity
* information
......