Skip to content

GitLab

D
dorie
Commit 0507c49d authored by Lukas Riedel's avatar Lukas Riedel
Browse files
Options

Rebuild grid creation structures and routines 

* Add 'GridCreator'  class for building the grid
* Add 'GridMapper' class for managing load-balanced mapping
* Update 'H5File' to read dataset with arbitrary data types

*This breaks everything else!* [ci skip]
parent 43cc4607
Hide whitespace changes
Inline Side-by-side
Showing with 423 additions and 156 deletions
dune/dorie/solver/util_grid_creator.hh
View file @ 0507c49d
......@@ -5,7 +5,12 @@
#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>
......@@ -19,148 +24,257 @@
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include "util_grid_mapper.hh"
#include "util_h5tools.hh"
namespace Dune {
namespace Dorie {
/// Check if parallel execution is allowed
/** Currently, we only allow parallel execution on YaspGrid due to
* constraints imposed by the linear solver
*/
template<typename GridType>
void check_parallel_allowed (const Dune::MPIHelper& helper)
{
if(helper.size() > 1
&& !std::is_same<YaspGrid<GridType::dimension>,GridType>::value) {
DUNE_THROW(Dune::Exception,"DORiE does not support parallel execution with this grid configuration!");
}
}
/// Specialized constructor for parallel YaspGrid construction
/** \param extensions Extensions of the grid
* \param elements Number of elements in each direction
* \return Shared pointer to the grid
/// 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
regular //!< 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
*/
template<typename GridType>
std::enable_if_t<
std::is_same<YaspGrid<GridType::dimension>,GridType>::value,
std::shared_ptr<YaspGrid<GridType::dimension>> >
grid_cube_construction (const Dune::FieldVector<double,GridType::dimension>& extensions,
const std::array<unsigned int,GridType::dimension>& elements)
template<typename Grid>
class GridCreator
{
std::array<int,GridType::dimension> elements_;
std::copy(elements.cbegin(),elements.cend(),elements_.begin());
private:
//! Spatial dimensions of the grid
static constexpr int _dim = Grid::dimension;
//! Dune MPI helper for querying processor ranks
const Dune::MPIHelper& _helper;
//! Config file tree
const Dune::ParameterTree& _config;
//! Verbosity of this class
const int _verbose;
//! The mode for reading the grid.
GridMode _grid_mode;
return std::make_shared<YaspGrid<GridType::dimension>>(extensions,elements_);
}
public:
explicit GridCreator (const Dune::MPIHelper& helper,
const Dune::ParameterTree& config):
_helper(helper),
_config(config),
_verbose(_config.get<int>("output.verbose"))
{
const auto grid_type = _config.get<std::string>("grid.type");
if (grid_type == "rectangular")
_grid_mode = GridMode::regular;
else if (grid_type == "gmsh")
_grid_mode = GridMode::gmsh;
else
DUNE_THROW(IOError, "Unsupported grid type!");
/// Default rectangular grid constructor. Call StructuredGridFactory
/** \param extensions Extensions of the grid
* \param elements Number of elements in each direction
* \return Shared pointer to the grid
*/
template<typename GridType>
std::enable_if_t<
!std::is_same<YaspGrid<GridType::dimension>,GridType>::value,
std::shared_ptr<GridType> >
grid_cube_construction (const Dune::FieldVector<double,GridType::dimension>& extensions,
const std::array<unsigned int,GridType::dimension>& elements)
{
const Dune::FieldVector<double,GridType::dimension> origin(.0);
// this->check_parallel_allowed();
}
return Dune::StructuredGridFactory<GridType>::createCubeGrid(origin,extensions,elements);
}
/// Return the creared grid and the parameter mapping
std::tuple<std::shared_ptr<Grid>, GridMapper<Grid>> get_grid_and_mapper ()
{
std::shared_ptr<Grid> grid;
std::vector<int> element_index_map;
std::vector<int> boundary_index_map;
if (_grid_mode == GridMode::gmsh) {
std::tie(grid,
element_index_map,
boundary_index_map) = create_gmsh_grid();
}
else {
std::tie(grid,
element_index_map,
boundary_index_map) = create_rect_grid();
}
/// Read simplex Gmsh file in 2D or 3D, create a grid, and return a shared pointer to the grid
/** Refine the grid globally and call default load balancing for multi-processing support.
* \tparam GridType Type of the grid. Must be supported by Dune::GridFactory
* \param inifile Parameter file parser
* \return Shared pointer to grid
*/
template<class GridType>
std::shared_ptr<GridType> build_grid_gmsh (const Dune::ParameterTree& inifile, const Dune::MPIHelper& helper)
{
check_parallel_allowed<GridType>(helper);
Dune::Timer timer;
ios_base_all_saver restore(std::cout);
enum { dim = GridType::dimension };
const int verbose = inifile.get<int>("output.verbose");
const int level = inifile.get<int>("grid.initialLevel");
const std::string meshfilename = inifile.get<std::string>("grid.gridFile");
if(verbose>0 && helper.rank()==0){
std::cout << "GRID SETUP:" << std::endl;
std::cout << " Reading Gmsh file " << meshfilename;
if(verbose<3) std::cout << " quietly" << std::endl;
else std::cout << std::endl << "--- GMSH READER OUTPUT ---" << std::endl;
}
Dune::GridFactory<GridType> factory;
if(helper.rank() == 0){
std::vector<int> boundary_index_map;
std::vector<int> element_index_map;
Dune::GmshReader<GridType>::read(factory,meshfilename,boundary_index_map,element_index_map,verbose>2?true:false,false);
}
auto grid = std::shared_ptr<GridType>(factory.createGrid());
grid->loadBalance();
grid->globalRefine(level);
if(verbose>0 && helper.rank()==0){
if (verbose>2) std::cout << "--- GMSH READER OUTPUT END ---" << std::endl;
std::cout << " Global Refinement level: " << level << std::endl;
if(verbose>1){
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: grid setup time" << std::setw(12) << timer.elapsed() << std::endl;
}
}
return grid;
}
/// Read parameter file specifications, create a grid, and return a shared pointer to the grid
/** Refine the grid globally and call default load balancing for multi-processing support.
* \tparam GridType Type of the grid. Must be supported by Dune::StructuredGridFactory
* \param inifile Parameter file parser
* \return Shared pointer to grid
*/
template<class GridType>
std::shared_ptr<GridType> build_grid_cube (const Dune::ParameterTree& inifile, const Dune::MPIHelper& helper)
{
check_parallel_allowed<GridType>(helper);
Dune::Timer timer;
ios_base_all_saver restore(std::cout);
enum { dim = GridType::dimension };
const int verbose = inifile.get<int>("output.verbose");
const int level = inifile.get<int>("grid.initialLevel");
const Dune::FieldVector<double,dim> lowerLeft(0.);
const Dune::FieldVector<double,dim> upperRight = inifile.get<Dune::FieldVector<double,dim>>("grid.extensions");
const std::array<unsigned int,dim> elements = inifile.get<std::array<unsigned int,dim>>("grid.cells");
auto grid = grid_cube_construction<GridType>(upperRight,elements);
grid->globalRefine(level);
grid->loadBalance();
if(verbose>0 && helper.rank()==0){
std::cout << "GRID SETUP:" << std::endl;
std::cout << " Creating Rectangular Grid in " << dim << " Dimensions." << std::endl;
std::cout << " Grid Extensions: ";
for (unsigned int i=0; i<dim; i++)
std::cout << upperRight[i] << " ";
std::cout << std::endl;
std::cout << " Grid Cells: ";
for (unsigned int i=0; i<dim; i++)
std::cout << elements[i] << " ";
std::cout << std::endl;
std::cout << " Global Refinement level: " << level << std::endl;
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: grid setup time" << std::setw(12) << timer.elapsed() << std::endl;
}
return grid;
}
return std::make_tuple(grid,
GridMapper(grid,
element_index_map,
boundary_index_map));
}
private:
/// Check if parallel execution is allowed
/** Currently, we only allow parallel execution on YaspGrid due to
* constraints imposed by the linear solver
*/
void check_parallel_allowed () const
{
if (_helper.size() > 1
&& not std::is_same_v<YaspGrid<Grid::dimension>, Grid>) {
DUNE_THROW(Dune::Exception,"DORiE does not support parallel "
"execution with this grid configuration!");
}
}
/// Create a GMSH grid from a file and the element and boundary index maps.
std::tuple<std::shared_ptr<Grid>, std::vector<int>, std::vector<int>>
create_gmsh_grid () const
{
Dune::Timer timer;
ios_base_all_saver restore(std::cout);
const int level = _config.get<int>("grid.initialLevel");
const std::string meshfilename = _config.get<std::string>("grid.gridFile");
if (_verbose > 0 && _helper.rank()==0){
std::cout << "GRID SETUP:" << std::endl;
std::cout << " Reading Gmsh file " << meshfilename;
if (_verbose < 3)
std::cout << " quietly" << std::endl;
else
std::cout << std::endl << "--- GMSH READER OUTPUT ---" << std::endl;
}
Dune::GridFactory<Grid> factory;
std::vector<int> boundary_index_map;
std::vector<int> element_index_map;
if (_helper.rank() == 0) {
Dune::GmshReader<Grid>::read(factory,
meshfilename,
boundary_index_map,
element_index_map,
_verbose > 2 ? true : false
);
}
auto grid = std::shared_ptr<Grid>(factory.createGrid());
grid->globalRefine(level);
if(_verbose>0 && _helper.rank()==0){
if (_verbose > 2)
std::cout << "--- GMSH READER OUTPUT END ---" << std::endl;
std::cout << " Global Refinement level: " << level << std::endl;
if (_verbose > 1) {
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: grid setup time" << std::setw(12) << timer.elapsed() << std::endl;
}
}
return std::make_tuple(grid, element_index_map, boundary_index_map);
}
/// Read parameter file specifications, create a grid, and return a shared pointer to the grid
std::tuple<std::shared_ptr<Grid>, std::vector<int>, std::vector<int>>
create_rect_grid () const
{
// check_parallel_allowed<GridType>(helper);
Dune::Timer timer;
ios_base_all_saver restore(std::cout);
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");
auto grid = grid_cube_construction<Grid>(upperRight, elements);
grid->globalRefine(level);
if(_verbose > 0 && _helper.rank()==0){
std::cout << "GRID SETUP:" << std::endl;
std::cout << " Creating Rectangular Grid in " << _dim << " Dimensions." << std::endl;
std::cout << " Grid Extensions: ";
for (unsigned int i=0; i<_dim; i++)
std::cout << upperRight[i] << " ";
std::cout << std::endl;
std::cout << " Grid Cells: ";
for (unsigned int i=0; i<_dim; i++)
std::cout << elements[i] << " ";
std::cout << std::endl;
std::cout << " Global Refinement level: " << level << std::endl;
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: grid setup time" << std::setw(12) << timer.elapsed() << std::endl;
}
// get the element index map from an H5 file
auto element_index_map = 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);
}
/// Default rectangular grid constructor. Calls StructuredGridFactory.
/** \param extensions Extensions of the grid
* \param elements Number of elements in each direction
* \return Shared pointer to the grid
*/
template<typename T=Grid>
std::enable_if_t<
not std::is_same_v<YaspGrid<_dim>, T>,
std::shared_ptr<Grid>>
grid_cube_construction (const Dune::FieldVector<double, _dim>& extensions,
const std::array<unsigned int, _dim>& elements) const
{
const Dune::FieldVector<double, _dim> origin(.0);
return Dune::StructuredGridFactory<Grid>::createCubeGrid(origin,
extensions,
elements);
}
/// Specialized constructor for parallel YaspGrid construction.
/** This is necessary because of a bug in StructuredGridFactory which does
* not create a grid with overlap.
* \todo The bug was fixed. This function can be removed after the next
* DUNE release.
* \param extensions Extensions of the grid
* \param elements Number of elements in each direction
* \return Shared pointer to the grid
*/
template<typename T=Grid>
std::enable_if_t<
std::is_same_v<YaspGrid<_dim>, T>,
std::shared_ptr<Grid>>
grid_cube_construction (const Dune::FieldVector<double, _dim>& extensions,
const std::array<unsigned int, _dim>& elements) const
{
std::array<int, _dim> elements_;
std::copy(elements.cbegin(),elements.cend(),elements_.begin());
return std::make_shared<YaspGrid<_dim>>(extensions, elements_);
}
/// Read the element mapping from an H5 file
std::vector<int> read_element_index_map_from_file () const
{
// read necessary values from config
const auto file_path = _config.get<std::string>("grid.mappingFile");
auto dataset_name = _config.get<std::string>("grid.mappingFileDataset");
// create H5 file reader
H5File file(file_path, _verbose);
// split group from dataset name if applicable, and open
const auto pos = dataset_name.rfind('/');
if (pos != std::string::npos) {
const auto group = dataset_name.substr(0, pos);
dataset_name = dataset_name.substr(pos+1);
file.open(group);
}
else {
file.open();
}
// Read the data
std::vector<int> element_index_map;
std::vector<int> cells;
file.read_dataset(element_index_map, cells, dataset_name, H5T_NATIVE_INT);
// verify number of cells
const auto elements = _config.get<std::array<unsigned int, _dim>>("grid.cells");
for (size_t i = 0; i < _dim; ++i) {
assert(cells[i] == elements[i]);
}
return element_index_map;
}
};
} // namespace Dorie
} // namespace Dune
......
dune/dorie/solver/util_grid_mapper.hh 0 → 100644
View file @ 0507c49d
#ifndef DUNE_DORIE_UTIL_GRID_MAPPER_HH
#define DUNE_DORIE_UTIL_GRID_MAPPER_HH
#include <vector>
#include <memory>
#include <map>
/// A class for managing local grid mappings and data load balancing.
template<typename Grid>
class GridMapper
{
public:
using GridView = typename Grid::LevelGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Index = typename Mapper::Index;
using IdSet = typename Grid::GlobalIdSet;
using IdType = typename IdSet::IdType;
private:
//! A shared pointer to the grid
std::shared_ptr<Grid> _grid;
//! The grid view to operate on.
GridView _gv;
//! The mapper for local grid indices of entities
Mapper _mapper;
//! The global persistent ID set of grid entities.
const IdSet& _idset;
//! Maps the index of a grid entity to an arbitrary index.
std::vector<int> _element_index_map;
//! Maps the index of a grid boundary to an arbitrary index.
std::vector<int> _boundary_index_map;
public:
/// Create the grid mapper from a *unbalanced* grid.
/** All data structures received must be located on a single process.
* After constructing this GridMapper, the grid and the maps are properly
* partitioned and can be retrieved using the getter functions.
* \param grid Shared pointer to the grid.
* \param element_index_map The element index map for the entire grid.
* \param boundary_index_map The boundary index map for the entire grid.
*/
explicit GridMapper (const std::shared_ptr<Grid> grid,
const std::vector<int>& element_index_map,
const std::vector<int>& boundary_index_map) :
_grid(grid),
_gv(_grid->levelGridView(0)),
_mapper(_gv, Dune::mcmgElementLayout()),
_idset(_grid->globalIdSet()),
_element_index_map(element_index_map),
_boundary_index_map(boundary_index_map)
{
// create a mapping using global persistent IDs
auto element_id_map = build_element_id_map();
// broadcast this map
broadcast_element_id_map(element_id_map);
// load-balance the grid to all processors
_grid->loadBalance();
_mapper.update();
// rebuild the index map for local grid partition
rebuild_element_index_map(element_id_map);
}
/// Return a reference to the element index map
const std::vector<int>& get_element_index_map () const
{
return _element_index_map;
}
/// Return a reference to the boundary index map
const std::vector<int>& get_boundary_index_map () const
{
return _boundary_index_map;
}
private:
/// Build a (global) id map from an element index map
/** \todo Ensure that grid is not load balanced
* \return Global ID map with the same mapping as element_index_map.
*/
std::map<IdType, int> build_element_id_map () const
{
// iterate over grid elements
std::map<IdType, int> ret;
for (auto&& elem : elements(_gv)) {
const auto index = _mapper.index(elem);
const auto id = _idset.id(elem);
ret.emplace(id, _element_index_map.at(index));
}
return ret;
}
/// Build a element_index_map from a element_id_map
/** \param element_id_map The ID map to build the index map from
*/
void rebuild_element_index_map (const std::map<IdType, int>& element_id_map)
{
_element_index_map.resize(_mapper.size());
for (auto&& elem : elements(_gv)) {
const auto index = _mapper.index(elem);
const auto id = _idset.id(elem);
_element_index_map.at(index) = element_id_map.at(id);
}
}
/// Broadcast the element id map to all processors
/** \param element_id_map The map to be broadcasted (and hence manipulated)
*/
void broadcast_element_id_map (std::map<IdType, int>& element_id_map) const
{
// broadcast the size of the map
auto size = element_id_map.size();
MPI_Bcast(&size, 1,
MPI_UNSIGNED_LONG, 0, Dune::MPIHelper::getCommunicator()
);
// build two vectors from keys and valuses and fill them
std::vector<IdType> ids;
std::vector<int> values;
for (auto&& [id, value] : element_id_map) {
ids.push_back(id);
values.push_back(value);
}
// extend to correct size on all processes
ids.resize(size);
values.resize(size);
// broadcast the data
MPI_Bcast(ids.data(), size,
MPI_LONG, 0, Dune::MPIHelper::getCommunicator()
);
MPI_Bcast(values.data(), size,
MPI_INT, 0, Dune::MPIHelper::getCommunicator()
);
// clear the map and fill it with the values for this process
element_id_map.clear();
std::transform(ids.begin(), ids.end(), values.begin(),
std::inserter(element_id_map, element_id_map.end()),
[](const auto key, const auto val) {
return std::make_pair(key, val);
}
);
}
};
#endif // DUNE_DORIE_UTIL_GRID_MAPPER_HH
\ No newline at end of file
dune/dorie/solver/util_h5tools.hh
View file @ 0507c49d
......@@ -36,7 +36,7 @@ namespace Dune {
DUNE_THROW(Exception,"File '" << filePath << "' is not readable");
}
void open(const std::string& group_name)
void open(const std::string& group_name="./")
{
if(fileOpen)
close();
......@@ -68,9 +68,12 @@ namespace Dune {
*
*/
template<typename VEC, typename CELLVEC>
void read_dataset(VEC& local_data, CELLVEC& local_cells, const std::string& dataset_name)
void read_dataset(VEC& local_data,
CELLVEC& local_cells,
const std::string& dataset_name,
const hid_t data_type)
{
if(!fileOpen)
if (not fileOpen)
DUNE_THROW(Exception,"File must be opened before reading data");
// open the dataset
......@@ -89,31 +92,29 @@ namespace Dune {
arr_dim = H5Sget_simple_extent_ndims(dataspace_id);
// get the size of the problem
std::vector<hsize_t> dims;
dims.resize(arr_dim);
status = H5Sget_simple_extent_dims(dataspace_id , &(dims[0]), 0);
std::vector<hsize_t> dims(arr_dim, 0);
status = H5Sget_simple_extent_dims(dataspace_id , dims.data(), 0);
assert(status > -1);
typedef std::vector<hsize_t>::size_type Index;
Index local_size=1;
hsize_t offset[arr_dim];
for(Index i=0; i<arr_dim; i++)
{
local_cells[arr_dim-i-1]=dims[i];
local_size*=dims[i];
offset[i]=0;
}