Parallelism in b2000++

b2000++ can exploit shared memory parallelism as well as distributed memory. For the distributed parallelism `MPI`_ is utilized, and there are two different approaches:

• A centralized approach using workers to solve the linear equation systems, while the root process acts as the main process that has to hold the complete problem and does everything outside the linear equation system solving on its own.

• A partitioned mode, where all processes operate on a part of the overall domain and, thus, also the computations outside the linear equation solving can make use of the distributed memory parallelism.

Currently only MUMPS is supported as MPI parallel solver.

For shared memory parallelism, b2000++ uses `Intel's Threading Building Blocks (oneTBB)`_. But the deployed linear eation solvers like `MUMPS`_ may also make use of `OpenMP`_.

See also the section on parallelism in the user manual.

MPI parallelization

The default parallelization strategy is the centralized approach with the root process following the usual code path and doing all the work outside the linear equation solving for the complete domain. This only allows for a fairly limited utilization of distributed parallelism, as all the parts outside the linear equation solver do not benefit from more MPI processes and the complete problem itself has to fit into the memory of the root process. All other processes follow a different code path and essentially just stand-by waiting on commands from the root process.

The partitioned mode has to be activated by passing the -partitioned option to the b2000++ execution on the command line. In this mode the domain is split across processes, with each process only working on a subset of the elements in the domain.

To compile b2000++ with MPI support, provide -DUSE_MPI=ON to the cmake configuration step.

The interaction with MPI is encapsulated in a communicator class, which allows for replacing the actual communication infrastructure beneath, if needed. Note however, that the interface is basically mirroring the MPI interface.

MPI Encapsulation

The communication layer in b2000++ is provided by the abstract b2000::Communicator class in src/utils/b2comm.H. It encapsulates MPI if available and otherwise provides a serial implementation of the respective functions to allow for the code to operate the same way in either case (if possible). The derived class for the MPI implementation is provided in src/utils/b2comm_MPI.H as b2000::B2mpi.

In its main routine b2000++ creates a pointer of the abstract class and constructs the actual communicator object with b2000::CreateComm, which will yield a b2000::B2mpi object, if MPI is supported. This constructor will invoke MPI_Init and hence has to be invoked on all processes, but only once. Alternative constructors exist. For example, you can also create a B2mpi communicator object out of a provided MPI_Comm communicator. Note however, that a global reference on the communicator to use is kept in the SolverSlave class in form of a weak_ptr. This reference is established in the main program directly after creating the communicator.

The implemented interface is not complete and only provides the functions currently needed by b2000++.

Main-Worker approach

In this approach only the main process (rank 0) runs the actual main program of b2000++. All other processes will enter the dedicated command receiving routine implemented in the b2000::linalg::SolverSlave class. The main process of the communicator does not hold an instance of this slave, but instead relies on static properties of the class itself. All other processes act as workers and instantiate an object of this class and invoke its Run() method.

SolverSlave: Run

The Run() method of the SolverSlave represents, after an initialization with the debug level to use, an endless loop that listens on commands from the master. These commands are given in the form of broadcasts of a small buffer of uint_64 that encodes the instruction to perform on all Slaves. The first entry in this buffer array indicates which instruction to perform and may be:

• 0 to terminate the execution,

• 1 to perform a MUMPS operation,

The main process

Rank 0 is used as main process, holds the global matrix and performs the actual computation, driving the execution of the simulation. It assembles the global matrix in the CSC format (compressed sparse columns).

MUMPS worker

b2000++ may use the linear algebra solution for matrices of different structure. The corresponding MUMPS_STRUC data structures are kept separately on the worker processes. To allow this to be equally handled on both the main and the worker processes, a MumpsWorkers class is employed (implemented in src/utils/linear_algebra/b2mumps_workers.H). There are two maps utilized in the MumpsWorkers class to allow for different solver invocations. One for real-numbered matrices (dmumps_object_map) and one for complex numbered matrices (zmumps_object_map). The key in the map is obtained by incrementing a counter upon worker creation on the main process.

Partitioned approach

In the partitioned approach (activated via the command-line flag -partitioned), all processes follow the main code path, but only take care of a subset of the domains elements. model.partitioned provides a flag to indicate whether we are operating in partitioned mode or not.

