Commit aefbe074 authored by wenqing's avatar wenqing
Browse files

Merge branch 'parallel_comp_T_H_element_LF' into 'master'

Parallel computing with Taylor-Hood elements: implementation and tests with HydroMechanics process

See merge request !3866
parents 0ad843c8 78e0d481
Pipeline #7128 failed with stages
in 68 minutes and 35 seconds
......@@ -15,13 +15,12 @@
#include "NodePartitionedMeshReader.h"
#include "BaseLib/Logging.h"
#ifdef USE_PETSC
#include <mpi.h>
#endif
#include <numeric>
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/RunTime.h"
#include "MeshLib/Elements/Elements.h"
#include "MeshLib/MeshEnums.h"
......@@ -381,10 +380,55 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::newMesh(
std::vector<MeshLib::Element*> const& mesh_elems,
MeshLib::Properties const& properties) const
{
std::vector<std::size_t> gathered_n_active_base_nodes(_mpi_comm_size);
MPI_Allgather(&_mesh_info.active_base_nodes,
1,
MPI_UNSIGNED_LONG,
gathered_n_active_base_nodes.data(),
1,
MPI_UNSIGNED_LONG,
_mpi_comm);
std::vector<std::size_t> gathered_n_base_node_before_rank;
gathered_n_base_node_before_rank.reserve(
gathered_n_active_base_nodes.size());
for (std::size_t i = 0; i < gathered_n_active_base_nodes.size(); i++)
{
gathered_n_base_node_before_rank.push_back(
std::transform_reduce(begin(gathered_n_active_base_nodes),
begin(gathered_n_active_base_nodes) + i, 0.0,
std::plus{}, [](auto const n) { return n; }));
}
std::vector<std::size_t> gathered_n_active_high_order_nodes(_mpi_comm_size);
std::size_t const n_active_high_order_nodes =
_mesh_info.active_nodes - _mesh_info.active_base_nodes;
MPI_Allgather(&n_active_high_order_nodes,
1,
MPI_UNSIGNED_LONG,
gathered_n_active_high_order_nodes.data(),
1,
MPI_UNSIGNED_LONG,
_mpi_comm);
std::vector<std::size_t> gathered_n_high_order_node_before_rank;
gathered_n_high_order_node_before_rank.reserve(
gathered_n_active_high_order_nodes.size());
for (std::size_t i = 0; i < gathered_n_active_high_order_nodes.size(); i++)
{
gathered_n_high_order_node_before_rank.push_back(std::transform_reduce(
begin(gathered_n_active_high_order_nodes),
begin(gathered_n_active_high_order_nodes) + i, 0.0, std::plus{},
[](auto const n) { return n; }));
}
return new MeshLib::NodePartitionedMesh(
mesh_name, mesh_nodes, glb_node_ids, mesh_elems, properties,
_mesh_info.global_nodes, _mesh_info.active_base_nodes,
_mesh_info.active_nodes);
_mesh_info.global_base_nodes, _mesh_info.global_nodes,
_mesh_info.active_base_nodes, _mesh_info.active_nodes,
std::move(gathered_n_base_node_before_rank),
std::move(gathered_n_high_order_node_before_rank));
}
void NodePartitionedMeshReader::setNodes(
......
......@@ -28,8 +28,13 @@ public:
/// Construct a mesh subset from vector of nodes on the given mesh.
/// \param msh Mesh
/// \param vec_items Vector of Node pointers.
MeshSubset(const Mesh& msh, std::vector<Node*> const& vec_items)
: _msh(msh), _nodes(vec_items)
/// \param use_TaylorHood_elements Flag to indicate whether the Taylor-Hood
/// elements are used.
MeshSubset(const Mesh& msh, std::vector<Node*> const& vec_items,
const bool use_TaylorHood_elements = false)
: _msh(msh),
_nodes(vec_items),
_use_TaylorHood_elements(use_TaylorHood_elements)
{
// If the mesh nodes and the given nodes point to the same vector, they
// must be equal.
......@@ -82,6 +87,8 @@ public:
return _nodes[i]->getID();
}
bool useTaylorHoodElements() const { return _use_TaylorHood_elements; }
std::vector<Element*>::const_iterator elementsBegin() const
{
return _msh.getElements().cbegin();
......@@ -99,5 +106,6 @@ public:
private:
Mesh const& _msh;
std::vector<Node*> const& _nodes;
bool const _use_TaylorHood_elements;
};
} // namespace MeshLib
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
* Created on October 1, 2021, 2:13 PM
*/
#include "NodePartitionedMesh.h"
namespace MeshLib
{
std::vector<int> getEndNodeIDRanks(
std::size_t const n_global_nodes,
std::vector<std::size_t> const& gathered_n_active_base_node_at_rank,
std::vector<std::size_t> const& gathered_n_active_high_order_node_at_rank)
{
std::vector<int> data;
std::transform(gathered_n_active_base_node_at_rank.begin() + 1,
gathered_n_active_base_node_at_rank.end(),
gathered_n_active_high_order_node_at_rank.begin() + 1,
std::back_inserter(data), std::plus<int>());
data.push_back(n_global_nodes);
return data;
}
NodePartitionedMesh::NodePartitionedMesh(
const std::string& name,
const std::vector<Node*>& nodes,
const std::vector<std::size_t>& glb_node_ids,
const std::vector<Element*>& elements,
Properties properties,
const std::size_t n_global_base_nodes,
const std::size_t n_global_nodes,
const std::size_t n_active_base_nodes,
const std::size_t n_active_nodes,
std::vector<std::size_t>&& gathered_n_base_node_at_rank,
std::vector<std::size_t>&& gathered_n_high_order_node_at_rank)
: Mesh(name, nodes, elements, properties),
_global_node_ids(glb_node_ids),
_n_global_base_nodes(n_global_base_nodes),
_n_global_nodes(n_global_nodes),
_n_active_base_nodes(n_active_base_nodes),
_n_active_nodes(n_active_nodes),
_gathered_n_active_base_node_at_rank(
std::move(gathered_n_base_node_at_rank)),
_gathered_n_active_high_order_node_at_rank(
std::move(gathered_n_high_order_node_at_rank)),
_gathered_end_node_ids_ranks(getEndNodeIDRanks(
n_global_nodes, _gathered_n_active_base_node_at_rank,
_gathered_n_active_high_order_node_at_rank)),
_is_single_thread(false)
{
}
bool NodePartitionedMesh::isGhostNode(const std::size_t node_id) const
{
return node_id >= _n_active_nodes;
}
std::size_t NodePartitionedMesh::getMaximumNConnectedNodesToNode() const
{
auto const& nodes_connections =
MeshLib::calculateNodesConnectedByElements(*this);
auto const max_connections = std::max_element(
nodes_connections.cbegin(), nodes_connections.cend(),
[](auto const& connections_node_a, auto const& connections_node_b)
{ return (connections_node_a.size() < connections_node_b.size()); });
// Return the number of connected nodes +1 for the node itself.
return max_connections->size() + 1;
}
int NodePartitionedMesh::getPartitionID(const int global_node_id) const
{
return std::upper_bound(std::cbegin(_gathered_end_node_ids_ranks),
std::cend(_gathered_end_node_ids_ranks),
global_node_id) -
_gathered_end_node_ids_ranks.begin();
}
} // namespace MeshLib
......@@ -34,6 +34,7 @@ public:
explicit NodePartitionedMesh(const Mesh& mesh)
: Mesh(mesh),
_global_node_ids(mesh.getNumberOfNodes()),
_n_global_base_nodes(mesh.getNumberOfBaseNodes()),
_n_global_nodes(mesh.getNumberOfNodes()),
_n_active_base_nodes(mesh.getNumberOfBaseNodes()),
_n_active_nodes(mesh.getNumberOfNodes()),
......@@ -56,25 +57,32 @@ public:
\param elements Vector for elements. Ghost elements are stored
after regular (non-ghost) elements.
\param properties Mesh property.
\param n_global_base_nodes Number of the base nodes of the global mesh.
\param n_global_nodes Number of all nodes of the global mesh.
\param n_active_base_nodes Number of the active base nodes.
\param n_active_nodes Number of all active nodes.
\param gathered_n_base_node_at_rank Number of base nodes up to each
rank.
\param gathered_n_high_order_node_at_rank Number of high order nodes up
to each rank.
*/
NodePartitionedMesh(const std::string& name,
const std::vector<Node*>& nodes,
const std::vector<std::size_t>& glb_node_ids,
const std::vector<Element*>& elements,
Properties properties,
const std::size_t n_global_nodes,
const std::size_t n_active_base_nodes,
const std::size_t n_active_nodes)
: Mesh(name, nodes, elements, properties),
_global_node_ids(glb_node_ids),
_n_global_nodes(n_global_nodes),
_n_active_base_nodes(n_active_base_nodes),
_n_active_nodes(n_active_nodes),
_is_single_thread(false)
NodePartitionedMesh(
const std::string& name,
const std::vector<Node*>& nodes,
const std::vector<std::size_t>& glb_node_ids,
const std::vector<Element*>& elements,
Properties properties,
const std::size_t n_global_base_nodes,
const std::size_t n_global_nodes,
const std::size_t n_active_base_nodes,
const std::size_t n_active_nodes,
std::vector<std::size_t>&& gathered_n_base_node_at_rank,
std::vector<std::size_t>&& gathered_n_high_order_node_at_rank);
/// Get the number of nodes of the global mesh for linear elements.
std::size_t getNumberOfGlobalBaseNodes() const
{
return _n_global_base_nodes;
}
/// Get the number of all nodes of the global mesh.
......@@ -85,29 +93,10 @@ public:
return _global_node_ids[node_id];
}
/// Get the number of the active nodes of the partition for linear elements.
std::size_t getNumberOfActiveBaseNodes() const
{
return _n_active_base_nodes;
}
/// Get the number of all active nodes of the partition.
std::size_t getNumberOfActiveNodes() const { return _n_active_nodes; }
/// Check whether a node with ID of node_id is a ghost node
bool isGhostNode(const std::size_t node_id) const
{
if (node_id < _n_active_base_nodes)
{
return false;
}
if (!isBaseNode(*_nodes[node_id],
getElementsConnectedToNode(node_id)) &&
node_id < getLargestActiveNodeID())
{
return false;
}
return true;
}
bool isGhostNode(const std::size_t node_id) const;
/// Get the largest ID of active nodes for higher order elements in a
/// partition.
......@@ -116,19 +105,27 @@ public:
return getNumberOfBaseNodes() + _n_active_nodes - _n_active_base_nodes;
}
std::size_t getGatheredNumberOfNonGhostBaseNodesAtRank(
int const partition_id) const
{
return _gathered_n_active_base_node_at_rank[partition_id];
}
std::size_t getGatheredNumberOfNonGhostHighOrderNodesAtRank(
int const partition_id) const
{
return _gathered_n_active_high_order_node_at_rank[partition_id];
}
// TODO I guess that is a simplified version of computeSparsityPattern()
/// Get the maximum number of connected nodes to node.
std::size_t getMaximumNConnectedNodesToNode() const
std::size_t getMaximumNConnectedNodesToNode() const;
int getPartitionID(const int global_node_id) const;
int getNumberOfPartitions() const
{
auto const& nodes_connections =
MeshLib::calculateNodesConnectedByElements(*this);
auto const max_connections = std::max_element(
nodes_connections.cbegin(), nodes_connections.cend(),
[](auto const& connections_node_a, auto const& connections_node_b) {
return (connections_node_a.size() < connections_node_b.size());
});
// Return the number of connected nodes +1 for the node itself.
return max_connections->size() + 1;
return _gathered_n_active_base_node_at_rank.size();
}
bool isForSingleThread() const { return _is_single_thread; }
......@@ -137,6 +134,9 @@ private:
/// Global IDs of nodes of a partition
std::vector<std::size_t> _global_node_ids;
/// Number of the nodes of the global mesh linear interpolations.
std::size_t _n_global_base_nodes;
/// Number of all nodes of the global mesh.
std::size_t _n_global_nodes;
......@@ -146,6 +146,19 @@ private:
/// Number of the all active nodes.
std::size_t _n_active_nodes;
/// Gathered the numbers of non-ghost base nodes before each rank. Its size
/// is the
/// rank size.
std::vector<std::size_t> _gathered_n_active_base_node_at_rank;
/// Gathered the numbers of non-ghost high order nodes before each rank.
/// Its size is
/// the rank size.
std::vector<std::size_t> _gathered_n_active_high_order_node_at_rank;
/// Gathered the end node id of each rank.
std::vector<int> _gathered_end_node_ids_ranks;
const bool _is_single_thread;
};
......
......@@ -14,6 +14,7 @@
#include "BaseLib/Error.h"
#include "MeshLib/MeshSubset.h"
#include "MeshLib/Node.h" // for bool isBaseNode(Node const& node);
#ifdef USE_PETSC
#include "MeshLib/NodePartitionedMesh.h"
......@@ -58,7 +59,7 @@ MeshComponentMap MeshComponentMap::getSubset(
std::vector<int> const& new_global_component_ids) const
{
{ // Testing first an assumption met later in the code that the meshes for
// all the bulk_mesh_subsets are equal.
// all the bulk_mesh_subsets are equal.
auto const first_mismatch =
std::adjacent_find(begin(bulk_mesh_subsets), end(bulk_mesh_subsets),
[](auto const& a, auto const& b)
......@@ -322,6 +323,61 @@ void MeshComponentMap::createSerialMeshComponentMap(
}
#ifdef USE_PETSC
GlobalIndexType getGlobalIndexWithTaylorHoodElement(
MeshLib::NodePartitionedMesh const& partitioned_mesh,
std::size_t const global_node_id,
int const number_of_components_at_base_node,
int const number_of_components_at_high_order_node,
int const component_id_at_base_node,
int const component_id_at_high_order_node,
bool const is_base_node)
{
int const partition_id = partitioned_mesh.getPartitionID(global_node_id);
auto const n_total_active_base_nodes_before_this_rank =
partitioned_mesh.getGatheredNumberOfNonGhostBaseNodesAtRank(
partition_id);
auto const n_total_active_high_order_nodes_before_this_rank =
partitioned_mesh.getGatheredNumberOfNonGhostHighOrderNodesAtRank(
partition_id);
auto const node_id_offset =
n_total_active_base_nodes_before_this_rank +
n_total_active_high_order_nodes_before_this_rank;
auto const index_offset = n_total_active_base_nodes_before_this_rank *
number_of_components_at_base_node +
n_total_active_high_order_nodes_before_this_rank *
number_of_components_at_high_order_node;
if (is_base_node)
{
return static_cast<GlobalIndexType>(
index_offset + (global_node_id - node_id_offset) *
number_of_components_at_base_node) +
component_id_at_base_node;
}
int const n_active_base_nodes_of_this_partition =
(partition_id == partitioned_mesh.getNumberOfPartitions() - 1)
? partitioned_mesh.getNumberOfGlobalBaseNodes() -
n_total_active_base_nodes_before_this_rank
: partitioned_mesh.getGatheredNumberOfNonGhostBaseNodesAtRank(
partition_id + 1) -
n_total_active_base_nodes_before_this_rank;
return static_cast<GlobalIndexType>(
index_offset +
n_active_base_nodes_of_this_partition *
number_of_components_at_base_node +
(global_node_id - node_id_offset -
n_active_base_nodes_of_this_partition) *
number_of_components_at_high_order_node) +
component_id_at_high_order_node;
}
void MeshComponentMap::createParallelMeshComponentMap(
std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
{
......@@ -335,17 +391,36 @@ void MeshComponentMap::createParallelMeshComponentMap(
"the order type of ComponentOrder::BY_LOCATION");
}
// get number of unknowns
const MeshLib::NodePartitionedMesh& partitioned_mesh =
static_cast<const MeshLib::NodePartitionedMesh&>(
components[0].getMesh());
//
// get the number of unknowns and the number of components at extra node
//
GlobalIndexType num_unknowns = 0;
int components_at_high_order_nodes = 0;
for (auto const& c : components)
{
const MeshLib::NodePartitionedMesh& partitioned_mesh =
static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
num_unknowns += partitioned_mesh.getNumberOfGlobalNodes();
if (partitioned_mesh.getNumberOfNodes() == c.getNumberOfNodes())
{
components_at_high_order_nodes++;
num_unknowns += partitioned_mesh.getNumberOfGlobalNodes();
}
else
{
num_unknowns += partitioned_mesh.getNumberOfGlobalBaseNodes();
}
}
// construct dict (and here we number global_index by component type)
//
// If !has_TaylorHood_element, components_at_base_nodes are the number of
// components at any nodes, either base nodes or high order nodes.
int const components_at_base_nodes = components.size();
int comp_id = 0;
int comp_id_at_high_order_node = 0;
_num_global_dof = 0;
_num_local_dof = 0;
for (auto const& c : components)
......@@ -355,15 +430,32 @@ void MeshComponentMap::createParallelMeshComponentMap(
std::size_t const mesh_id = c.getMeshID();
const MeshLib::NodePartitionedMesh& partitioned_mesh =
static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
const auto& sub_mesh_nodes = c.getNodes();
bool const use_whole_nodes =
(partitioned_mesh.getNumberOfNodes() == c.getNumberOfNodes());
// mesh items are ordered first by node, cell, ....
for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
{
GlobalIndexType global_index = static_cast<GlobalIndexType>(
components.size() * partitioned_mesh.getGlobalNodeID(j) +
comp_id);
const bool is_ghost = partitioned_mesh.isGhostNode(
partitioned_mesh.getNode(j)->getID());
const auto node_id = c.getNodeID(j);
const bool is_base_node =
MeshLib::isBaseNode(*sub_mesh_nodes[j],
partitioned_mesh.getElementsConnectedToNode(
*sub_mesh_nodes[j]));
const auto global_node_id =
partitioned_mesh.getGlobalNodeID(node_id);
GlobalIndexType global_index =
(!c.useTaylorHoodElements())
? static_cast<GlobalIndexType>(
components_at_base_nodes * global_node_id + comp_id)
: getGlobalIndexWithTaylorHoodElement(
partitioned_mesh, global_node_id,
components_at_base_nodes,
components_at_high_order_nodes, comp_id,
comp_id_at_high_order_node, is_base_node);
const bool is_ghost = partitioned_mesh.isGhostNode(node_id);
if (is_ghost)
{
_ghosts_indices.push_back(global_index);
......@@ -380,11 +472,21 @@ void MeshComponentMap::createParallelMeshComponentMap(
_num_local_dof++;
}
_dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Node, j),
comp_id, global_index));
_dict.insert(
Line(Location(mesh_id, MeshLib::MeshItemType::Node, node_id),
comp_id, global_index));
}
if (use_whole_nodes)
{
_num_global_dof += partitioned_mesh.getNumberOfGlobalNodes();
comp_id_at_high_order_node++;
}
else
{
_num_global_dof += partitioned_mesh.getNumberOfGlobalBaseNodes();
}
_num_global_dof += partitioned_mesh.getNumberOfGlobalNodes();
comp_id++;
}
}
......
......@@ -180,6 +180,12 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
&mesh);
const bool use_TaylorHood_elements =
variable_p->getShapeFunctionOrder() !=
variable_u->getShapeFunctionOrder()
? true
: false;
HydroMechanicsProcessData<DisplacementDim> process_data{
materialIDs(mesh),
std::move(media_map),
......@@ -188,7 +194,8 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
specific_body_force,
mass_lumping,
hydraulic_process_id,
mechanics_related_process_id};
mechanics_related_process_id,
use_TaylorHood_elements};
SecondaryVariableCollection secondary_variables;
......
......@@ -90,12 +90,14 @@ template <int DisplacementDim>
void HydroMechanicsProcess<DisplacementDim>::constructDofTable()
{
// Create single component dof in every of the mesh's nodes.
_mesh_subset_all_nodes =
std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
_mesh_subset_all_nodes = std::make_unique<MeshLib::MeshSubset>(
_mesh, _mesh.getNodes(), _process_data.use_TaylorHood_elements);
// Create single component dof in the mesh's base nodes.
_base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
_mesh_subset_base_nodes =
std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
_mesh_subset_base_nodes = std::make_unique<MeshLib::MeshSubset>(
_mesh, _base_nodes, _process_data.use_TaylorHood_elements);
// TODO move the two data members somewhere else.
// for extrapolation of secondary variables of stress or strain
......
......@@ -64,6 +64,8 @@ struct HydroMechanicsProcessData
/// ID of the processes that contains mechanical process.
int const mechanics_related_process_id;
const bool use_TaylorHood_elements;
MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
std::array<MeshLib::PropertyVector<double>*, 3> principal_stress_vector = {
nullptr, nullptr, nullptr};
......
......@@ -586,3 +586,91 @@ AddTest(
flow_gravity_ts_16_t_40000000.000000.vtu flow_gravity_ts_16_t_40000000.000000.vtu pressure pressure 1e-9 0
flow_gravity_ts_16_t_40000000.000000.vtu flow_gravity_ts_16_t_40000000.000000.vtu velocity velocity 1e-9 0
)
#Parallel computing with PETSc enabled.
AddTest(
NAME ParallelFEM_SimpleHM_SingleMesh
PATH HydroMechanics/ParallelComputing/SimpleHM/SingleMesh
EXECUTABLE ogs
EXECUTABLE_ARGS drainage.prj
WRAPPER mpirun
WRAPPER_ARGS -np 2
TESTER vtkdiff
REQUIREMENTS OGS_USE_MPI
DIFF_DATA
drainage_ts_10_t_10_000000_0.vtu drainage_ts_10_t_10_000000_0.vtu pressure pressure 1e-10 1e-9
drainage_ts_10_t_10_000000_1.vtu drainage_ts_10_t_10_000000_1.vtu pressure pressure 1e-10 1e-9
drainage_ts_10_t_10_000000_0.vtu drainage_ts_10_t_10_000000_0.vtu q q 1e-10 1e-9
drainage_ts_10_t_10_000000_1.vtu drainage_ts_10_t_10_000000_1.vtu q q 1e-10 1e-9
drainage_ts_10_t_10_000000_0.vtu drainage_ts_10_t_10_000000_0.vtu displacement displacement 1e-10 1e-9
drainage_ts_10_t_10_000000_1.vtu drainage_ts_10_t_10_000000_1.vtu displacement displacement 1e-10 1e-9