Select Git revision
CRSMatrix.hpp
-
Stéphane Del Pino authored
The new class is meant to be much faster - it requires an estimation of the number values per row - if the number of stored values per line exceeds the reserved one, an a special storage zone (per row) is used. It is important to give a precise estimation of the number of values per rows: - if the provided expected number of values is too small with regard to the real number of values stored for a given row, then building performances can be noticeably affected - if the provided expected number of values per row corresponds exactly to the effective number of stored values, then the construction of the CRSMatrix costs nothing. In the later case, one benefits from a significant CPU and memory costs reduction.
Stéphane Del Pino authoredThe new class is meant to be much faster - it requires an estimation of the number values per row - if the number of stored values per line exceeds the reserved one, an a special storage zone (per row) is used. It is important to give a precise estimation of the number of values per rows: - if the provided expected number of values is too small with regard to the real number of values stored for a given row, then building performances can be noticeably affected - if the provided expected number of values per row corresponds exactly to the effective number of stored values, then the construction of the CRSMatrix costs nothing. In the later case, one benefits from a significant CPU and memory costs reduction.
ConnectivityDispatcher.cpp 40.22 KiB
#include <mesh/ConnectivityDispatcher.hpp>
#include <mesh/ItemOfItemType.hpp>
#include <utils/CRSGraph.hpp>
#include <utils/Partitioner.hpp>
#include <iostream>
#warning remove
#include <utils/Demangle.hpp>
template <size_t Dimension>
template <ItemType item_type>
void
ConnectivityDispatcher<Dimension>::_buildNewOwner()
{
std::cout << "#### " << rang::fgB::yellow << itemName(item_type) << rang::fg::reset << " _buildNewOwner() ####\n";
if constexpr (item_type == ItemType::cell) {
CRSGraph connectivity_graph = m_connectivity.ownCellToCellGraph();
Partitioner P;
Array new_owner_array = P.partition(connectivity_graph);
std::cout << "new_owner_array = " << new_owner_array << '\n';
const auto& cell_is_owned = m_connectivity.cellIsOwned();
std::cout << "cell_is_owned = " << cell_is_owned << '\n';
CellValue<int> cell_new_owner(m_connectivity);
size_t i = 0;
for (CellId cell_id = 0; cell_id < m_connectivity.numberOfCells(); ++cell_id) {
if (cell_is_owned[cell_id]) {
cell_new_owner[cell_id] = new_owner_array[i++];
}
}
synchronize(cell_new_owner);
std::cout << "cell_new_owner = " << cell_new_owner << '\n';
std::cout << "cell_number = " << m_connectivity.template number<item_type>() << '\n';
this->_dispatchedInfo<ItemType::cell>().m_new_owner = cell_new_owner;
} else {
const auto& item_to_cell_matrix = m_connectivity.template getItemToItemMatrix<item_type, ItemType::cell>();
const auto& cell_number = m_connectivity.cellNumber();
const auto& cell_new_owner = this->_dispatchedInfo<ItemType::cell>().m_new_owner;
using ItemId = ItemIdT<item_type>;
ItemValue<int, item_type> item_new_owner(m_connectivity);
parallel_for(
item_new_owner.numberOfItems(), PUGS_LAMBDA(const ItemId& l) {
const auto& item_to_cell = item_to_cell_matrix[l];
CellId Jmin = item_to_cell[0];
for (size_t j = 1; j < item_to_cell.size(); ++j) {
const CellId J = item_to_cell[j];
if (cell_number[J] < cell_number[Jmin]) {
Jmin = J;
}
}
item_new_owner[l] = cell_new_owner[Jmin];
});
synchronize(item_new_owner);
std::cout << itemName(item_type) << "_new_owner = " << item_new_owner << '\n';
std::cout << itemName(item_type) << "_number = " << m_connectivity.template number<item_type>() << '\n';
this->_dispatchedInfo<item_type>().m_new_owner = item_new_owner;
}
}
template <size_t Dimension>
template <ItemType item_type>
void
ConnectivityDispatcher<Dimension>::_buildItemToExchangeLists()
{
this->_buildItemListToSend<item_type>();
this->_buildNumberOfItemToExchange<item_type>();
if constexpr (item_type == ItemType::cell) {
this->_buildCellNumberIdMap();
}
this->_buildRecvItemIdCorrespondanceByProc<item_type>();
}
template <size_t Dimension>
template <ItemType item_type>
void
ConnectivityDispatcher<Dimension>::_buildItemListToSend()
{
if constexpr (item_type == ItemType::cell) {
const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
const auto& cell_to_node_matrix = m_connectivity.cellToNodeMatrix();
const auto& cell_is_owned = m_connectivity.cellIsOwned();
const auto& cell_new_owner = this->_dispatchedInfo<ItemType::cell>().m_new_owner;
const auto& cell_number = m_connectivity.cellNumber();
std::cout << "- cell_new_owner = " << cell_new_owner << '\n';
std::vector<std::vector<CellId>> cell_vector_to_send_by_proc(parallel::size());
Array<bool> send_to_rank(parallel::size());
for (CellId cell_id = 0; cell_id < m_connectivity.numberOfCells(); ++cell_id) {
send_to_rank.fill(false);
const auto& cell_to_node = cell_to_node_matrix[cell_id];
std::cout << rang::fgB::cyan << " looking for cell " << cell_id << "(" << cell_number[cell_id] << ")"
<< rang::fg::reset << "\n";
if (cell_is_owned[cell_id]) {
for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
const NodeId& node_id = cell_to_node[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t i_node_cell = 0; i_node_cell < node_to_cell.size(); ++i_node_cell) {
const CellId& node_cell = node_to_cell[i_node_cell];
send_to_rank[cell_new_owner[node_cell]] = true;
std::cout << rang::fgB::green << " storing cell " << node_cell << "(" << cell_number[node_cell] << ")"
<< rang::fg::reset << "\n";
}
}
for (size_t i_rank = 0; i_rank < send_to_rank.size(); ++i_rank) {
if (send_to_rank[i_rank]) {
cell_vector_to_send_by_proc[i_rank].push_back(cell_id);
}
}
}
}
auto& cell_list_to_send_by_proc = this->_dispatchedInfo<ItemType::cell>().m_list_to_send_by_proc;
cell_list_to_send_by_proc.resize(parallel::size());
for (size_t i = 0; i < parallel::size(); ++i) {
cell_list_to_send_by_proc[i] = convert_to_array(cell_vector_to_send_by_proc[i]);
}
for (size_t i = 0; i < parallel::size(); ++i) {
std::cout << itemName(item_type) << " send to " << i << ": " << cell_list_to_send_by_proc[i] << '\n';
}
} else {
const auto& cell_list_to_send_by_proc = this->_dispatchedInfo<ItemType::cell>().m_list_to_send_by_proc;
const auto& item_is_owned = m_connectivity.template isOwned<item_type>();
using ItemId = ItemIdT<item_type>;
const auto& cell_to_sub_item_matrix = m_connectivity.template getItemToItemMatrix<ItemType::cell, item_type>();
auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
item_list_to_send_by_proc.resize(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
Array<bool> tag(m_connectivity.template numberOf<item_type>());
tag.fill(false);
std::vector<ItemId> item_id_vector;
for (size_t j = 0; j < cell_list_to_send_by_proc[i_rank].size(); ++j) {
const CellId& cell_id = cell_list_to_send_by_proc[i_rank][j];
const auto& cell_sub_item_list = cell_to_sub_item_matrix[cell_id];
for (size_t r = 0; r < cell_sub_item_list.size(); ++r) {
const ItemId& item_id = cell_sub_item_list[r];
if (item_is_owned[item_id] and (not tag[item_id])) {
item_id_vector.push_back(item_id);
tag[item_id] = true;
}
}
}
item_list_to_send_by_proc[i_rank] = convert_to_array(item_id_vector);
}
for (size_t i = 0; i < parallel::size(); ++i) {
std::cout << itemName(item_type) << " send to " << i << ": " << item_list_to_send_by_proc[i] << '\n';
}
}
}
template <size_t Dimension>
template <ItemType item_type>
void
ConnectivityDispatcher<Dimension>::_buildNumberOfItemToExchange()
{
const auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
Array<unsigned int> nb_item_to_send_by_proc(parallel::size());
for (size_t i = 0; i < parallel::size(); ++i) {
nb_item_to_send_by_proc[i] = item_list_to_send_by_proc[i].size();
}
this->_dispatchedInfo<item_type>().m_list_to_send_size_by_proc = nb_item_to_send_by_proc;
std::cout << rang::fgB::yellow << itemName(item_type) << " nb_item_to_send_by_proc = " << nb_item_to_send_by_proc
<< rang::fg::reset << '\n';
this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc = parallel::allToAll(nb_item_to_send_by_proc);
std::cout << rang::fgB::blue << itemName(item_type)
<< " m_list_to_recv_size_by_proc = " << this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc
<< rang::fg::reset << '\n';
}
template <size_t Dimension>
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void
ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
Array<std::remove_const_t<DataType>>& gathered_array)
{
std::cout << rang::style::underline << "-- _gatherFrom[" << itemName(item_type) << "][" << demangle<DataType>()
<< "] --" << rang::style::reset << "\n";
std::vector<Array<const DataType>> recv_item_data_by_proc = this->exchange(data_to_gather);
if constexpr (not(std::is_same_v<std::decay_t<DataType>, CellType>)) {
for (size_t i_rank = 0; i_rank < recv_item_data_by_proc.size(); ++i_rank) {
std::cout << "- from " << rang::fgB::cyan << i_rank << ": recv_" << itemName(item_type)
<< "_data_by_proc=" << recv_item_data_by_proc[i_rank] << rang::fg::reset << '\n';
}
}
const auto& recv_id_correspondance_by_proc = this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
if constexpr (not(std::is_same_v<std::decay_t<DataType>, CellType>)) {
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
std::cout << rang::fgB::green << "* recv_id_correspondance_by_proc[" << i_rank
<< "] = " << recv_id_correspondance_by_proc[i_rank] << rang::fg::reset << '\n';
}
}
Assert(recv_id_correspondance_by_proc.size() == parallel::size());
gathered_array = Array<std::remove_const_t<DataType>>(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
size_t shift = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
Assert(recv_id_correspondance_by_proc[i_rank].size() == recv_item_data_by_proc[i_rank].size());
for (size_t r = 0; r < recv_id_correspondance_by_proc[i_rank].size(); ++r) {
const auto& item_id = recv_id_correspondance_by_proc[i_rank][r];
gathered_array[item_id + shift] = recv_item_data_by_proc[i_rank][r];
}
shift += recv_item_data_by_proc[i_rank].size();
}
if constexpr (not(std::is_same_v<std::decay_t<DataType>, CellType>)) {
std::cout << rang::style::reversed << rang::fgB::yellow << "=> gathered_array = " << gathered_array
<< rang::style::reset << rang::fg::reset << "\n";
}
}
template <size_t Dimension>
template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
void
ConnectivityDispatcher<Dimension>::_gatherFrom(
const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
Array<std::remove_const_t<DataType>>& gathered_array)
{
using MutableDataType = std::remove_const_t<DataType>;
constexpr ItemType item_type = ItemOfItem::item_type;
using ItemId = ItemIdT<item_type>;
const auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
std::vector<Array<MutableDataType>> data_to_send_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
std::vector<MutableDataType> data_by_item_vector;
for (size_t j = 0; j < item_list_to_send_by_proc[i_rank].size(); ++j) {
const ItemId& item_id = item_list_to_send_by_proc[i_rank][j];
const auto& item_data = data_to_gather.itemArray(item_id);
for (size_t l = 0; l < item_data.size(); ++l) {
data_by_item_vector.push_back(item_data[l]);
}
}
data_to_send_by_proc[i_rank] = convert_to_array(data_by_item_vector);
}
const auto& number_of_sub_item_per_item_to_recv_by_proc =
this->_dispatchedInfo<ItemOfItem>().m_number_of_sub_item_per_item_to_recv_by_proc;
std::vector<Array<MutableDataType>> recv_data_to_gather_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
recv_data_to_gather_by_proc[i_rank] =
Array<MutableDataType>(sum(number_of_sub_item_per_item_to_recv_by_proc[i_rank]));
}
parallel::exchange(data_to_send_by_proc, recv_data_to_gather_by_proc);
const size_t recv_array_size = [&] {
size_t size = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
size += recv_data_to_gather_by_proc[i_rank].size();
}
return size;
}();
gathered_array = Array<std::remove_const_t<DataType>>(recv_array_size);
{
size_t l = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
for (size_t j = 0; j < recv_data_to_gather_by_proc[i_rank].size(); ++j) {
gathered_array[l++] = recv_data_to_gather_by_proc[i_rank][j];
}
}
Assert(gathered_array.size() == l);
}
}
template <size_t Dimension>
void
ConnectivityDispatcher<Dimension>::_buildCellNumberIdMap()
{
const auto recv_cell_number_by_proc = this->exchange(m_connectivity.template number<ItemType::cell>());
auto& cell_number_id_map = this->_dispatchedInfo<ItemType::cell>().m_number_to_id_map;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
CellId cell_id = 0;
for (size_t i = 0; i < recv_cell_number_by_proc[i_rank].size(); ++i) {
const int cell_number = recv_cell_number_by_proc[i_rank][i];
auto [iterator, inserted] = cell_number_id_map.insert(std::make_pair(cell_number, cell_id));
if (inserted)
++cell_id;
}
}
}
template <size_t Dimension>
template <typename ItemOfItemT>
void
ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
{
static_assert(ItemOfItemT::item_type == ItemType::cell, "Dispatcher requires to be built using cell as master "
"entities");
const auto& cell_sub_item_number_to_recv_by_proc =
this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
auto& sub_item_number_id_map = this->_dispatchedInfo<ItemOfItemT::sub_item_type>().m_number_to_id_map;
int sub_item_id = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
for (size_t i = 0; i < cell_sub_item_number_to_recv_by_proc[i_rank].size(); ++i) {
int sub_item_number = cell_sub_item_number_to_recv_by_proc[i_rank][i];
std::cout << "looking for " << sub_item_number;
auto [iterator, inserted] = sub_item_number_id_map.insert(std::make_pair(sub_item_number, sub_item_id));
if (inserted) {
std::cout << ": inserted with item_id " << sub_item_id << '\n';
++sub_item_id;
} else {
std::cout << ": *exists* with item_id " << iterator->second << '\n';
}
}
}
}
template <size_t Dimension>
template <typename SubItemOfItemT>
void
ConnectivityDispatcher<Dimension>::_buildNumberOfSubItemPerItemToRecvByProc()
{
const auto& item_to_sub_item_matrix =
m_connectivity.template getItemToItemMatrix<SubItemOfItemT::item_type, SubItemOfItemT::sub_item_type>();
ItemValue<int, SubItemOfItemT::item_type> number_of_sub_item_per_item(m_connectivity);
using ItemId = ItemIdT<SubItemOfItemT::item_type>;
parallel_for(
number_of_sub_item_per_item.numberOfItems(),
PUGS_LAMBDA(const ItemId& j) { number_of_sub_item_per_item[j] = item_to_sub_item_matrix[j].size(); });
this->_dispatchedInfo<SubItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc =
this->exchange(number_of_sub_item_per_item);
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
std::cout << "rcv from " << i_rank << ": " << rang::fgB::red << itemName(SubItemOfItemT::item_type) << "->"
<< itemName(SubItemOfItemT::sub_item_type) << ": "
<< this->_dispatchedInfo<SubItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc[i_rank]
<< rang::fg::reset << '\n';
}
}
template <size_t Dimension>
template <typename SubItemOfItemT>
void
ConnectivityDispatcher<Dimension>::_buildSubItemNumbersToRecvByProc()
{
const std::vector<Array<const int>> sub_item_numbers_to_send_by_proc = [&]() {
const auto& item_to_sub_item_matrix =
m_connectivity.template getItemToItemMatrix<SubItemOfItemT::item_type, SubItemOfItemT::sub_item_type>();
const auto& sub_item_number = m_connectivity.template number<SubItemOfItemT::sub_item_type>();
using ItemId = ItemIdT<SubItemOfItemT::item_type>;
using SubItemId = ItemIdT<SubItemOfItemT::sub_item_type>;
std::vector<Array<const int>> mutable_sub_item_numbers_to_send_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
const auto& item_list_to_send_by_proc = this->_dispatchedInfo<SubItemOfItemT::item_type>().m_list_to_send_by_proc;
std::vector<int> sub_item_numbers_by_item_vector;
for (size_t j = 0; j < item_list_to_send_by_proc[i_rank].size(); ++j) {
const ItemId& item_id = item_list_to_send_by_proc[i_rank][j];
const auto& sub_item_list = item_to_sub_item_matrix[item_id];
for (size_t r = 0; r < sub_item_list.size(); ++r) {
const SubItemId& sub_item_id = sub_item_list[r];
sub_item_numbers_by_item_vector.push_back(sub_item_number[sub_item_id]);
}
}
mutable_sub_item_numbers_to_send_by_proc[i_rank] = convert_to_array(sub_item_numbers_by_item_vector);
}
return mutable_sub_item_numbers_to_send_by_proc;
}();
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
std::cout << "sub_item_nb_to_send " << i_rank << ": " << rang::fgB::magenta << itemName(SubItemOfItemT::item_type)
<< "->" << itemName(SubItemOfItemT::sub_item_type) << ": " << sub_item_numbers_to_send_by_proc[i_rank]
<< rang::fg::reset << '\n';
}
const auto& number_of_sub_item_per_item_to_recv_by_proc =
this->_dispatchedInfo<SubItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc;
std::vector<Array<int>> sub_item_numbers_to_recv_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
sub_item_numbers_to_recv_by_proc[i_rank] = Array<int>(sum(number_of_sub_item_per_item_to_recv_by_proc[i_rank]));
}
parallel::exchange(sub_item_numbers_to_send_by_proc, sub_item_numbers_to_recv_by_proc);
auto& const_sub_item_numbers_to_recv_by_proc =
this->_dispatchedInfo<SubItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
const_sub_item_numbers_to_recv_by_proc.resize(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
const_sub_item_numbers_to_recv_by_proc[i_rank] = sub_item_numbers_to_recv_by_proc[i_rank];
}
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
std::cout << "sub_item_nb_to_recv " << i_rank << ": " << rang::fgB::magenta << itemName(SubItemOfItemT::item_type)
<< "->" << itemName(SubItemOfItemT::sub_item_type) << ": "
<< this->_dispatchedInfo<SubItemOfItemT>().m_sub_item_numbers_to_recv_by_proc[i_rank] << rang::fg::reset
<< '\n';
}
}
template <size_t Dimension>
template <typename ItemOfItemT>
void
ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
{
constexpr ItemType item_type = ItemOfItemT::item_type;
constexpr ItemType sub_item_type = ItemOfItemT::sub_item_type;
constexpr const bool print = (item_type == ItemType::cell) and (sub_item_type == ItemType::node);
const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
const auto& number_of_sub_item_per_item_to_recv_by_proc =
this->_dispatchedInfo<ItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc;
const auto& sub_item_number_id_map = this->_dispatchedInfo<sub_item_type>().m_number_to_id_map;
std::cout << "sub_item_number_id_map = ";
std::cout << rang::fgB::magenta << rang::style::reversed;
for (auto&& i_sub : sub_item_number_id_map) {
std::cout << ' ' << i_sub.first << "->" << i_sub.second;
}
std::cout << rang::fg::reset << rang::style::reset;
std::cout << '\n';
const auto& recv_item_of_item_numbers_by_proc =
this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
std::cout << rang::fgB::blue << "recv_item_of_item_numbers_by_proc[" << i_rank
<< "] = " << recv_item_of_item_numbers_by_proc[i_rank] << rang::fg::reset << '\n';
}
std::vector<std::vector<unsigned int>> item_to_subitem_legacy;
size_t number_of_node_by_cell = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
int l = 0;
for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
std::vector<unsigned int> sub_item_vector;
for (int k = 0; k < number_of_sub_item_per_item_to_recv_by_proc[i_rank][i]; ++k) {
const auto& searched_sub_item_id = sub_item_number_id_map.find(recv_item_of_item_numbers_by_proc[i_rank][l++]);
Assert(searched_sub_item_id != sub_item_number_id_map.end());
sub_item_vector.push_back(searched_sub_item_id->second);
}
number_of_node_by_cell += sub_item_vector.size();
item_to_subitem_legacy.emplace_back(sub_item_vector);
}
}
if constexpr (print) {
for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
std::cout << rang::fgB::magenta << "item_to_subitem_legacy[" << i
<< "] = " << convert_to_array(item_to_subitem_legacy[i]) << rang::fg::reset << '\n';
}
}
Array<unsigned int> item_to_subitem_row_map(item_to_subitem_legacy.size() + 1);
Array<unsigned int> item_to_subitem_list(number_of_node_by_cell);
item_to_subitem_row_map[0] = 0;
for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
item_to_subitem_row_map[i + 1] = item_to_subitem_row_map[i] + item_to_subitem_legacy[i].size();
}
size_t l = 0;
for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
const auto& subitem_list = item_to_subitem_legacy[i];
for (size_t j = 0; j < subitem_list.size(); ++j, ++l) {
item_to_subitem_list[l] = subitem_list[j];
}
}
m_new_descriptor.itemOfItemVector<ItemOfItemT>() = ConnectivityMatrix(item_to_subitem_row_map, item_to_subitem_list);
}
template <size_t Dimension>
template <ItemType item_type>
void
ConnectivityDispatcher<Dimension>::_buildRecvItemIdCorrespondanceByProc()
{
const auto& item_list_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
using ItemId = ItemIdT<item_type>;
std::vector<Array<const ItemId>> recv_item_id_correspondance_by_proc(parallel::size());
const ItemValue<const int, item_type>& item_number = m_connectivity.template number<item_type>();
std::vector<Array<const int>> send_item_number_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
Array<int> send_item_number(item_list_to_send_by_proc[i_rank].size());
const Array<const ItemId> send_item_id = item_list_to_send_by_proc[i_rank];
parallel_for(
send_item_number.size(), PUGS_LAMBDA(size_t j) { send_item_number[j] = item_number[send_item_id[j]]; });
send_item_number_by_proc[i_rank] = send_item_number;
}
const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
std::vector<Array<int>> recv_item_number_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
recv_item_number_by_proc[i_rank] = Array<int>(item_list_to_recv_size_by_proc[i_rank]);
}
parallel::exchange(send_item_number_by_proc, recv_item_number_by_proc);
const auto& item_number_to_id_map = this->_dispatchedInfo<item_type>().m_number_to_id_map;
for (size_t i_rank = 0; i_rank < item_list_to_recv_size_by_proc.size(); ++i_rank) {
Array<ItemId> item_id_correspondance(item_list_to_recv_size_by_proc[i_rank]);
for (size_t l = 0; l < item_list_to_recv_size_by_proc[i_rank]; ++l) {
const int& recv_item_number = recv_item_number_by_proc[i_rank][l];
const auto& searched_item_id = item_number_to_id_map.find(recv_item_number);
Assert(searched_item_id != item_number_to_id_map.end());
item_id_correspondance[l] = searched_item_id->second;
}
recv_item_id_correspondance_by_proc[i_rank] = item_id_correspondance;
}
this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc = recv_item_id_correspondance_by_proc;
}
template <size_t Dimension>
template <ItemType item_type>
void
ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
{
using ItemId = ItemIdT<item_type>;
// Getting references
Array<const size_t> number_of_item_ref_list_per_proc =
parallel::allGather(m_connectivity.template numberOfRefItemList<item_type>());
const size_t number_of_item_list_sender = [&]() {
size_t mutable_number_of_item_list_sender = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
mutable_number_of_item_list_sender += (number_of_item_ref_list_per_proc[i_rank] > 0);
}
return mutable_number_of_item_list_sender;
}();
if (number_of_item_list_sender > 0) {
if (number_of_item_list_sender > 1) {
std::cerr << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red
<< "need to check that knowing procs know the same item_ref_lists!" << rang::fg::reset << '\n';
}
if (number_of_item_list_sender < parallel::size()) {
const size_t sender_rank = [&]() {
size_t i_rank = 0;
for (; i_rank < parallel::size(); ++i_rank) {
if (number_of_item_ref_list_per_proc[i_rank] > 0) {
break;
}
}
return i_rank;
}();
Assert(number_of_item_list_sender < parallel::size());
// sending is boundary property
Array<RefItemListBase::Type> ref_item_list_type{number_of_item_ref_list_per_proc[sender_rank]};
if (parallel::rank() == sender_rank) {
for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
++i_item_ref_list) {
auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
ref_item_list_type[i_item_ref_list] = item_ref_list.type();
}
}
parallel::broadcast(ref_item_list_type, sender_rank);
// sending references tags
Array<RefId::TagNumberType> ref_tag_list{number_of_item_ref_list_per_proc[sender_rank]};
if (parallel::rank() == sender_rank) {
for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
++i_item_ref_list) {
auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
ref_tag_list[i_item_ref_list] = item_ref_list.refId().tagNumber();
}
}
parallel::broadcast(ref_tag_list, sender_rank);
// sending references name size
Array<size_t> ref_name_size_list{number_of_item_ref_list_per_proc[sender_rank]};
if (parallel::rank() == sender_rank) {
for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
++i_item_ref_list) {
auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
ref_name_size_list[i_item_ref_list] = item_ref_list.refId().tagName().size();
}
}
parallel::broadcast(ref_name_size_list, sender_rank);
// sending references name size
Array<RefId::TagNameType::value_type> ref_name_cat{sum(ref_name_size_list)};
if (parallel::rank() == sender_rank) {
size_t i_char = 0;
for (size_t i_item_ref_list = 0; i_item_ref_list < m_connectivity.template numberOfRefItemList<item_type>();
++i_item_ref_list) {
auto item_ref_list = m_connectivity.template refItemList<item_type>(i_item_ref_list);
for (auto c : item_ref_list.refId().tagName()) {
ref_name_cat[i_char++] = c;
}
}
}
parallel::broadcast(ref_name_cat, sender_rank);
std::vector<RefId> ref_id_list = [&]() {
std::vector<RefId> mutable_ref_id_list;
mutable_ref_id_list.reserve(ref_name_size_list.size());
size_t begining = 0;
for (size_t i_ref = 0; i_ref < ref_name_size_list.size(); ++i_ref) {
const size_t size = ref_name_size_list[i_ref];
mutable_ref_id_list.emplace_back(ref_tag_list[i_ref], std::string{&(ref_name_cat[begining]), size});
begining += size;
}
return mutable_ref_id_list;
}();
using block_type = int32_t;
constexpr size_t block_size = sizeof(block_type);
const size_t nb_block = ref_id_list.size() / block_size + (ref_id_list.size() % block_size != 0);
for (size_t i_block = 0; i_block < nb_block; ++i_block) {
ItemValue<block_type, item_type> item_references(m_connectivity);
item_references.fill(0);
if (m_connectivity.template numberOfRefItemList<item_type>() > 0) {
const size_t max_i_ref = std::min(ref_id_list.size(), block_size * (i_block + 1));
for (size_t i_ref = block_size * i_block, i = 0; i_ref < max_i_ref; ++i_ref, ++i) {
block_type ref_bit{1 << i};
auto item_ref_list = m_connectivity.template refItemList<item_type>(i_ref);
const auto& item_list = item_ref_list.list();
for (size_t i_item = 0; i_item < item_list.size(); ++i_item) {
const ItemId& item_id = item_list[i_item];
item_references[item_id] |= ref_bit;
}
}
}
const auto& nb_item_to_send_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_size_by_proc;
const auto& send_item_id_by_proc = this->_dispatchedInfo<item_type>().m_list_to_send_by_proc;
std::vector<Array<const block_type>> send_item_refs_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
Array<block_type> send_item_refs(nb_item_to_send_by_proc[i_rank]);
const Array<const ItemId> send_item_id = send_item_id_by_proc[i_rank];
parallel_for(
send_item_id.size(), PUGS_LAMBDA(size_t l) {
const ItemId& item_id = send_item_id[l];
send_item_refs[l] = item_references[item_id];
});
send_item_refs_by_proc[i_rank] = send_item_refs;
}
std::vector<Array<block_type>> recv_item_refs_by_proc(parallel::size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
recv_item_refs_by_proc[i_rank] =
Array<block_type>(this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc[i_rank]);
}
parallel::exchange(send_item_refs_by_proc, recv_item_refs_by_proc);
const auto& recv_item_id_correspondance_by_proc =
this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
std::vector<block_type> item_refs(m_new_descriptor.template itemNumberVector<item_type>().size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
for (size_t r = 0; r < recv_item_refs_by_proc[i_rank].size(); ++r) {
const ItemId& item_id = recv_item_id_correspondance_by_proc[i_rank][r];
item_refs[item_id] = recv_item_refs_by_proc[i_rank][r];
}
}
const size_t max_i_ref = std::min(ref_id_list.size(), block_size * (i_block + 1));
for (size_t i_ref = block_size * i_block, i = 0; i_ref < max_i_ref; ++i_ref, ++i) {
block_type ref_bit{1 << i};
std::vector<ItemId> item_id_vector;
for (uint32_t i_item = 0; i_item < item_refs.size(); ++i_item) {
const ItemId item_id{i_item};
if (item_refs[item_id] & ref_bit) {
item_id_vector.push_back(item_id);
}
}
Array<const ItemId> item_id_array = convert_to_array(item_id_vector);
RefItemListBase::Type type = ref_item_list_type[i_ref];
m_new_descriptor.addRefItemList(RefItemList<item_type>(ref_id_list[i_ref], item_id_array, type));
}
}
}
}
}
template <size_t Dimension>
void
ConnectivityDispatcher<Dimension>::_dispatchEdges()
{
if constexpr (Dimension > 2) {
this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfCell>();
this->_buildSubItemNumbersToRecvByProc<EdgeOfCell>();
this->_buildSubItemNumberToIdMap<EdgeOfCell>();
this->_buildItemToExchangeLists<ItemType::edge>();
m_new_descriptor.setEdgeNumberVector([&] {
Array<int> edge_number_vector;
this->_gatherFrom(m_connectivity.template number<ItemType::edge>(), edge_number_vector);
return edge_number_vector;
}());
this->_buildItemToSubItemDescriptor<EdgeOfCell>();
this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfEdge>();
this->_buildSubItemNumbersToRecvByProc<NodeOfEdge>();
this->_buildItemToSubItemDescriptor<NodeOfEdge>();
this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfFace>();
this->_buildSubItemNumbersToRecvByProc<EdgeOfFace>();
this->_buildItemToSubItemDescriptor<EdgeOfFace>();
m_new_descriptor.setFaceEdgeIsReversed([&] {
Array<bool> face_edge_is_reversed;
this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), face_edge_is_reversed);
return face_edge_is_reversed;
}());
m_new_descriptor.setEdgeOwnerVector([&] {
Array<int> edge_owner_vector;
this->_gatherFrom(this->_dispatchedInfo<ItemType::edge>().m_new_owner, edge_owner_vector);
return edge_owner_vector;
}());
this->_buildItemReferenceList<ItemType::edge>();
}
}
template <size_t Dimension>
void
ConnectivityDispatcher<Dimension>::_dispatchFaces()
{
if constexpr (Dimension > 1) {
this->_buildNumberOfSubItemPerItemToRecvByProc<FaceOfCell>();
this->_buildSubItemNumbersToRecvByProc<FaceOfCell>();
this->_buildSubItemNumberToIdMap<FaceOfCell>();
this->_buildItemToExchangeLists<ItemType::face>();
this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfFace>();
this->_buildSubItemNumbersToRecvByProc<NodeOfFace>();
this->_buildItemToSubItemDescriptor<NodeOfFace>();
m_new_descriptor.setFaceNumberVector([&] {
Array<int> face_number_vector;
this->_gatherFrom(m_connectivity.template number<ItemType::face>(), face_number_vector);
return face_number_vector;
}());
this->_buildItemToSubItemDescriptor<FaceOfCell>();
m_new_descriptor.setCellFaceIsReversed([&] {
Array<bool> cell_face_is_reversed;
this->_gatherFrom(m_connectivity.cellFaceIsReversed(), cell_face_is_reversed);
return cell_face_is_reversed;
}());
m_new_descriptor.setFaceOwnerVector([&] {
Array<int> face_owner_vector;
this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, face_owner_vector);
return face_owner_vector;
}());
this->_buildItemReferenceList<ItemType::face>();
}
}
template <size_t Dimension>
ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType& connectivity)
: m_connectivity(connectivity)
{
{
Array connectivity_id_list = parallel::allGather(connectivity.id());
if (min(connectivity_id_list) != max(connectivity_id_list)) {
throw UnexpectedError("connectivity ids diverged in parallel");
}
}
std::cout << "--- initial mesh ---\n";
std::cout << "--- cells\n";
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
std::cout << "cell " << cell_id << "(" << connectivity.cellNumber()[cell_id] << "):" << rang::style::bold
<< ((connectivity.cellIsOwned()[cell_id]) ? 'O' : 'G') << rang::style::reset << ": ";
auto node_list = connectivity.cellToNodeMatrix()[cell_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
std::cout << ' ' << node_id << "(" << connectivity.nodeNumber()[node_id] << "):" << rang::style::bold
<< ((connectivity.nodeIsOwned()[node_id]) ? 'O' : 'G') << rang::style::reset;
}
std::cout << '\n';
}
if constexpr (Dimension > 1) {
std::cout << "--- faces\n";
for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
std::cout << "cell " << face_id << "(" << connectivity.faceNumber()[face_id] << "):" << rang::style::bold
<< ((connectivity.faceIsOwned()[face_id]) ? 'O' : 'G') << rang::style::reset << ": ";
auto node_list = connectivity.faceToNodeMatrix()[face_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
std::cout << ' ' << node_id << "(" << connectivity.nodeNumber()[node_id] << "):" << rang::style::bold
<< ((connectivity.nodeIsOwned()[node_id]) ? 'O' : 'G') << rang::style::reset;
}
std::cout << '\n';
}
}
std::cout << "--- nodes\n";
for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
std::cout << "node " << node_id << "(" << connectivity.nodeNumber()[node_id] << "):" << rang::style::bold
<< ((connectivity.nodeIsOwned()[node_id]) ? 'O' : 'G') << rang::style::reset << ": ";
std::cout << '\n';
}
this->_buildNewOwner<ItemType::cell>();
if constexpr (Dimension > 1) {
this->_buildNewOwner<ItemType::face>();
}
if constexpr (Dimension > 2) {
this->_buildNewOwner<ItemType::edge>();
}
this->_buildNewOwner<ItemType::node>();
std::cout << rang::style::reversed << "--- cell to exchange ---" << rang::style::reset << '\n';
this->_buildItemToExchangeLists<ItemType::cell>();
std::cout << rang::style::reversed << "--- number of node per cell ---" << rang::style::reset << '\n';
this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfCell>();
std::cout << rang::style::reversed << "--- node by cell to recv by proc ---" << rang::style::reset << '\n';
this->_buildSubItemNumbersToRecvByProc<NodeOfCell>();
m_new_descriptor.setCellNumberVector([&] {
Array<int> cell_number_vector;
this->_gatherFrom(m_connectivity.template number<ItemType::cell>(), cell_number_vector);
return cell_number_vector;
}());
std::cout << rang::fgB::green << "cell_number_vector = " << m_new_descriptor.cellNumberVector() << rang::fg::reset
<< "\n";
std::cout << rang::style::reversed << "--- node of cell to id map ---" << rang::style::reset << '\n';
this->_buildSubItemNumberToIdMap<NodeOfCell>();
std::cout << rang::style::reversed << "--- node to exchange ---" << rang::style::reset << '\n';
this->_buildItemToExchangeLists<ItemType::node>();
m_new_descriptor.setCellTypeVector([&] {
Array<CellType> cell_type_vector;
this->_gatherFrom(m_connectivity.cellType(), cell_type_vector);
return cell_type_vector;
}());
m_new_descriptor.setCellOwnerVector([&] {
Array<int> cell_owner_vector;
this->_gatherFrom(this->_dispatchedInfo<ItemType::cell>().m_new_owner, cell_owner_vector);
return cell_owner_vector;
}());
m_new_descriptor.setNodeNumberVector([&] {
Array<int> node_number_vector;
this->_gatherFrom(m_connectivity.template number<ItemType::node>(), node_number_vector);
return node_number_vector;
}());
m_new_descriptor.setNodeOwnerVector([&] {
Array<int> node_owner_vector;
this->_gatherFrom(this->_dispatchedInfo<ItemType::node>().m_new_owner, node_owner_vector);
return node_owner_vector;
}());
this->_buildItemToSubItemDescriptor<NodeOfCell>();
this->_buildItemReferenceList<ItemType::cell>();
this->_dispatchFaces();
this->_dispatchEdges();
this->_buildItemReferenceList<ItemType::node>();
std::cout << "m_new_descriptor.nodeNumberVector() = " << m_new_descriptor.nodeNumberVector() << '\n';
std::cout << rang::fgB::blue << m_new_descriptor.nodeNumberVector() << rang::fg::reset << '\n';
std::cout << "m_new_descriptor.nodeOwnerVector() = " << m_new_descriptor.nodeOwnerVector() << '\n';
std::cout << rang::fgB::blue << m_new_descriptor.nodeOwnerVector() << rang::fg::reset << '\n';
std::cout << "m_new_descriptor.cellOwnerVector().size() = " << m_new_descriptor.cellOwnerVector().size()
<< '\n';
std::cout << rang::fgB::blue << m_new_descriptor.cellOwnerVector() << rang::fg::reset << '\n';
std::cout << "m_new_descriptor.cellTypeVector().size() = " << m_new_descriptor.cellTypeVector().size()
<< '\n';
std::cout << "m_new_descriptor.cellToNodeMatrix().numberOfRows() = "
<< m_new_descriptor.cellToNodeMatrix().numberOfRows() << '\n';
for (CellId cell_id = 0; cell_id < m_new_descriptor.cellToNodeMatrix().numberOfRows(); ++cell_id) {
std::cout << "- cell [" << cell_id << "](" // << m_new_descriptor.cellNumberVector()[cell_id]
<< ")" << m_new_descriptor.cellToNodeMatrix()[cell_id] << '\n';
}
m_dispatched_connectivity = ConnectivityType::build(m_new_descriptor);
{
Array connectivity_id_list = parallel::allGather(m_dispatched_connectivity->id());
if (min(connectivity_id_list) != max(connectivity_id_list)) {
throw UnexpectedError("connectivity ids diverged in parallel");
}
}
}
template ConnectivityDispatcher<1>::ConnectivityDispatcher(const ConnectivityType&);
template ConnectivityDispatcher<2>::ConnectivityDispatcher(const ConnectivityType&);
template ConnectivityDispatcher<3>::ConnectivityDispatcher(const ConnectivityType&);