When we partition the mesh, we need to consider elements in branches and processes with global element ids.

Branches may be removed eventually, as it is not really used anymore. For now we distribute elements in each branch across all processes. For each branch we store the offset in the global list of elements as partition_first. In addition, each element maintains a local internal id, which is a consecutive numbering of all elements across all branches. Each element also maintains a global internal id, which indicates the consecutive numbering across branches, but over all processes, hence, it is not consecutive across branches on individual processes, when there is more than one process.

See the implementation in io/b2000_db/b2000_db_v3/b2domain_database.C and its header. The Branch struct there in the Domain class now has a split_partition method to distribute the elements within each branch across processes. This is a simple splitting along the list of elements, without considering adjacency and communication needs.

We also need to take care of the dofs, their count and numbering. Each process needs to have all dofs for each of its elements, but on its surfaces those dofs overlap with the neighbors, requiring a scheme to identify global ids for each dof.

Where all dofs are used, we want to partition this and only use those dofs that are present on the local partition. To this end, many of the counts in the domain description have been duplicated to include both, the overall number aswell as the counts on the local partition. For the naming of these we use the convention to append _all for variables referring to the overall system across all partitions and _partition for those that refer to local partition variables. The _partition suffix may also be omitted, so names without a suffix usually refer to the partition local quantities.

Nomenclature:

To avoid confusion between the naming of element-local and global systems and the distribution across processes we utilize the postfixes _partition and _all as explained above. We also need to distinguish between the total system with all dofs, including known ones and the effective system, which only enumerates the unknowns. In the laid out logic that means we have total_dof_partition and we have a mapping from that to eff_dof_partition. We also need a mapping from the total_dof_partition to the total_dof_all. And a mapping from the eff_dof_partition to the eff_dof_all. Where the eff_dof_all first has to be found.

Handling of Nodes

Some nodes will be found on partition boundaries. These may have different node types assigned to them from the adjacent nodes, and need to be merged. This requires a communication between partitions to determine the node type to be used. The preliminary approach to deal with nodes is that all processes hold a list for potentially all nodes of the mesh, but only those get filled and evaluated that are adjacent to elements in the partition. This allows direct access to nodes via their global node identification.

Node types are inferred from the adjacent elements, and two neighboring elements may assign different node types to the same node. Thus, we need to have a method to “merge” those nodes and find the proper node type to use. This happens in a first step when reading the nodes on a per element basis. At this point in time, nothing about the connectivity in the mesh is known. In the partitioned approach the nodes on the partition boundaries need to be communicated to achieve the merging.

As we do not know who the neighbors are and who has to communicate with whom, we instead make use of the global node list that we maintain on each process, and perform a reduction with the merge operation for all nodes. This implies a relatively large reduction operation for large numbers of global nodes, and is not the most scalable approach. It also means that each process needs to hold data on all nodes. For a reduced communication effort we will probably have to delay this merging until after building up the connectivity in the mesh and then only communicate with the neighbors. The allreduce for the node merging imposes logarithmic runtime complexity over number of processes.

Nodes that are associated to any elements on the local process, are left as null pointers in the global list of nodes.

Handling of Degrees of Freedom

The degrees of freedom are split into element internal degrees of freedoms and degrees of freedom on the nodes. The element internal degrees of freedom are numbered first, and the nodal DoFs are appended to them. In the partitioned mode, the element internal numberings are performed locally at first and once the number of degrees of freedom in the respective partitions are known, are corrected after computing the prefix sum via MPI_Exscan to reflect the actual global dof numbering across all partitions.

Once the element internal degrees of freedom have been accounted for, the node dofs are defined. First only the local nodes are considered. That means all those which are not null pointers in the list of all nodes on the current process. Then the number of degrees of freedom is collected by an allreduce for the list of all nodes, and the numbering is adapted accordingly with the partial sums of all nodal degrees of freedom.

The node id refers to the local node numbering. .. _Intel’s Threading Building Blocks (oneTBB): https://www.intel.com/content/www/us/en/developer/tools/oneapi/onetbb.html .. _MUMPS: https://mumps-solver.org/index.php .. _OpenMP: https://www.openmp.org/ .. _MPI: https://www.mpi-forum.org/