Skip to content

GitLab

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

Overhaul grid creation and mapping process 

* GridMapper is now called internally by the constructor of GridCreator.
    The object is not returned by the creator anymore.
* Specialize GridMappers for 'gmsh' and 'rectangular' grids and
    derive them from GridMapperBase
* GridMapper now reads the mapping H5 file
* The grid is now retrieved from GridCreator
* Local element and boundary index maps are now retrieved from
    GridCreator
* Adapt simulations to expect GridCreator as argument

Minor changes:
* Store file name in H5File
* Add free function for lowercasing strings
* Add logger to GridMapper
parent fbb7f2fa
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 484 additions and 398 deletions
dune/dorie/common/grid_creator.hh
View file @ d99c069a
......@@ -5,7 +5,7 @@
#include <array>
#include <vector>
#include <iostream>
#include <iomanip>
#include <iomanip>
#include <type_traits>
#include <map>
#include <string>
......@@ -31,11 +31,6 @@
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
......@@ -55,10 +50,18 @@ 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:
explicit GridCreator (const Dune::ParameterTree& config,
......@@ -67,38 +70,38 @@ public:
_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.get_element_index_map();
_boundary_index_map = mapper.get_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.get_element_index_map();
_boundary_index_map = mapper.get_boundary_index_map();
}
else {
_log->error("DORiE expects to work on UGGrid for "
......@@ -107,18 +110,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& get_element_index_map () const { return _element_index_map; }
/// Return the local boundary index map
const Map& get_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,7 +180,7 @@ 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;
......@@ -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())
{
const auto global_id = _config.get<int>("grid.globalIndex");
const auto total_size = std::accumulate(elements.begin(),
elements.end(),
1,
std::multiplies<int>());
// fill index map with the global id and copy over the cells
_log->debug("Mapping global medium index '{}' to entire grid",
global_id);
element_index_map.resize(total_size, global_id);
std::copy(elements.begin(), elements.end(), cells.begin());
return std::make_tuple(element_index_map, cells);
}
// load heterogeneous mapping
auto dataset_name = _config.get<std::string>("grid.mappingFileDataset");
_log->debug("Loading medium mapping from file. "
"Filename: {}, Dataset: {}",
file_path, dataset_name);
// create H5 file reader
H5File file(file_path);
// Read the data
file.read_dataset(dataset_name,
H5T_NATIVE_INT,
element_index_map,
cells
);
// shape has to be reversed
std::reverse(begin(cells), end(cells));
// check if cells has correct size and extensions
if (cells.size() != Grid::dimension) {
_log->error("Mapping dataset dimensions do not match "
"grid dimensions. Expected: {}, Read: {}",
Grid::dimension, cells.size());
DUNE_THROW(IOError, "Mapping dataset has wrong dimensions!");
}
for (size_t i = 0; i < _dim; ++i) {
if (cells[i] != elements[i]) {
_log->error("Mapping dataset extensions do not match grid "
"extensions. Expected: {}, Read: {}",
to_string(elements), to_string(cells));
DUNE_THROW(IOError, "Mapping dataset extensions don't match "
" the grid");
}
}
return std::make_tuple(element_index_map, cells);
}
};
} // namespace Dorie
......
dune/dorie/common/grid_mapper.hh
View file @ d99c069a
This diff is collapsed.
dune/dorie/common/h5file.hh
View file @ d99c069a
......@@ -23,6 +23,8 @@ private:
using Logger = std::shared_ptr<spdlog::logger>;
//! The logger of this instance
const Logger _log;
//! The name of this file
const std::string _file_path;
//! ID of the H5 file
hid_t _file_id;
//! ID of the currently opened group
......@@ -38,6 +40,7 @@ public:
explicit H5File(const std::string& file_path,
const Logger log=get_logger(log_base)):
_log(log),
_file_path(file_path),
_file_id(-1),
_group_id(-1)
{
......@@ -80,6 +83,9 @@ public:
assert(status > -1);
}
/// Return the file path of this object
std::string path () const { return _file_path; }
/// Open a group
/** \param group_path Full internal path to the group.
* Defaults to the base group.
......@@ -156,7 +162,7 @@ public:
dims.data(),
0);
assert(status > -1);
// log dataset extensions as read
_log->debug(" Dataset extensions: {}", to_string(dims));
......
dune/dorie/dorie.cc
View file @ d99c069a
......@@ -90,7 +90,6 @@ int main(int argc, char** argv)
{
if (gtype == "gmsh"){
Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
case 1:{
Sim<Simplex<2,1>> sim(richards_config, grid_mapper, helper);
......@@ -117,7 +116,6 @@ int main(int argc, char** argv)
else if (gtype == "rectangular"){
if (adapt_policy != Dune::Dorie::AdaptivityPolicy::None) {
Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
case 1:{
Sim<CubeAdaptive<2,1>> sim(richards_config, grid_mapper, helper);
......@@ -143,7 +141,6 @@ int main(int argc, char** argv)
}
else{ // no adaptivity
Dune::Dorie::GridCreator<Dune::YaspGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
case 1:{
Sim<Cube<2,1>> sim(richards_config, grid_mapper, helper);
......@@ -173,7 +170,6 @@ int main(int argc, char** argv)
{
if (gtype == "gmsh"){
Dune::Dorie::GridCreator<Dune::UGGrid<3>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
case 1:{
Sim<Simplex<3,1>> sim(richards_config, grid_mapper, helper);
......@@ -200,7 +196,6 @@ int main(int argc, char** argv)
else if (gtype == "rectangular"){
if(adapt_policy != Dune::Dorie::AdaptivityPolicy::None){
Dune::Dorie::GridCreator<Dune::UGGrid<3>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
case 1:{
Sim<CubeAdaptive<3,1>> sim(richards_config, grid_mapper, helper);
......@@ -226,7 +221,6 @@ int main(int argc, char** argv)
}
else{ // no adaptivity
Dune::Dorie::GridCreator<Dune::YaspGrid<3>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
case 1:{
Sim<Cube<3,1>> sim(richards_config, grid_mapper, helper);
......
dune/dorie/model/richards/richards.cc
View file @ d99c069a
......@@ -6,7 +6,7 @@ namespace Dorie{
template<typename Traits>
RichardsSimulation<Traits>::RichardsSimulation (
Dune::ParameterTree& _inifile,
std::shared_ptr<GridMapper> _grid_mapper,
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
):
SimulationBase(log_richards,
......@@ -16,7 +16,7 @@ RichardsSimulation<Traits>::RichardsSimulation (
inifile(_inifile),
output_type(_inifile.get<bool>("output.asciiVtk") ?
Dune::VTK::OutputType::ascii : Dune::VTK::OutputType::base64),
grid(_grid_mapper->get_grid()),
grid(_grid_creator.grid()),
gv(grid->leafGridView()),
mbe_slop(estimate_mbe_entries<typename MBE::size_type>(
Traits::dim,Traits::GridGeometryType)),
......@@ -32,7 +32,7 @@ RichardsSimulation<Traits>::RichardsSimulation (
Dune::PDELab::constraints(*gfs,*cc,false);
// --- Create the new parameter class
auto element_map = _grid_mapper->get_element_index_map();
auto element_map = _grid_creator.get_element_index_map();
fparam = std::make_shared<FlowParameters>(inifile, grid, element_map);
// --- Operator Helper Classes ---
......
dune/dorie/model/richards/richards.hh
View file @ d99c069a
......@@ -17,7 +17,7 @@
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/interpolator.hh>
#include <dune/dorie/common/time_controller.hh>
#include <dune/dorie/common/grid_mapper.hh>
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/common/grid_function_writer.hh>
#include <dune/dorie/model/richards/adaptivity/handler.hh>
......@@ -73,7 +73,7 @@ struct RichardsSimulationTraits : public BaseTraits
// -- DORiE Class Definitions -- //
/// Wrapper for grid and volume & boundary mappings
using GridMapper = Dune::Dorie::GridMapper<typename BaseTraits::Grid>;
using GridCreator = Dune::Dorie::GridCreator<typename BaseTraits::Grid>;
/// Class defining the soil parameterization
using FlowParameters = Dune::Dorie::FlowParameters<BaseTraits>;
/// Class defining boundary conditions
......@@ -173,7 +173,7 @@ private:
// -- DORiE Class Definitions -- //
/// Wrapper for grid and volume & boundary mappings
using GridMapper = typename Traits::GridMapper;
using GridCreator = typename Traits::GridCreator;
/// Class defining the soil parameterization
using FlowParameters = typename Traits::FlowParameters;
/// Class defining boundary conditions
......@@ -290,7 +290,7 @@ public:
*/
RichardsSimulation (
Dune::ParameterTree& _inifile,
std::shared_ptr<GridMapper> _grid_mapper,
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
);
......
dune/dorie/model/transport/transport.cc
View file @ d99c069a
......@@ -6,13 +6,13 @@ namespace Dorie{
template<typename Traits>
TransportSimulation<Traits>::TransportSimulation(
Dune::ParameterTree& _inifile,
std::shared_ptr<GridMapper> _grid_mapper,
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
) : SimulationBase(log_transport,
_inifile.get<std::string>("output.logLevel"),
_helper)
, helper(_helper)
, grid(_grid_mapper->get_grid())
, grid(_grid_creator.grid())
, inifile(_inifile)
, gv(grid->leafGridView())
, mbe_slop(estimate_mbe_entries<typename MBE::size_type>(Traits::dim,Traits::GridGeometryType))
......
dune/dorie/model/transport/transport.hh
View file @ d99c069a
......@@ -74,7 +74,7 @@ struct TransportSimulationTraits : public BaseTraits
// -- DORiE Class Definitions -- //
/// Wrapper for grid and volume & boundary mappings
using GridMapper = Dune::Dorie::GridMapper<typename BaseTraits::Grid>;
using GridCreator = Dune::Dorie::GridCreator<typename BaseTraits::Grid>;
/// Class defining boundary conditions
using SoluteBoundary = Dune::Dorie::SoluteBoundary<BaseTraits>;
/// Class defining the initial condition
......@@ -166,7 +166,7 @@ protected:
// -- D/*ORiE Class Definitions -- //
/// Wrapper for grid and volume & boundary mappings
using GridMapper = typename Traits::GridMapper;
using GridCreator = typename Traits::GridCreator;
/// Class defining boundary conditions
using SoluteBoundary = typename Traits::SoluteBoundary;
/// Class defining the initial condition
......@@ -263,7 +263,7 @@ public:
*/
TransportSimulation (
Dune::ParameterTree& _inifile,
std::shared_ptr<GridMapper> _grid_mapper,
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
);
......
dune/dorie/test/test-grid-creation.cc
View file @ d99c069a
......@@ -24,6 +24,9 @@ int main(int argc, char** argv)
auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv);
log.reset(); // silence compiler warning
if (inifile.get<bool>("misc.debugMode"))
Dune::Dorie::Setup::debug_hook();
const auto dimensions = inifile.get<int>("grid.dimensions");
if (dimensions == 2) {
constexpr int dim = 2;
......
dune/dorie/test/test-grid-creation.hh
View file @ d99c069a
......@@ -21,11 +21,10 @@ void create_grid_and_test_mapping (const Dune::ParameterTree& config,
{
// create the grid mapper
Dune::Dorie::GridCreator<Grid> gc(config, helper);
auto grid_mapper = gc.get_mapper();
// get grid and data
auto grid = grid_mapper->get_grid();
auto element_map = grid_mapper->get_element_index_map();
auto grid = gc.grid();
auto element_map = gc.get_element_index_map();
// check mapping
auto gv = grid->levelGridView(0);
......
dune/dorie/test/test-parameterization-scaled.cc
View file @ d99c069a
......@@ -76,9 +76,8 @@ int main (int argc, char** argv)
// create the grid
Dune::Dorie::GridCreator<Traits<2>::Grid> gc(inifile, helper);
auto mapper = gc.get_mapper();
auto grid = mapper->get_grid();
const auto index_map = mapper->get_element_index_map();
auto grid = gc.grid();
const auto index_map = gc.get_element_index_map();
inifile = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
......
dune/dorie/test/test-parameterization.cc
View file @ d99c069a
......@@ -164,9 +164,8 @@ int main (int argc, char** argv)
// create the grid
Dune::Dorie::GridCreator<Traits<2>::Grid> gc(inifile, helper);
auto mapper = gc.get_mapper();
auto grid = mapper->get_grid();
const auto index_map = mapper->get_element_index_map();
auto grid = gc.grid();
const auto index_map = gc.get_element_index_map();
inifile = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
......
test/test-mass-conservation.cc
View file @ d99c069a
......@@ -62,8 +62,7 @@ int main (int argc, char** argv)
using TestSim = Dune::Dorie::TestSimulation<CubeAdaptive<2, 1>>;
Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
TestSim sim(richards_config, grid_mapper, helper);
TestSim sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
result = sim.run_test();
......@@ -71,9 +70,8 @@ int main (int argc, char** argv)
else {
using TestSim = Dune::Dorie::TestSimulation<Cube<2, 1>>;
Dune::Dorie::GridCreator<Dune::YaspGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
TestSim sim(richards_config, grid_mapper, helper);
TestSim sim(richards_config, grid_creator, helper);
result = sim.run_test();
}
}
......@@ -81,9 +79,8 @@ int main (int argc, char** argv)
else if (grid_type == "gmsh") {
using TestSim = Dune::Dorie::TestSimulation<Simplex<2, 1>>;
Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
TestSim sim(richards_config, grid_mapper, helper);
TestSim sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
result = sim.run_test();
}
......
test/test-mass-conservation.hh
View file @ d99c069a
......@@ -191,7 +191,7 @@ protected:
using Grid = typename Traits::Grid;
using U = typename Traits::U;
using GridMapper = typename Traits::GridMapper;
using GridCreator = typename Traits::GridCreator;
using FlowParameters = typename Traits::FlowParameters;
using FlowBoundary = typename Traits::FlowBoundary;
......@@ -217,9 +217,9 @@ public:
*/
TestSimulation(
Dune::ParameterTree& _inifile,
std::shared_ptr<GridMapper> _grid_mapper,
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
) : RichardsSimulation<Traits>(_inifile, _grid_mapper, _helper)
) : RichardsSimulation<Traits>(_inifile, _grid_creator, _helper)
, acc(0.)
, acc_square(0.)
{
......
test/test-simulation-dummy-transport.cc
View file @ d99c069a
......@@ -92,10 +92,9 @@ int main(int argc, char** argv) {
if (gtype == "rectangular")
{
Dune::Dorie::GridCreator<Dune::YaspGrid<2>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
using Traits = Cube<2>;
Sim<Traits> sim(inifile_transport,grid_mapper,helper);
Sim<Traits> sim(inifile_transport,grid_creator,helper);
typename Traits::Vector flux;
typename Traits::Scalar sat;
......@@ -106,7 +105,7 @@ int main(int argc, char** argv) {
using GFFlux = typename Traits::GFWaterFlux;
using GFSat = typename Traits::GFWaterContent;
auto grid = grid_mapper->get_grid();
auto grid = grid_creator.grid();
auto gf_flux = std::make_shared<GFFlux>(grid->leafGridView(),flux);
auto gf_theta = std::make_shared<GFSat>(grid->leafGridView(),sat);
......@@ -126,10 +125,9 @@ int main(int argc, char** argv) {
if (gtype == "rectangular")
{
Dune::Dorie::GridCreator<Dune::YaspGrid<3>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
using Traits = Cube<3>;
Sim<Traits> sim(inifile_transport,grid_mapper,helper);
Sim<Traits> sim(inifile_transport,grid_creator,helper);
typename Traits::Vector flux;
typename Traits::Scalar theta;
......@@ -137,7 +135,7 @@ int main(int argc, char** argv) {
flux[dim-1] = inifile.get<double>("_constantWater.flux",-5.5e-8);
theta = inifile.get<double>("_constantWater.waterContent",1.);
auto grid = grid_mapper->get_grid();
auto grid = grid_creator.grid();
using GFFlux = typename Traits::GFWaterFlux;
using GFSat = typename Traits::GFWaterContent;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment