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 parallel workload is constructed from the actual element-to-global-degree-of-freedom mapping used by the solver. Two elements are considered connected if they share at least one global degree of freedom. The initial mesh partition is generated with the dual-mesh partitioner of METIS. If METIS is not available or does not return a valid partition, a deterministic weighted fallback partitioner is used.

The initial partition alone is not sufficient for race-free direct assembly because degrees of freedom on a partition boundary are shared by elements assigned to different partitions. numgeo therefore classifies the elements into two groups:

Partition interiors: elements whose global degrees of freedom occur in one partition only.
Partition interface: elements touching at least one degree of freedom that is also used by another partition.

Every interface element is removed from its initial partition and added to a separate serial-interface workload. Consequently, 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:

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.
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 the complete element-to-global-degree-of-freedom connectivity. It is rebuilt together with the sparsity pattern and whenever the stored workload is no longer compatible with the current model or requested number of threads. 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;
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.