Skip to content

Parallel assembly

For developers

In the finite element method, contributions evaluated on the element level must be assembled into global matrices and vectors. For an element \(e\) with local-to-global degree-of-freedom map \(\mathcal{I}_e\), the elemental matrix \(\boldsymbol{K}_e\) and vectors \(\boldsymbol{f}_e\) are added to the corresponding entries of the global system:

\[ K_{\mathcal{I}_e(i)\mathcal{I}_e(j)} \mathrel{+}= K_{e,ij}, \qquad f_{\mathcal{I}_e(i)} \mathrel{+}= f_{e,i}. \]

The same principle applies to the internal-force vector, dynamic-force vector, body-force vector and contributions associated with prescribed boundary conditions. Since neighbouring elements usually share degrees of freedom, a direct parallel loop over all elements would allow several threads to update the same global matrix row or vector entry simultaneously. Such concurrent read-modify-write operations constitute data races and may result in lost contributions and non-reproducible solutions.

numgeo therefore uses a race-free interior/interface partitioning strategy. The strategy avoids atomic operations during the assembly of the finite element contributions while retaining direct insertion into the global compressed sparse row (CSR) matrix and the global vectors.

Workload construction

The workload is constructed in two distinct stages. First, the geometrical element-to-node mesh is partitioned with the dual-mesh partitioner of METIS. Passing the compact geometrical connectivity is important because this is the representation for which the METIS mesh-partitioning algorithm is designed. If METIS is unavailable, incompatible with the integer interface used by numgeo or does not return a valid partition, a deterministic weighted fallback partitioner is used.

The geometrical partition alone is not sufficient for race-free direct assembly. The actual write locations are governed by the element-to-global-degree-of-freedom maps, which may contain several fields per node and may also include non-geometrical couplings. numgeo therefore evaluates the exact degree-of-freedom connectivity after the initial partition has been obtained. Each global degree of freedom is assigned deterministically to one owner partition. An element remains in a parallel partition interior only if its initial partition owns every global degree of freedom addressed by that element. All other elements are moved to a separate serial-interface workload.

This owner-based extraction avoids serialising both sides of every partition boundary. At the same time, the remaining interiors of any two partitions have disjoint sets of global degrees of freedom:

\[ \mathcal{D}_{p} \cap \mathcal{D}_{q} = \varnothing \qquad \text{for} \qquad p \neq q, \]

where \(\mathcal{D}_{p}\) denotes the set of global degrees of freedom addressed by the interior elements of partition \(p\). This condition guarantees that two concurrently processed partitions cannot write to the same global matrix row or vector entry.

Assembly sequence

The assembly is performed in two stages:

  1. The partition interiors are processed in parallel using OpenMP. One OpenMP work item corresponds to one complete partition. The elements within a partition are evaluated sequentially by the thread assigned to that partition.
  2. After the parallel region has finished, the interface elements are evaluated and assembled sequentially.

The procedure may be summarised as

construct or validate workload

parallel over partition interiors
    for every element in the partition
        gather global solution variables
        evaluate the element response
        scatter the element matrix and vectors directly to the global system
    end
end parallel

for every interface element
    gather global solution variables
    evaluate the element response
    scatter the element matrix and vectors directly to the global system
end

No atomic update is required for the interior contributions because the partition interiors are degree-of-freedom-disjoint. The serial treatment of the interface guarantees the same property for the remaining contributions. The approach is used for the assembly of the element matrix and the associated internal, dynamic, body-force and boundary-condition vectors. Operations based on the same fixed element connectivity, such as gathering the global internal-force vector, use the same workload decomposition.

Local element evaluation

Each element is evaluated independently. Before calling the element formulation, numgeo gathers the current global solution variables according to the element map \(\mathcal{I}_e\). Element-local arrays are then constructed and the constitutive and element routines are evaluated. If the element requests a material, radius or volume cutback, its contribution is discarded and is not scattered to the global system.

Otherwise, the local matrix and vectors are inserted directly into the known CSR structure. The position of a local column degree of freedom in a global matrix row is determined by a binary search in the sorted column-index array. This avoids temporary thread-local global matrices and the subsequent reduction of such matrices.

Reduction of auxiliary quantities

Some quantities calculated during element evaluation are not assembled additively. Examples are cutback flags and critical time-step criteria. For these quantities, every partition stores a local result. The partition results and the result of the serial interface are combined after the parallel region.

The reduction reproduces the ordering rules of serial assembly. In particular, when several elements yield the same critical value, the element occurring first in the global element order is retained. This explicit tie-breaking prevents the selected controlling element from depending on OpenMP scheduling.

Validity of the workload

The workload depends on both the geometrical mesh connectivity used for the initial partition and the complete element-to-global-degree-of-freedom connectivity used for race-free interface extraction. It is rebuilt together with the sparsity pattern and whenever the mesh dimensions, requested number of threads or generation of the local-to-global maps changes. The generation check is constant-time during repeated element evaluations. A complete connectivity comparison is retained in the one-time workload validation and in debug mode. The validation checks that:

  • every element occurs exactly once, either in one partition interior or in the serial interface;
  • all stored element and degree-of-freedom indices are valid;
  • every parallel interior element belongs to the owner partition of all its global degrees of freedom;
  • no global degree of freedom occurs in two different partition interiors; and
  • the stored connectivity agrees with the current local-to-global element maps.

These checks are important after changes to the mesh, active degrees of freedom or element mappings. A workload generated for an outdated connectivity must never be reused.

Scope and limitations

The partitioning strategy applies to assembly operations whose write locations follow the fixed finite element connectivity. Contributions with connectivity that can change during the analysis, particularly contact contributions, require their own synchronisation strategy and are not assigned to the static element-interface workload.

The achievable speed-up depends on the ratio between interior and interface elements, the computational cost of the element formulations and the balance of the partitions. A very small model, or a partitioning that places a large fraction of the elements on the interface, may provide only limited benefit. Conversely, computationally expensive constitutive models and meshes with comparatively small interfaces generally offer favourable parallel efficiency.