Skip to content

Commit

Permalink
Allow data sharing in StructuredDataLookup::reinit().
Browse files Browse the repository at this point in the history
  • Loading branch information
@bangerth
bangerth committed Jun 29, 2021
1 parent 84876b9 commit b9589aa
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 39 deletions.
37 changes: 32 additions & 5 deletions include/aspect/structured_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,14 +87,41 @@ namespace aspect
*
* The data in @p data_table consists of a Table for each of the @p n_components components.
*
* The last two arguments are rvalue references, and the function will
* move the data so provided into another storage location. In other
* words, after the call, the variables passed as the last two
* arguments may be empty or otherwise altered.
* The `coordinate_values` and `data_table` arguments are rvalue references,
* and the function will move the data so provided into another storage
* location. In other words, after the call, the variables passed as the
* last two arguments may be empty or otherwise altered.
*
#if DEAL_II_VERSION_GTE(9,3,0)
* If ASPECT is built on deal.II version 9.3 or higher, this class
* is able to share data between processes located on the same
* machine. In this case, if the last argument, `root_process` is
* set to anything other than `numbers::invalid_unsigned_int`, then only the
* indicated root process within the given MPI communicator needs to
* pass in valid arguments whereas on all other processes the values
* of `column_names`, `coordinate_values`, and `data_table` are
* ignored. Instead, the indicated root process will make sure that
* the other processes obtain valid data, for example by sharing
* data tables in shared memory spaces. This reduces both the
* computational effort in reading data from disk and parsing it on
* all processes, and can vastly reduce the amount of memory required
* on machines where many processor cores have access to shared memory
* resources.
*
* If `root_process` equals `numbers::invalid_unsigned_int`, then
* every process needs to pass data for all arguments and the
* `mpi_communicator` argument is ignored.
#else
* If ASPECT is built on a version of deal.II older than 9.3,
* then all processes need to pass in valid values for the first
* three arguments and the last two arguments are ignored.
#endif
*/
void reinit(const std::vector<std::string> &column_names,
std::vector<std::vector<double>> &&coordinate_values,
std::vector<Table<dim,double> > &&data_table);
std::vector<Table<dim,double> > &&data_table,
const MPI_Comm &mpi_communicator,
const unsigned int root_process);

