Commit bca7b615 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch '123-add-input-of-boundary-segmentation' into 'master'

Resolve "Add input of boundary segmentation"

Closes #123

See merge request !119
parents c3fd7c01 b3272895
......@@ -25,6 +25,7 @@
grid functions of water content and velocity.
* Infrastructure for the input of Miller scaling fields. !110
* Add logging framework 'spdlog' as submodule !106
* Support input of boundary segmentations !119
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......
......@@ -45,6 +45,7 @@ dune_cmake_sphinx_doc(SPHINX_CONF ${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in
man-installation.rst
man-config-file.rst
man-parameter-file.rst
man-grid.rst
MODULE_ONLY)
if(TARGET sphinx_html)
......
......@@ -76,28 +76,80 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> int </values>
<suggestion> 0 </suggestion>
</parameter>
</category>
<parameter name="mappingFile">
<definition> The H5 file containing the mapping from cell to medium index.
Specify the dataset inside the file with the next key. Leave empty or
set to ``None`` for a global homogeneous medium. In this case,
``grid.globalIndex`` has to be specified.
<category name="grid.mapping">
<parameter name="file">
<definition> The H5 file containing all mapping datasets.
Leave empty or set to ``None`` for global mappings.
</definition>
<values> path </values>
<suggestion> None </suggestion>
<comment> This category is only used for 'rectangular' grids </comment>
</parameter>
<parameter name="mappingFileDataset">
<definition> The H5 dataset inside ``grid.mappingFile`` containing the
mapping from cell to medium index.
<parameter name="volume">
<definition> The H5 dataset containing the
mapping from cell to medium index. May specify a global index for the
entire volume if its value can be parsed as ``int``.
</definition>
<values> path </values>
<suggestion> 0 </suggestion>
<values> path or int </values>
</parameter>
<parameter name="globalIndex">
<definition> The medium index to use for all grid cells, if
``mappingFile`` is unset or ``None``.
<parameter name="boundaryLower">
<definition> The H5 dataset mapping the lower boundary faces to
boundary condition indices. May specify a global index for the
boundary if its value can be parsed as ``int``.
</definition>
<values> int </values>
<suggestion> 0 </suggestion>
<values> path or int </values>
</parameter>
<parameter name="boundaryUpper">
<definition> The H5 dataset mapping the upper boundary faces to
boundary condition indices. May specify a global index for the
boundary if its value can be parsed as ``int``.
</definition>
<suggestion> 1 </suggestion>
<values> path or int </values>
</parameter>
<parameter name="boundaryLeft">
<definition> The H5 dataset mapping the left boundary faces to
boundary condition indices. May specify a global index for the
boundary if its value can be parsed as ``int``.
</definition>
<suggestion> 2 </suggestion>
<values> path or int </values>
</parameter>
<parameter name="boundaryRight">
<definition> The H5 dataset mapping the right boundary faces to
boundary condition indices. May specify a global index for the
boundary if its value can be parsed as ``int``.
</definition>
<suggestion> 3 </suggestion>
<values> path or int </values>
</parameter>
<parameter name="boundaryFront">
<definition> The H5 dataset mapping the front boundary faces to
boundary condition indices (3D only). May specify a global index for
the boundary if its value can be parsed as ``int``.
</definition>
<suggestion> 4 </suggestion>
<values> path or int </values>
<comment> Only in 3D </comment>
</parameter>
<parameter name="boundaryBack">
<definition> The H5 dataset mapping the back boundary faces to
boundary condition indices (3D only). May specify a global index for
the boundary if its value can be parsed as ``int``.
</definition>
<suggestion> 5 </suggestion>
<values> path or int </values>
</parameter>
</category>
......
......@@ -48,6 +48,7 @@ Manual
man-installation
python-dorie-wrapper
man-config-file
man-grid
man-bcfile
man-cookbook
man-parameter-file
......
Grid Creation and Mapping
=========================
To guarantee numeric accuracy, :doc:`boundary conditions <man-bcfile>`
and :doc:`parameterizations <man-parameter-file>` must exactly map to grid
boundary faces and grid cells, respectively. DORiE therefore not only requires
a specification of the actual grid, but also of these mappings.
DORiE supports two types of grid input, which are controlled via the config
file parameter ``grid.gridType``:
1. Building a rectangular grid on the fly (``rectangular``).
Depending on whether grid adaptivity is enabled or not, DORiE will use
a structured or an unstructured grid manager. This does not change the input
scheme. The user has to specify the config parameters for ``[grid]``
``dimensions``, spatial ``extensions``, and grid ``cells`` in each
direction. For the mapping of boundary conditions and parameters,
:ref:`mapping datasets <man-grid_mapping-dataset>` are required.
2. Building a grid according to a GMSH grid file (``gmsh``).
This mode always uses an unstructured grid manager. The user only has to
specify the ``[grid]`` ``dimensions`` and the GMSH ``gridFile``.
The mapping is read from this file and must be considered when
:ref:`building the mesh <man-grid_gmsh>`.
.. _man-grid_mapping-dataset:
Mapping Datasets
----------------
Mappings are realized as HDF5_ datasets, where every value is mapped to a
single grid entity. Changing the grid configuration requires adapting these
datasets! All config keys of this part refer to the config section
``[grid.mapping]``. There is one dataset for the mapping of the soil
architecture (``volume``) and one dataset for each segment of the boundary
(``boundaryXYZ``). The datasets contain *integer* indices, which are then
identified with certain boundary conditions or parameterizations given by
the respective YAML input files. All mapping datasets must be stored inside
the same H5 file (``file``). Use the Python module ``h5py`` for easily
creating these data files.
The dataset layout follows the C memory layout, with the primary dimension
running fastest. The coordinate system is (right-hand) cartesian, with the
x-axis running to the right in both 2D and 3D.
All datasets follow the general layout ``array[z, y, x]``, where unused
dimensions are simply omitted. The first dimension is always the vertical axis.
The ``volume`` dataset must have the same dimension as the grid itself
(``grid.dimensions``), whereas the ``boundaryXYZ`` datasets are interpreted as
projections of the boundary in their respective normal direction, reducing the
dimension by one. The ``boundaryLower`` dataset would have the layout
``array[x]`` for a 2D grid, and ``array[y, x]`` for a 3D grid. Therefore,
``boundaryLower`` and ``boundaryUpper`` are "seen from above"
(3D-layout: ``array[y, x]``), and ``boundaryLeft`` and ``boundaryRight`` are
"seen from the right" (3D-layout: ``array[z, y]``).
The following Python code creates a volume mapping dataset ``layered`` in a
file ``mapping.h5``, where the lower half of the domain maps to medium 0,
and the upper half to medium 1. The 3D domain contains 10 cells in each
direction.
.. code-block:: python
import numpy as np
import h5py
# Create dataset. Mind the data type!
size = 10
layered = np.zeros((size, size, size), dtype=np.int_)
layered[5:, ...] = 1
# Write dataset to file in 'append' mode
with h5py.File("mapping.h5", 'a') as file:
file.create_dataset("layered", data=layered)
Homogeneous Mappings
^^^^^^^^^^^^^^^^^^^^
If the entire volume, or a complete boundary, should be mapped to a single
index, no dataset is required. Instead, you can set the value for the
respective config file parameter to the desired index. If no dataset should be
read at all, set ``file`` to ``none``.
Even with ``file`` set to a valid H5 file, DORiE will *always* try to
interpret the input values for ``grid.mapping`` as integers. If this succeeds,
the value is interpreted as "global" mapping index for the respective part
of the grid. Therefore, **do not use dataset names starting with digits!**
.. _man-grid_gmsh:
Mapping Soil Layers with GMSH
-----------------------------
GMSH_ supports mapping lines, surfaces, and volumes to *physical* entities.
These entities may combine multiple of the aforementioned *geometric*
entities. Physical entities are assigned a non-negative index upon creation.
These indices can be set by the user in any of the GMSH frontends.
Indices may not occur multiple times, even if they are assigned to different
types of physical entities.
The following ``.geo`` GMSH input file creates a rectangular domain that is
split in half. The lower part is mapped to index 1, the upper part to index 2.
Additionally, a different index for every boundary orientation is set. Notice
that the left and right boundaries consist of two *geometric* lines each.
This code can directly be transferred to a Python script using
the pygmsh_ module. It writes a ``.geo`` file while checking for a correct
syntax within your script. It is readily installed into the
:doc:`virtual environment <python-dorie-wrapper>`.
.. code-block:: default
// define geometric entities
Point(1) = {0, 0, 0, 0.1};
Point(2) = {2, 0, 0, 0.1};
Point(3) = {2, 2, 0, 0.1};
Point(4) = {0, 2, 0, 0.1};
Point(5) = {0, 1, 0, 0.1};
Point(6) = {2, 1, 0, 0.1};
Line(1) = {1, 2};
Line(2) = {2, 6};
Line(3) = {6, 3};
Line(4) = {3, 4};
Line(5) = {4, 5};
Line(6) = {5, 1};
Line(7) = {5, 6};
Line Loop(1) = {1, 2, -7, 6};
Plane Surface(1) = {1};
Line Loop(2) = {7, 3, 4, 5};
Plane Surface(2) = {2};
// define physical entities, index in round brackets
Physical Surface(1) = {1}; // lower
Physical Surface(2) = {2}; // upper
// entire set of physical entities must always be defined!
Physical Line(3) = {1}; // bottom
Physical Line(4) = {2, 3}; // right
Physical Line(5) = {4}; // top
Physical Line(6) = {5, 6}; // left
A ``.geo`` file is the basis for creating the actual mesh in GMSH. You can
load it into the GMSH GUI, or perform the meshing directly using the
`GMSH command line interface
<http://gmsh.info/doc/texinfo/gmsh.html#Non_002dinteractive-mode>`_:
gmsh <geo-file> -<dim>
Replace ``<geo-file>`` with the appropriate file, and ``dim`` with the
spatial dimension of the intended mesh.
.. _HDF5: https://www.h5py.org/
.. _GMSH: http://gmsh.info/
.. _pygmsh: https://pypi.org/project/pygmsh/
Soil Parameters and Architecture
================================
Soil Parameterization
=====================
Specifying the domain properties requires input of the soil architecture
Specifying the domain properties requires input of the
:doc:`soil architecture <man-grid>`
(in relation to the grid) and of the soil parameters. The latter incorporate
the *subscale physics*, the representation of the soil hydraulic properties
below the REV scale. DORiE requires several input files for retrieving this
......@@ -9,10 +10,8 @@ information, depending on the type of grid used for the computation.
A YAML_ parameter file is always required. It specifies a set of soil
parameterizations, along with a medium index. This index is used to reference
the medium specified in the architecture files. There are two types of these:
If you insert an unstructured grid with GMSH, you have to incorporate the
medium information. If you create a rectangular grid, you need to add a mapping
dataset, specifying the medium index for each (initial) grid cell.
the medium specified in the architecture files, or with the "global" indices
given in the config file.
Additionally, you can specify settings for the small-scale heterogeneities
of the soil. They require an input via H5 datasets that contain appropriate
......@@ -40,8 +39,8 @@ of scaling.
Every ``scaling_section`` has the same keys: The H5 filepath is given by
``file``, and the internal dataset path by ``dataset``. You can then choose an
``interpolation`` method for the dataset, which requires you to insert values
for the physical ``extensions`` of the dataset and the ``offset`` of its (front)
lower left corner from the grid origin.
for the physical ``extensions`` of the dataset and the ``offset`` of its
(front) lower left corner from the grid origin.
.. code-block:: yaml
......@@ -73,98 +72,6 @@ The parameterization ``type`` must match the static ``type`` member of one of
them. Likewise, the parameters must match the names of parameters
these objects use.
Medium Mapping Dataset
----------------------
This dataset is required when creating a rectangular grid, whether adaptive or
not. It has to be realized as dataset inside an HDF5_ file. The dataset must
have the same dimensions as the target grid, and must specify the medium index
for every cell of this grid. As all medium indices are integers, the dataset
must contain **integers, not floats**! An HDF5 file with such a dataset can
easily be created with the Python module ``h5py``.
The dataset layout follows the C memory layout, with the primary dimension
running fastest. A 3D HDF5 dataset or a 3D ``numpy`` array have the layout
``array[z][y][x]``, and 2D arrays have the layout ``array[y][x]``. The first
dimension is always the vertical axis.
The following Python code creates a dataset ``layered`` in a file
``mapping.h5``, where the lower half of the domain maps to medium 0, and the
upper half to medium 1. The 3D domain contains 10 cells in each direction.
.. code-block:: python
import numpy as np
import h5py
# Create dataset. Mind the data type!
size = 10
layered = np.zeros((size, size, size), dtype=np.int_)
layered[5:, ...] = 1
# Write dataset to file
with h5py.File("mapping.h5", 'a') as file:
f.create_dataset("layered", data=layered)
The dataset and the file containing it are specified in the
:doc:`config file <man-config-file>` via the keys ``grid.mappingFile`` and
``grid.mappingFileDataset``.
Homogeneous Media
^^^^^^^^^^^^^^^^^
In the special case of a globally homogeneous medium, you can omit the mapping
dataset. This is indicated by leaving ``grid.mappingFile`` undefined or setting
its value to ``None``. You can then define the medium index for the entire
domain with the key ``grid.globalIndex``, which is unused otherwise.
Mapping Soil Layers with GMSH
-----------------------------
GMSH_ supports mapping lines, surfaces, and volumes to *physical* entities.
These entities may combine multiple of the aforementioned *geometrical*
entities. Physical entities are assigned a non-negative index upon creation.
These indices can be set by the user in any of the GMSH frontends.
Indices may not occur multiple times, even if they are assigned to different
types of physical entities.
The following ``.geo`` GMSH input file creates a rectangular domain that is
split in half. The lower part is mapped to index 1, the upper part to index 2.
This code can directly tranferred to a Python script using the pygmsh_ module.
It writes a ``.geo`` file while checking for a correct syntax within your
script. It is readily installed into the
:doc:`virtual environment <python-dorie-wrapper>`.
.. code-block:: default
// define geometric entities
Point(1) = {0, 0, 0, 0.1};
Point(2) = {2, 0, 0, 0.1};
Point(3) = {2, 2, 0, 0.1};
Point(4) = {0, 2, 0, 0.1};
Point(5) = {0, 1, 0, 0.1};
Point(6) = {2, 1, 0, 0.1};
Line(1) = {1, 2};
Line(2) = {2, 6};
Line(3) = {6, 3};
Line(4) = {3, 4};
Line(5) = {4, 5};
Line(6) = {5, 1};
Line(7) = {5, 6};
Line Loop(1) = {1, 2, -7, 6};
Plane Surface(1) = {1};
Line Loop(2) = {7, 3, 4, 5};
Plane Surface(2) = {2};
// define physical entities, index in round brackets
Physical Surface(1) = {1}; // lower
Physical Surface(2) = {2}; // upper
// entire set of physical entities must always be defined!
Physical Line(3) = {1}; // bottom
Physical Line(4) = {2, 3}; // right
Physical Line(5) = {4}; // top
Physical Line(6) = {5, 6}; // left
Scaling Field Datasets
----------------------
One or more datasets in possibly multiple HDF5_ files are required for creating
......@@ -295,5 +202,3 @@ Their ``type`` and the names of their member parameters must match.
.. _YAML: https://gettaurus.org/docs/YAMLTutorial/
.. _HDF5: https://www.h5py.org/
.. _GMSH: http://gmsh.info/
.. _pygmsh: https://pypi.org/project/pygmsh/
......@@ -4,16 +4,9 @@
#include <memory>
#include <array>
#include <vector>
#include <iostream>
#include <iomanip>
#include <type_traits>
#include <map>
#include <string>
#include <mpi.h>
#include <dune/common/timer.hh>
#include <dune/common/ios_state.hh>
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
......@@ -26,22 +19,19 @@
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/dorie/common/grid_mapper.hh>
#include <dune/dorie/common/h5file.hh>
namespace Dune {
namespace Dorie {
/// Type of grid loading mode. Defines how parameters are read from the grid.
enum class GridMode {
gmsh, //!< Read a GMSH file containing parameter and boundary mappings
rectangular //!< Build a grid internally and read parameter mappings from files
};
/// This class builds grids by reading GMSH files or using the Dune GridFactory.
/** Additionally, it provides the mapping for associating grid elements with
* soil layers (parameter sets) via the GridMapper.
*
* \tparam Grid The type of the grid to build
* \ingroup Grid
* \author Lukas Riedel
* \date 2018
*/
template<typename Grid>
class GridCreator
......@@ -55,50 +45,64 @@ private:
const Dune::ParameterTree& _config;
//! Verbosity of this class
const std::shared_ptr<spdlog::logger> _log;
//! The mode for reading the grid.
//! The mode for reading the grid.
GridMode _grid_mode;
//! Typedef for the created grid mapper
using GridMapper = typename Dune::Dorie::GridMapper<Grid>;
using GridMapper = typename Dune::Dorie::GridMapperBase<Grid>;
/// Type of the maps returned by the GridMapper
using Map = typename GridMapper::Map;
/// The grid
std::shared_ptr<Grid> _grid;
/// The element mapping
Map _element_index_map;
/// Boundary mapping
Map _boundary_index_map;
public:
/// Build the grid and retrieve mapping information
/**
* \param config User config settings
* \param helper Instance of the MPI helper
*/
explicit GridCreator (const Dune::ParameterTree& config,
const Dune::MPIHelper& helper):
_helper(helper),
_config(config),
_log(Dorie::get_logger(log_base))
{
check_parallel_allowed();
const auto grid_type = _config.get<std::string>("grid.gridType");
if (grid_type == "rectangular")
{
_grid_mode = GridMode::rectangular;
else if (grid_type == "gmsh")
_grid_mode = GridMode::gmsh;
else {
_log->error("Unsupported grid type: {}", grid_type);
DUNE_THROW(IOError, "Unsupported grid type!");
}
check_parallel_allowed();
}
_grid = create_rect_grid();
/// Return the creared grid and the parameter mapping
/** Building the grid and subsequently the GridMapper ensures that the grid
* and the mapping data are properly load-balanced. The GridMapper manages
* load-balancing and data broadcasting in its constructor.
* \return Shared pointer to Dorie::GridMapper.
*/
std::shared_ptr<GridMapper> get_mapper ()
{
if (_grid_mode == GridMode::gmsh)
Dorie::GridMapper<Grid, GridMode::rectangular> mapper(_grid, _config);
_element_index_map = mapper.element_index_map();
_boundary_index_map = mapper.boundary_index_map();
}
else if (grid_type == "gmsh")
{
if constexpr (std::is_same_v<Grid, Dune::UGGrid<2>> or
std::is_same_v<Grid, Dune::UGGrid<3>>)
{
_grid_mode = GridMode::gmsh;
auto [grid,
element_index_map,
boundary_index_map] = create_gmsh_grid();
return
std::make_shared<GridMapper>(
grid, element_index_map, boundary_index_map);
element_index_map,
boundary_index_map] = create_gmsh_grid();
_grid = grid;
Dorie::GridMapper<Grid, GridMode::gmsh> mapper(_grid,
element_index_map,
boundary_index_map);
_element_index_map = mapper.element_index_map();
_boundary_index_map = mapper.boundary_index_map();
}
else {
_log->error("DORiE expects to work on UGGrid for "
......@@ -107,18 +111,19 @@ public:
<< "algorithms");
}
}
else
{
auto [grid,
element_index_map,
boundary_index_map,
cells] = create_rect_grid();
return
std::make_shared<GridMapper>(
grid, element_index_map, boundary_index_map, cells);
else {
_log->error("Unsupported grid type: {}", grid_type);
DUNE_THROW(IOError, "Unsupported grid type!");
}
}
/// Return a shared pointer to the created grid
std::shared_ptr<Grid> grid () const { return _grid; }
/// Return the local element index map
const Map& element_index_map () const { return _element_index_map; }
/// Return the local boundary index map
const Map& boundary_index_map () const { return _boundary_index_map; }
/// Return true if parallel execution is allowed. Issue a warning if not.
/** Currently, we only allow parallel execution on YaspGrid due to
* constraints imposed by the linear solver.
......@@ -176,10 +181,9 @@ private:
* \return Tuple: Shared pointer to grid; Element index map;
* Boundary index map; Number of cells in each direction
*/
std::tuple<std::shared_ptr<Grid>, std::vector<int>, std::vector<int>, std::vector<int>>
std::shared_ptr<Grid>
create_rect_grid () const
{
Dune::Timer timer;
const auto level = _config.get<int>("grid.initialLevel");
const auto upperRight = _config.get<Dune::FieldVector<double, _dim>>("grid.extensions");
const auto elements = _config.get<std::array<unsigned int, _dim>>("grid.cells");
......@@ -193,11 +197,7 @@ private:
_log->debug(" Applying global refinement of level: {}", level);
grid->globalRefine(level);
// get the element index map from an H5 file
auto [element_index_map, cells] = read_element_index_map_from_file();
std::vector<int> boundary_index_map; // dummy variable
return std::make_tuple(grid, element_index_map, boundary_index_map, cells);
return grid;
}
/// Default rectangular grid constructor. Calls StructuredGridFactory.
......@@ -240,85 +240,6 @@ private:
return std::make_shared<YaspGrid<_dim>>(extensions, elements_);
}
/// Read the element mapping from an H5 file
/** This reads the "grid.mappingFile" if given and returns the resulting
* map. If not given, "grid.globalIndex" is read and broadcast to the map.
* \return Tuple: Element index map; number of cells in each direction.
*/
std::tuple<std::vector<int>, std::vector<int>>
read_element_index_map_from_file () const
{
// get number of cells
const auto elements = _config.get<
std::array<int, _dim>>("grid.cells");
// check if mapping file is given
const auto file_path = _config.get<std::string>("grid.mappingFile", "");
// instantiate return variables
std::vector<int> element_index_map;
std::vector<int> cells(_dim);
// check for global homogeneous grid index
if (file_path == "None"
or file_path == "none"
or file_path.empty())
{