/**
* Loads a data text file. Throws an exception if the file does not
Expand Down
118 changes: 87 additions & 31 deletions source/structured_data.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,43 +169,97 @@ namespace aspect
void
StructuredDataLookup<dim>::reinit(const std::vector<std::string> &column_names,
std::vector<std::vector<double>> &&coordinate_values_,
std::vector<Table<dim,double> > &&data_table)
std::vector<Table<dim,double> > &&data_table,
const MPI_Comm &mpi_communicator,
const unsigned int root_process)
{
Assert(coordinate_values_.size()==dim, ExcMessage("Invalid size of coordinate_values."));
for (unsigned int d=0; d<dim; ++d)
#if DEAL_II_VERSION_GTE(9,3,0)
const bool supports_shared_data = true;
#else
const bool supports_shared_data = false;
#endif

// If this is the root process, or if we don't support data sharing,
// then set up the various member variables we need to compute
// from the input data
if ((supports_shared_data == false)
||
((supports_shared_data == true)
&&
(root_process != numbers::invalid_unsigned_int)
&&
(Utilities::MPI::this_mpi_process(mpi_communicator) == root_process)))
{
this->coordinate_values[d] = std::move(coordinate_values_[d]);
AssertThrow(this->coordinate_values[d].size()>1,
ExcMessage("Error: At least 2 entries per coordinate direction are required."));
table_points[d] = this->coordinate_values[d].size();
Assert(coordinate_values_.size()==dim, ExcMessage("Invalid size of coordinate_values."));
for (unsigned int d=0; d<dim; ++d)
{
this->coordinate_values[d] = std::move(coordinate_values_[d]);
AssertThrow(this->coordinate_values[d].size()>1,
ExcMessage("Error: At least 2 entries per coordinate direction are required."));
table_points[d] = this->coordinate_values[d].size();
}

components = column_names.size();
data_component_names = column_names;
Assert(data_table.size() == components,
ExcMessage("Error: Incorrect number of columns specified."));
for (unsigned int c=0; c<components; ++c)
Assert(data_table[c].size() == table_points,
ExcMessage("Error: One of the data tables has an incorrect size."));

// compute maximum_component_value for each component:
maximum_component_value = std::vector<double>(components,-std::numeric_limits<double>::max());
for (unsigned int c=0; c<components; ++c)
{
const unsigned int n_elements = data_table[c].n_elements();
for (unsigned int idx=0; idx<n_elements; ++idx)
maximum_component_value[c] = std::max(maximum_component_value[c], data_table[c](
compute_table_indices(table_points, idx)));
}

// In case the data is specified on a grid that is equidistant
// in each coordinate direction, we only need to store
// (besides the data) the number of intervals in each direction and
// the begin- and end-points of the coordinates.
// In case the grid is not equidistant, we need to keep
// all the coordinates in each direction, which is more costly.
coordinate_values_are_equidistant = data_is_equidistant<dim> (coordinate_values);
}

components = column_names.size();
data_component_names = column_names;
Assert(data_table.size() == components,
ExcMessage("Error: Incorrect number of columns specified."));
for (unsigned int c=0; c<components; ++c)
Assert(data_table[c].size() == table_points,
ExcMessage("Error: One of the data tables has an incorrect size."));

// compute maximum_component_value for each component:
maximum_component_value = std::vector<double>(components,-std::numeric_limits<double>::max());
for (unsigned int c=0; c<components; ++c)
// If deal.II is new enough to support sharing data, and if the
// caller of this function has actually requested this, then we have
// set up member variables on the root process, but not on any of
// the other processes. Broadcast the data to the remaining
// processes
if ((supports_shared_data == true)
&&
(root_process != numbers::invalid_unsigned_int))
{
const unsigned int n_elements = data_table[c].n_elements();
for (unsigned int idx=0; idx<n_elements; ++idx)
maximum_component_value[c] = std::max(maximum_component_value[c], data_table[c](
compute_table_indices(table_points, idx)));
#if DEAL_II_VERSION_GTE(9,3,0)
coordinate_values = Utilities::MPI::broadcast (mpi_communicator,
coordinate_values,
root_process);
components = Utilities::MPI::broadcast (mpi_communicator,
components,
root_process);
data_component_names = Utilities::MPI::broadcast (mpi_communicator,
data_component_names,
root_process);
maximum_component_value = Utilities::MPI::broadcast (mpi_communicator,
maximum_component_value,
root_process);
coordinate_values_are_equidistant = Utilities::MPI::broadcast (mpi_communicator,
coordinate_values_are_equidistant,
root_process);

// We can then also prepare the data tables for sharing between
// processes
for (unsigned int c = 0; c < components; ++c)
data_table[c].replicate_across_communicator (mpi_communicator,
root_process);
#endif
}

// In case the data is specified on a grid that is equidistant
// in each coordinate direction, we only need to store
// (besides the data) the number of intervals in each direction and
// the begin- and end-points of the coordinates.
// In case the grid is not equidistant, we need to keep
// all the coordinates in each direction, which is more costly.
coordinate_values_are_equidistant = data_is_equidistant<dim> (coordinate_values);

// For each data component, set up a GridData,
// its type depending on the read-in grid.
data.resize(components);
Expand Down Expand Up @@ -437,7 +491,9 @@ namespace aspect
// objects.
this->reinit(column_names,
std::move(coordinate_values),
std::move(data_tables));
std::move(data_tables),
comm,
numbers::invalid_unsigned_int);
}


Expand Down
9 changes: 6 additions & 3 deletions unit_tests/utilities.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,8 @@ TEST_CASE("Utilities::AsciiDataLookup manual dim=1")
raw_data[1](0) = 5.0;
raw_data[1](1) = 3.0;

lookup.reinit(column_names, std::move(coordinate_values), std::move(raw_data));
lookup.reinit(column_names, std::move(coordinate_values), std::move(raw_data),
MPI_COMM_SELF, numbers::invalid_unsigned_int);

INFO(lookup.get_data(Point<1>(1.5), 0));
INFO(lookup.get_data(Point<1>(1.5), 1));
Expand Down Expand Up @@ -111,7 +112,8 @@ TEST_CASE("Utilities::AsciiDataLookup manual dim=2")
raw_data[0](1,2) = 0.0;
raw_data[0](2,2) = 0.0;

lookup.reinit(column_names, std::move(coordinate_values), std::move(raw_data));
lookup.reinit(column_names, std::move(coordinate_values), std::move(raw_data),
MPI_COMM_SELF, numbers::invalid_unsigned_int);

REQUIRE(lookup.get_data(Point<2>(1.0,6.0),0) == Approx(5.0));
REQUIRE(lookup.get_data(Point<2>(2.0,6.0),0) == Approx(5.5));
Expand Down Expand Up @@ -142,7 +144,8 @@ TEST_CASE("Utilities::AsciiDataLookup manual dim=2 equid")
raw_data[0](1,2) = 0.0;
raw_data[0](2,2) = 0.0;

lookup.reinit(column_names, std::move(coordinate_values), std::move(raw_data));
lookup.reinit(column_names, std::move(coordinate_values), std::move(raw_data),
MPI_COMM_SELF, numbers::invalid_unsigned_int);

REQUIRE(lookup.get_data(Point<2>(1.0,6.0),0) == Approx(5.0));
REQUIRE(lookup.get_data(Point<2>(1.5,6.0),0) == Approx(5.5));
Expand Down

0 comments on commit b9589aa

Please sign in to comment.