Distributed-memory parallelization of the aggregated unfitted finite element method

TL;DR

This paper presents a scalable distributed-memory implementation of the aggregated unfitted finite element method (AgFEM), enabling efficient large-scale solutions with standard algebraic multigrid preconditioners on complex 3D problems.

Contribution

It introduces a parallel distributed-memory implementation of AgFEM that handles large problem sizes using standard solvers, contrasting with previous customized approaches.

Findings

Successfully solved 300 million degrees of freedom problems on 16,000 cores.

Achieved efficient weak scaling on complex 3D domains.

Demonstrated compatibility with standard algebraic multigrid preconditioners.

Abstract

The aggregated unfitted finite element method (AgFEM) is a methodology recently introduced in order to address conditioning and stability problems associated with embedded, unfitted, or extended finite element methods. The method is based on removal of basis functions associated with badly cut cells by introducing carefully designed constraints, which results in well-posed systems of linear algebraic equations, while preserving the optimal approximation order of the underlying finite element spaces. The specific goal of this work is to present the implementation and performance of the method on distributed-memory platforms aiming at the efficient solution of large-scale problems. In particular, we show that, by considering AgFEM, the resulting systems of linear algebraic equations can be effectively solved using standard algebraic multigrid preconditioners. This is in contrast with previous works that consider highly customized preconditioners in order to allow one the usage of iterative solvers in combination with unfitted techniques. Another novelty with respect to the methods available in the literature is the problem sizes that can be handled with the proposed approach. While most of previous references discussing linear solvers for unfitted methods are based on serial non-scalable algorithms, we propose a parallel distributed-memory method able to efficiently solve problems at large scales. This is demonstrated by means of a weak scaling test defined on complex 3D domains up to 300M degrees of freedom and one billion cells on 16K CPU cores in the Marenostrum-IV platform. The parallel implementation of the AgFEM method is available in the large-scale finite element package FEMPAR.

Tables

Table 1. Table 1 : Summary of the assembly operations related with the agFE space 𝒱 agg superscript 𝒱 agg \mathcal{V}^{\rm agg} . These assembly operations are related with the methodology introduced in [ 44 ] .

Case	Assembly operation
Right hand side vector
$g^{a,k}\in {ℐ}^{\mathrm{in}}$	${\left[𝐛\right]}_i\text{ += }{\left[{𝐛}_k\right]}_a$, $i=g^{a,k}$  (standard case)
$g^{a,k}\notin {ℐ}^{\mathrm{in}}$	${\left[𝐛\right]}_i\text{ += }{C}_{ni}{\left[{𝐛}_k\right]}_a$, $n=g^{a,k}$, $\forall i\in ℳ\left(n\right)$
System matrix
$g^{a,k}\in {ℐ}^{\mathrm{in}}$ and $g^{b,k}\in {ℐ}^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{\left[{𝐀}_k\right]}_{ab}$, $i=g^{a,k}$, $j=g^{b,k}$  (standard case)
$g^{a,k}\notin {ℐ}^{\mathrm{in}}$ and $g^{b,k}\in {ℐ}^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{C}_{ni}{\left[{𝐀}_k\right]}_{ab}$, $n=g^{a,k}$, $j=g^{b,k}$, $\forall i\in ℳ\left(n\right)$
$g^{a,k}\in {ℐ}^{\mathrm{in}}$ and $g^{b,k}\notin {ℐ}^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{C}_{mj}{\left[{𝐀}_k\right]}_{ab}$, $i=g^{a,k}$ , $m=g^{b,k}$, $\forall j\in ℳ\left(m\right)$
$g^{a,k}\notin {ℐ}^{\mathrm{in}}$ and $g^{b,k}\notin {ℐ}^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{C}_{ni}{C}_{mj}{\left[{𝐀}_k\right]}_{ab}$, $n=g^{a,k}$ , $m=g^{b,k}$, $\forall i\in ℳ\left(n\right)$, $\forall j\in ℳ\left(m\right)$

Table 2. Table 2 : Summary of the assembly operations related with the agFE space 𝒱 agg superscript 𝒱 agg \mathcal{V}^{\rm agg} in the distributed-memory implementation.

Case	Assembly operation
Right hand side vector
$j_{s,\ell }^a\in {𝒥}_s^{\mathrm{in}}$	${\left[𝐛\right]}_i\text{ += }{\left[{𝐛}_s^{\ell }\right]}_a$, $i=g_{s,\ell }^a$  (standard case)
$j_{s,\ell }^a\notin {𝒥}_s^{\mathrm{in}}$	${\left[𝐛\right]}_i\text{ += }{C}_s^{ni}{\left[{𝐛}_s^{\ell }\right]}_a$, $n=j_{s,\ell }^a$, $\forall i\in {ℳ}_s\left(n\right)$
System matrix
$j_{s,\ell }^a\in {𝒥}_s^{\mathrm{in}}$ and $j_{s,\ell }^b\in {𝒥}_s^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{\left[{𝐀}_s^{\ell }\right]}_{ab}$, $i=g_{s,\ell }^a$, $j=g_{s,\ell }^b$  (standard case)
$j_{s,\ell }^a\notin {𝒥}_s^{\mathrm{in}}$ and $j_{s,\ell }^b\in {𝒥}_s^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{C}_s^{ni}{\left[{𝐀}_s^{\ell }\right]}_{ab}$, $n=j_{s,\ell }^a$, $j=g_{s,\ell }^b$, $\forall i\in {ℳ}_s\left(n\right)$
$j_{s,\ell }^a\in {𝒥}_s^{\mathrm{in}}$ and $j_{s,\ell }^b\notin {𝒥}_s^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{C}_s^{mj}{\left[{𝐀}_s^{\ell }\right]}_{ab}$, $i=g_{s,\ell }^a$ , $m=j_{s,\ell }^b$, $\forall j\in {ℳ}_s\left(m\right)$
$j_{s,\ell }^a\notin {𝒥}_s^{\mathrm{in}}$ and $j_{s,\ell }^b\notin {𝒥}_s^{\mathrm{in}}$	${\left[𝐀\right]}_{ij}\text{ += }{C}_s^{ni}{C}_s^{mj}{\left[{𝐀}_s^{\ell }\right]}_{ab}$, $n=j_{s,\ell }^a$, $m=j_{s,\ell }^b$, $\forall i\in {ℳ}_s\left(n\right)$, $\forall j\in {ℳ}_s\left(m\right)$

Table 3. Table 3 : Number of sub-domains ( | 𝒮 | 𝒮 |\mathcal{S}| ) and number of total cells in the background mesh for the three families of partitioned meshes considered in the examples. In parentheses, the percentage of active cells (over the total number of background cells) for each of the three geometries is given. Within a parenthesis, the three given quantities correspond to the "pop-corn flake", the spiral, and the "Swiss cheese" geometries, respectively.

	Number of cells
$\left|𝒮\right|$	load 1	load 2	load 3
4	4,096 (34% 40%    -   )	32,768 (29% 29% 56%)	262,144 (26% 24% 47%)
32	32,768 (29% 29% 56%)	262,144 (26% 24% 47%)	2,097,152 (25% 22% 41%)
262	262,144 (26% 24% 47%)	2,097,152 (25% 22% 41%)	16,777,216 (24% 20% 38%)
2,097	2,097,152 (25% 22% 41%)	16,777,216 (24% 20% 38%)	134,217,728 (24% 20% 36%)
16,777	16,777,216 (24% 20% 38%)	134,217,728 (24% 20% 36%)	1,073,741,824 (24% 20% 35%)

Table 4. Table 4 : Summary of the main parameters and computational strategies used in the numerical examples.

Description	Considered methods/values
Model problem	Poisson equation
Problem geometry	"Popcorn flake", massive spiral, and "Swiss cheese"
Parallel mesh generation and partitioning	p4est library
FE spaces	Aggregated ${𝒱}^{\mathrm{agg}}$ and standard ${𝒱}^{\mathrm{std}}$
Cell type	Hexahedral cells
Interpolation	Piece-wise trilinear shape functions
Linear solver	Preconditioned conjugate gradients (PETSc)
Preconditioner	Smoothed-aggregation AMG (GAMG preconditioner from PETSc)
Weighting factor for space filling-curve	$w=10$
Coefficient used in Nitsche’s penalty terms	$\beta =10.0$ and $\beta =100.0$

Table 5. Table 5 : Weak scaling test for the Poisson equation on the unit cube: Number of sub-domains ( | 𝒮 | 𝒮 |\mathcal{S}| ) and number of cells of the considered meshes.

	Number of cells
$\left|𝒮\right|$	load 1	load 2	load 3
48	384,000	1,296,000	3,072,000
1,296	10,368,000	34,992,000	82,944,000
3,072	24,576,000	82,944,000	196,608,000
6,000	48,000,000	162,000,000	384,000,000
10,368	82,944,000	279,936,000	663,552,000
16,464	131,712,000	444,528,000	1,053,696,000

Equations

$${\mathcal{V}}^{\mathrm{std}}\doteq \{v\in C^0({\Omega }^{\mathrm{act}}):v{|}_T\in \mathcal{V}(T)\text{ for any }T\in \mathcal{T}\}.$$

$${\mathcal{V}}^{\mathrm{agg}}\doteq \{v\in {\mathcal{V}}^{\mathrm{std}}:v^i=\sum _{j\in \mathcal{M}(i)}{C}_{ij}{v}^j\text{ for any }i\in {\mathcal{I}}^{\mathrm{out}}\},$$

$$\tilde{v}\doteq \sum _{j\in {\mathcal{I}}^{\mathrm{in}}}v(x^j){\varphi }^j+\sum _{i\in {\mathcal{I}}^{\mathrm{out}}}\sum _{j\in \mathcal{M}(i)}{C}_{ij}v(x^j){\varphi }^i.$$

$${\tilde{\varphi }}^g\doteq \sum _{j\in {\mathcal{I}}^{\mathrm{in}}}{\varphi }^g(x^j){\varphi }^j+\sum _{i\in {\mathcal{I}}^{\mathrm{out}}}\sum _{j\in \mathcal{M}(i)}{C}_{ij}{\varphi }^g(x^j){\varphi }^i.$$

$${\tilde{\varphi }}^g\doteq {\varphi }^g+\sum _{i\in \mathcal{Z}(g)}{C}_{ig}{\varphi }^i,$$

$${\left[\mathbf{b}\right]}_i=\mathrm{b}({\varphi }^i)+\sum _{n\in \mathcal{Z}(i)}{C}_{ni}\mathrm{b}({\varphi }^n),\text{for }i\in {\mathcal{I}}^{\mathrm{in}},\text{ and}$$

$${\left[\mathbf{A}\right]}_{ij}=\mathrm{a}({\varphi }^i,{\varphi }^j)+\sum _{n\in \mathcal{Z}(i)}{C}_{ni}\mathrm{a}({\varphi }^n,{\varphi }^j)+\sum _{m\in \mathcal{Z}(j)}{C}_{mj}\mathrm{a}({\varphi }^i,{\varphi }^m)+\sum _{n\in \mathcal{Z}(i)}\sum _{m\in \mathcal{Z}(j)}{C}_{ni}{C}_{mj}\mathrm{a}({\varphi }^n,{\varphi }^m),$$

$$[\mathbf{b}{]}_i=\sum _{k\in \mathcal{B}(i)}[{\mathbf{b}}_k{]}_a+\sum _{n\in \mathcal{Z}(i)}\sum _{k\in \mathcal{B}(n)}{C}_{ni}[{\mathbf{b}}_k{]}_b\text{ for }i\in {\mathcal{I}}^{\mathrm{in}},$$

$${\mathcal{V}}_s^{\mathrm{agg}}=\left\{v_s\in {\mathcal{V}}_s^{\mathrm{std}}:v_s^j=\sum _{a\in \Lambda (T_{s^{\prime }}^{{\ell }^{\prime }})}{C}_s^{ja}{v}_{s^{\prime },{\ell }^{\prime }}^a\text{ for all }j\in {\mathcal{J}}_s^{\mathrm{out}}\right\},$$

$${\mathcal{V}}_s^{\mathrm{agg}}\doteq \{v_s\in {\mathcal{V}}_s^{\mathrm{std}}:v_s^j=\sum _{g\in {\mathcal{M}}_s(j)}{C}_s^{jg}[v^g{]}_s\text{ for any }j\in {\mathcal{J}}_s^{\mathrm{out}}\},$$

$$\left.\begin{array}[]{l}\mathcal{M}_{s}(j)\doteq\{g^{a}_{s,\ell}\}_{a\in\Lambda(T^{\ell}_{s})},\quad C^{jg}_{s}\doteq\phi^{b}_{s,\ell}(x^{j}_{s}),\quad[v^{g}]_{s}\doteq v^{b}_{s,\ell}\quad\text{ if }k\in\mathcal{K}_{s}\\ \mathcal{M}_{s}(j)\doteq\{\tilde{g}^{a}_{s,z}\}_{a\in\Lambda(\tilde{T}^{z}_{s})},\quad C^{jg}_{s}\doteq\tilde{\phi}^{b}_{s,z}(x^{j}_{s}),\quad[v^{g}]_{s}\doteq\tilde{v}^{b}_{s,z}\quad\text{ otherwise}\end{array}\right.$$

$$-\Delta u=f\text{in }\Omega ,u=g\text{on }\Gamma \doteq \partial \Omega ,$$

$$\begin{array}[]{l}\displaystyle\mathrm{a}(v,u)\doteq\int_{\Omega}\nabla v\cdot\nabla u\mathrm{\ d}\Omega+\int_{\Gamma}\left(\tau vu-v\left(n\cdot\nabla u\right)-u\left(n\cdot\nabla v\right)\right)\mathrm{\ d}{\Gamma},\text{ and}\\ \displaystyle\mathrm{b}(v)\doteq\int_{\Omega}vf\mathrm{\ d}\Omega+\int_{\Gamma}\left(\tau vg-\left(n\cdot\nabla v\right)g^{\mathrm{D}}\right)\mathrm{\ d}\Gamma,\end{array}$$

$$\tau =\frac{\beta }{h},$$

$${\int }_{T^k\cap \Omega }\nabla {\eta }_k\cdot \nabla {\xi }_k\mathrm{d}\Omega ={\lambda }_k{\int }_{T^k\cap \Gamma }(\nabla {\eta }_k\cdot n)(\nabla {\xi }_k\cdot n)\mathrm{d}\Gamma \text{ for all }{\xi }_k\in \mathcal{V}(T^k)\text{ for each cut cell id }k\in {\mathcal{K}}^{\mathrm{cut}}.$$

$$|u{|}_{H^1}\doteq {\left({\int }_{\Omega }\nabla u\cdot \nabla u\mathrm{d}\Omega \right)}^{1/2}\text{ and }\Vert u{\Vert }_{L^2}\doteq {\left({\int }_{\Omega }{u}^2\mathrm{d}\Omega \right)}^{1/2}$$

Full text

DISTRIBUTED-MEMORY PARALLELIZATION OF THE AGGREGATED UNFITTED FINITE ELEMENT METHOD

Francesc Verdugoa,***Corresponding author.

Emails: fverdugo@cimne.upc.edu (FV), amartin@cimne.upc.edu (AM), sbadia@cimne.upc.edu (SB)

, Alberto F. Martína,b, and Santiago Badiaa,b

aCIMNE – Centre Internacional de Mètodes Numèrics en Enginyeria

Esteve Terradas 5, 08860 Castelldefels, Spain.

bDepartment of Civil and Environmental Engineering, Universitat Politècnica de Catalunya

Jordi Girona 1-3, Edifici C1, 08034 Barcelona, Spain.

Abstract

The aggregated unfitted finite element method (AgFEM) is a methodology recently introduced in order to address conditioning and stability problems associated with embedded, unfitted, or extended finite element methods. The method is based on removal of basis functions associated with badly cut cells by introducing carefully designed constraints, which results in well-posed systems of linear algebraic equations, while preserving the optimal approximation order of the underlying finite element spaces. The specific goal of this work is to present the implementation and performance of the method on distributed-memory platforms aiming at the efficient solution of large-scale problems. In particular, we show that, by considering AgFEM, the resulting systems of linear algebraic equations can be effectively solved using standard algebraic multigrid preconditioners. This is in contrast with previous works that consider highly customized preconditioners in order to allow one the usage of iterative solvers in combination with unfitted techniques. Another novelty with respect to the methods available in the literature is the problem sizes that can be handled with the proposed approach. While most of previous references discussing linear solvers for unfitted methods are based on serial non-scalable algorithms, we propose a parallel distributed-memory method able to efficiently solve problems at large scales. This is demonstrated by means of a weak scaling test defined on complex 3D domains up to 300M degrees of freedom and one billion cells on 16K CPU cores in the Marenostrum-IV platform. The parallel implementation of the AgFEM method is available in the large-scale finite element package FEMPAR.

Keywords: Unfitted finite element methods $\cdot$ Algebraic multigrid $\cdot$ High performance scientific computing

1 Introduction

Unfitted finite element (FE) methods (also known as embedded or immersed FE methods) are useful techniques in order to perform numerical computations associated with geometrically complex domains. They are specially useful in multi-phase and multi-physics applications with moving interfaces (e.g., fracture mechanics [1], fluid–structure interaction [2], free surface flows [3]), in applications with varying domain topologies (e.g., shape or topology optimization [4], 3D printing simulations [5]), in applications where the geometry is not described by CAD data (e.g., medical simulations based on CT-scan images [6]), or in large-scale parallel computations, where generating and partitioning large unstructured meshes is particularly difficult. The main benefit of unfitted FE methods is that they do not require the generation of body-fitted meshes (i.e., meshes whose faces conform to the domain boundary). Instead, they embed the domain of interest in a geometrically simple background grid (usually a uniform or an adaptive Cartesian grid), which can be generated much more efficiently. The popularity of this approach is illustrated by the large number of different methodologies that have been proposed following this rationale (see, e.g., the cutFEM method [7], the cutIGA method [8], the Finite Cell Method [9], the AgFEM method [10], the cg-FEM method [11], the Immersed Boundary Method [12], the i-spline method [13], and some variants of the XFEM method [14]).

Unfortunately, unfitted FE methods have also well known drawbacks. One of the most notorious ones is the so-called small cut cell problem. The intersection of a background cell with the physical domain can be arbitrarily small and with arbitrarily high aspect ratios, which usually results in severely ill-conditioned systems of algebraic linear equations if no specific strategy is used to remedy it [15]. This flaw makes the solution of the underlying linear systems much more challenging than for standard FE methods based on body-fitted grids and it is still today one of the main limiting factors for the successful application of unfitted FE techniques in realistic large-scale applications. Sparse direct solvers [16] are usually considered when dealing with unfitted FE techniques (see, e.g., [9]) since they are robust enough to deal (up to a certain extent) with such ill-conditioned problems. However, its usage is prohibitive at large scales since their memory footprint and the algorithmic complexity scales supra-linearly with respect to the problem size. In FE analysis, the current way to effectively solve linear systems at large-scales is using iterative Krylov sub-space methods [17] in combination with parallel and scalable preconditioners. Unfortunately, the well known scalable preconditioners based on multigrid [18, 19] (either geometric [20] or algebraic [21, 22]), or multi-level domain decomposition [23, 24] are mainly designed for body-fitted meshes and cannot readily deal with the ill-conditioning associated with the small cut cell problem.

Specific preconditioners for unfitted methods have been proposed in the literature in different contexts. These techniques are generally based on the idea of splitting the system matrix into blocks associated with problematic degrees of freedom, i.e., the ones affected by small cuts, and unproblematic ones, i.e., the ones not affected (see, e.g., [25, 15, 26, 27, 28]). The methods differ in the particular techniques used to handle the block associated with non-problematic DOFs, the block of problematic DOFs, and the coupling between these two blocks. Even though these methodologies have been able to allow the usage of iterative solvers in combination with unfitted FE methods in particular settings, they are mainly serial non-scalable algorithms. Considered problems are typically academic examples that go up to few hundreds of thousands of DOFs. To the authors’ best knowledge, the only exceptions are the works in [29, 27]. Reference [27] is based on a single level Additive Schwarz preconditioner and reports results for problems up to almost 4 million DOFs solved with a hybrid OpenMP+MPI implementation with up to 16 MPI tasks and 64 OpenMP threads per MPI task. Even though this method is a step forward with respect to most previous works, it is well known that single-level Additive Schwarz preconditioners do not scale up algorithmically, and thus the method is not suitable for larger scales. On the other hand, the work in [29] considers a two-level balancing domain decomposition preconditioner [30] and reports results up to almost 6 million DOFs.

The goal of this work is to scale up to much larger problem sizes by considering an alternative approach. Instead of using a tailored preconditioner able to deal with linear systems affected by the small cut cell problem, we consider an enhanced FE formulation so that the resulting method leads to system matrices that are not affected by small cuts. This opens the door to consider well known scalable algorithms for conventional FE analysis such as algebraic multigrid (AMG) [21, 22], for which there are highly scalable parallel implementations in renowned scientific computing packages such as TRILINOS [31] or PETSc [32, 33]. This makes the method easy to use in parallel runs since the development of tailored linear solvers is not required. To this end, we consider the recently introduced aggregated unfitted finite element method (AgFEM) [10]. Other enhanced unfitted formulations could be also taken into account (e.g., the CutFEM method [7]), but this task is not in the scope of the current work. AgFEM is based on removal of basis functions associated with badly cut cells by introducing carefully designed constraints and its formulation shares the good properties of body-fitted FE methods such as stability, condition number bounds, optimal convergence, and continuity with respect to data. The detailed mathematical analysis of the method is included in [10] for elliptic problems and in [34] for the Stokes equation. Even though the method is potentially applicable to parallel large-scale computations, previous works [10, 34] focus mainly in the description of the methodology in a serial context and its mathematical analysis. The main contribution of the current work is this two-fold. On the one hand, we extend AgFEM to large-scale computations by presenting an efficient distributed-memory implementation, which is based on standard functionality provided by distributed-memory FE packages. In this case, we have implemented the method in the large-scale FE software package FEMPAR  [35]. The algorithmic strategies proposed for AgFEM in distributed-memory environments can also be used for the parallelization of other cell aggregation techniques for discontinuous Galerkin (DG) methods [36, 37, 38]. It will also be essential in CutFEM [7], when using substructuring domain decomposition techniques to end up with well-posed subdomain problems. On the other hand, we show that, by considering the parallel AgFEM, the resulting systems of linear algebraic equations can be effectively solved at large scales using standard AMG preconditioners. Here, we have considered the well-known PETSc’s native AMG framework, referred to as GAMG [39]. We have customized as less as possible the parameters of this solver in order to demonstrate that users can take it as it is given “out of the box” by the library. These properties are demonstrated with weak scaling tests up to 16K cores and up to nearly 300M DOFs in the Marenostrum-IV platform at the Barcelona Supercomputing Center [40] using the Poisson equation on complex 3D domains as a model problem.

The outline of the paper is as follows. In Sect. 2, we briefly present the main components of the AgFEM method in a serial context as a basis for the subsequent extension to distributed-memory platforms. In particular, we introduce the geometrical setup (Sects. 2.1 and 2.2), the FE spaces of the AgFEM method (Sect. 2.3), the associated shape function basis (Sect. 2.4), and the FE assembly (Sect. 2.5). In Sect. 3, we discuss the distributed-memory implementation. We detail the domain-decomposition setup (Sect. 3.1), the distributed FE spaces (Sect. 3.2), several important algorithmic components (Sects. 3.3, 3.4, 3.5, and 3.6), and the parallel FE assembly (Sect. 3.7). In Sect. 4, we report the numerical examples. We detail the HPC environment (Sect. 4.1), the considered model problem (Sect. 4.2), the considered linear solver (Sect. 4.3), the particular setup for the experiments (Sect. 4.4), and the results of the weak scaling test (Sects. 4.5 and 4.6). Finally, we close the paper with the conclusions in Sect. 5.

2 The aggregated unfitted finite element method

2.1 Immersed boundary setup

Let $\Omega\subset\mathbb{R}^{d}$ be an open bounded polygonal domain with $d\in\{2,3\}$ the number of spatial dimensions. Like in any other embedded boundary method, let us introduce an artificial domain $\Omega^{\rm art}$ such that it has a simple geometry (e.g., a cuboid) and it includes the physical domain $\Omega\subseteq\Omega^{\rm art}$ (see Fig.  1(a)). We denote by $\mathcal{T}^{\rm art}$ a partition of $\Omega^{\rm art}$ into cells with characteristic cell size $h$. We build $\mathcal{T}^{\rm art}$ as a Cartesian grid of $d$-dimensional parallelepipeds, i.e., hexahedra for $d=3$ and quadrilaterals for $d=2$ (see Fig.  1(b)). We consider Cartesian grids since they can be efficiently generated and partitioned in parallel with the p4est library [41]. Note, however, that the following rationale holds also for unstructured background meshes and other cell types such as triangles and tetrahedra.

We introduce the usual classification of cells in $\mathcal{T}^{\rm art}$ as internal, external or cut. Specifically, $T\in\mathcal{T}^{\rm art}$ is an internal cell if $T\subseteq\Omega$; $T$ is an external cell if $T\cap\Omega=\emptyset$; otherwise, $T$ is a cut cell (see Fig. 1(b)). For the sake of simplicity and without any loss of generality, we represent the domain boundary $\partial\Omega$ as the zero level-set of a given scalar function $\psi^{\rm ls}$, namely $\partial\Omega\doteq\{x\in\mathbb{R}^{d}:\psi^{\rm ls}(x)=0\}$. We define the physical domain as the set of points where the levelset function is negative, namely $\Omega\doteq\{x\in\mathbb{R}^{d}:\psi^{\rm ls}(x)<0\}$. In this setup, classifying the cells into internal, external or cut is straightforward by sampling the levelset function at cell nodes. We note that the problem geometry could be described by other means, e.g., using 3D CAD data, by providing techniques to perform the cell classification and to perform numerical integration in cut cells (see, e.g., [42]). In any case, the way the geometry is handled does not affect the following exposition. We represent the set of interior (resp., external and cut) cells with $\mathcal{T}^{\rm in}$ (resp., $\mathcal{T}^{\rm out}$ and $\mathcal{T}^{\rm cut}$). Furthermore, we define the set of active cells as the union of the sets of interior and cut cells, namely $\mathcal{T}^{\rm act}\doteq\mathcal{T}^{\rm in}\cup\mathcal{T}^{\rm cut}$. The domain associated with active cells, namely $\Omega^{\rm act}\doteq\left(\bigcup_{T\in\mathcal{T}^{\rm act}}T\right)\cup(\bigcup_{T,T^{\prime}\in\mathcal{T}^{\rm act}}{}\mkern 1.0mu\overline{\mkern-1.0muT}\cap{}\mkern 1.0mu\overline{\mkern-1.0muT}^{\prime})$ is referred to as the active domain, which satisfies $\Omega\subseteq\Omega^{\rm act}\subseteq\Omega^{\rm art}$. For simplicity, we refer to the set of active cells $\mathcal{T}^{\rm act}$ simply as $\mathcal{T}$ (i.e., we drop the label "act" in the notation). In the following developments, it is convenient to assign a unique global identifier (referred to as global id) to each cell in $\mathcal{T}$. We denote by $T^{k}$ the cell with the global id $k$, for $k=1,\ldots,|\mathcal{T}|$, where $|\cdot|$ is the number of elements in a set. With this notation, we can write $\mathcal{T}=\{T^{k}\}_{k\in\mathcal{K}}$, where $\mathcal{K}\doteq\{1,\ldots,|\mathcal{T}|\}$ is the index set of all global cell ids. We also introduce $\mathcal{K}^{\rm in}$ (resp. $\mathcal{K}^{\rm cut}$) the subset of $\mathcal{K}$ that contains the ids of all interior (resp. cut) cells, namely $\mathcal{K}^{\rm in}\doteq\{k\in\mathcal{K}:T^{k}\in\mathcal{T}^{\rm in}\}$ and $\mathcal{K}^{\rm cut}\doteq\{k\in\mathcal{K}:T^{k}\in\mathcal{T}^{\rm cut}\}$.

2.2 Cell aggregation

We briefly present here the (serial) cell aggregation strategy proposed in [10], which is the cornerstone of the AgFEM method. The extension of this cell aggregation method to parallel distributed-memory frameworks is discussed later in Sect. 3.4. The cell aggregation process consists in building a surjective map that for each cell in $\mathcal{T}$ returns an interior cell in $\mathcal{T}^{\rm in}$. For convenience, we represent this map as a transformation $R:\mathcal{K}\rightarrow\mathcal{K}^{\rm in}$ between cell ids. That is, $R$ takes a cell id $k\in\mathcal{K}$ and returns cell id $R(k)\in\mathcal{K}^{\rm in}$ corresponding to an interior cell. The transformation $R$ is called the root cell map since, for a given cell $T^{k}\in\mathcal{T}$, the associated interior cell $T^{R(k)}\in\mathcal{T}^{\rm in}$ is referred to as the root cell of $T^{k}$. The concept of root cell plays an important role in the construction of the so-called aggregated finite element (agFE) spaces, which will be discussed later in Sect 2.3.

The root cell map $R$ induces a partition of the cells in $\mathcal{T}$ into so-called cell aggregates. For each interior cell $T^{k}\in\mathcal{T}^{\rm in}$, we define a cell aggregate $\mathcal{A}_{k}$ as the set of cells that share $T^{k}$ as their root cell. Formally, the aggregate associated with cell $T^{k}$ is defined as $\mathcal{A}_{k}\doteq\{T^{k^{\prime}}\}_{k^{\prime}\in R^{-1}(k)}$, where the inverse of $R$ is defined as $R^{-1}(k)\doteq\{k^{\prime}\in\mathcal{K}:\ k=R(k^{\prime})\}$. Clearly, the set of aggregates $\{\mathcal{A}_{k}\}_{k\in\mathcal{K}^{\rm in}}$ is a partition of $\mathcal{T}$ since each cell in $\mathcal{T}$ belongs to one (and only one) aggregate.

As stated in [10], in order to construct AgFEM spaces, the map $R$ has to be such that the resulting cell aggregates $\mathcal{A}_{k}$ satisfy the following two properties for each $k\in\mathcal{K}^{\rm in}$: (a) the set $\mathcal{A}_{k}$ contains one and only one interior cell, namely $|\mathcal{A}_{k}\cap\mathcal{T}^{\rm in}|=1$; and (b) the (unique) interior cell $\tilde{T}_{0}\in\mathcal{A}_{k}\cap\mathcal{T}^{\rm in}$ can be reached from any cut cell $\tilde{T}_{n}\in\mathcal{A}_{k}\cap\mathcal{T}^{\rm cut}$ in the aggregate following a cell path $\tilde{T}_{0}$, $\tilde{T}_{1}$, $\ldots$, $\tilde{T}_{n}$, where all cells are in the aggregate, i.e., $\tilde{T}_{i}\in\mathcal{A}_{k}$, $i=0,\ldots,n$, and two consecutive cells, $\tilde{T}_{i-1}$ and $\tilde{T}_{i}$, share a non-external face (for $d=3$) or a non-external edge (for $d=2$) for any $i=1,\ldots,n$. More precisely, we say that $\tilde{T}_{i-1}$ and $\tilde{T}_{i}$ share a non-external face ($d=3$) or a non-external edge ($d=2$) if $\tilde{T}_{i-1}\cap\tilde{T}_{i}$ is a manifold in $\mathbb{R}^{d}$ with co-dimension 1 and $\tilde{T}_{i-1}\cap\tilde{T}_{i}\cap\Omega\neq\emptyset$. Note that sharing only a corner (or an edge for $d=3$) is not enough. For simplicity, we will use the words "edge" and "face" as completely equivalent synonyms in the following for the two-dimensional case ($d=2$).

In serial runs, a map $R$ fulfilling these properties is built with the strategy described in Alg. 1, and illustrated in Fig. 2 for a simple 2D case. The notation $\mathcal{K}^{\rm nei}(k)$ in line 1 of the algorithm denotes, for a given cell $T^{k}$ with id $k\in\mathcal{K}$, the set containing the ids of the cells that share a face with cell $T^{k}$. In the following, we often denote the cells with ids in $\mathcal{K}^{\rm nei}(k)$ as the face neighbors of cell $T^{k}$. Note that the algorithm consists in 3 basic steps (see Fig.  2). (a) We initialize the value of $R$ for interior cells by setting each interior cell as its own root cell, and we mark all interior cells as "touched". (b) We perform a loop in cut cells. For a given cut cell $T^{k}$, we define its root cell taking the root cell of one of the face neighbors of $T^{k}$ that have been already marked as "touched" (if such a neighbor exists). (c) After the loop, we mark all cells that have been assigned with a root cell id as "touched". Then, steps (b) and (c) are repeated as many times as needed until all the cells are assigned with a root cell id.

In line 1 of the algorithm, it is needed to choose one among the (possibly) several face neighbors that have already been "touched". Several criteria can be considered for this purpose. In our implementation, we chose the neighbor having the closest root cell (in terms of the euclidean distance between cell barycenters). In case of a tie between two or more neighbors, we chose one arbitrarily (e.g., the neighbor with smaller cell id). During the algorithm, we also populate a map $N:\mathcal{K}\rightarrow\mathcal{K}$ that stores the ids of the chosen face neighbor. More precisely, for a given cut cell $T^{k}$ with id $k\in\mathcal{K}^{\rm cut}$, $N(k)$ returns the id of face neighbor selected in line 1, whereas, for an interior cell $T^{k}$ with $k\in\mathcal{K}^{\rm in}$, $N(k)$ returns its own cell id, namely $N(k)=k$. The map $N$ plays an important role in the distributed-memory implementation. Note that the aggregation scheme can be easily applied to arbitrary spatial dimensions. It only requires to find face neighbors and a classification of cells as internal, external, or cut. These procedures are widely available in embedded domain FE codes.

2.3 Standard and aggregated finite element spaces

Using the root cell map $R$, the authors in [10] define agFE spaces, which are the core of the AgFEM method. The aim of this section is to present such agFE spaces in a serial context, for its later extension to the distributed-memory case. To this end, we first define the usual conforming FE space associated with the active part of the mesh, namely

[TABLE]

We denote by $\mathcal{V}(T)$ a vector space of functions defined on a generic cell $T\in\mathcal{T}$. In the numerical examples, we consider hexahedral cells and tri-linear shape functions. That is, we define $\mathcal{V}(T)\doteq\mathcal{Q}_{q}(T)$ as the space of polynomials that are of degree less or equal to $q$ ($q=1$ in the examples) with respect to each variable in $(x_{1},\ldots,x_{d})\in T$. Even though $d$-dimensional simplices are not considered in the numerical experiments, they can be also used in the presented method. In this case, we would define $\mathcal{V}(T)\doteq\mathcal{P}_{q}(T)$ as the space of polynomials of order less or equal to $q$ in the variables $(x_{1},\ldots,x_{d})\in T$. In order to simplify notation, we define the elemental functional spaces $\mathcal{V}(T)$ in the physical cell $T\subset\Omega^{\rm act}$ (even though in the code we use reference parametric spaces as usual). The space $\mathcal{V}^{\rm std}$ is the standard functional space used in unfitted FE methods (see, e.g., [29, 15, 9]). It is well known that it leads to arbitrary ill conditioned systems of linear algebraic equations (if no extra technique is used to remedy it). This is the main motivation in order to introduce the agFE spaces (see, e.g., [10, 34]), since they lead to linear systems whose condition number is not affected by small cuts.

In order to introduce the agFE spaces we need some extra notation. For the sake of simplicity in the exposition and without loss of generality, we restrict ourselves to scalar-valued Lagrangian FEs. The extension to vector-valued or tensor-valued Lagrangian elements is straightforward by applying the same approach component by component. The extension to other types of FEs (e.g., Nédélec elements [43]) is also possible. It only requires to work with the specific DOFs of the element instead of the Lagrangian nodal values.

Let $\{x^{a,k}\}_{a\in\Lambda(T^{k})}$ be the set containing the (physical) coordinates of the Lagrangian nodes associated with cell $T^{k}\in\mathcal{T}$. In particular, $x^{a,k}\in{}\mkern 1.0mu\overline{\mkern-1.0mu\Omega^{\rm act}}$ denotes the coordinate vector of node with id $a$ (in the local scope of cell $T^{k}$), and $\Lambda(T^{k})\doteq\{1,\ldots,|\mathcal{V}(T^{k})|\}$ is the set containing all local node ids in $T^{k}$. Since (for scalar-valued Lagrangian spaces) each local node is associated with one DOF of the space $\mathcal{V}(T^{k})$, we also refer to the local indices $a\in\Lambda(T^{k})$ as the local DOF ids of cell $T^{k}$ and $\Lambda(T^{k})$ as the set of local DOF ids in $T^{k}$. We also need to introduce the set $\{\phi^{a,k}\}_{a\in\Lambda(T^{k})}$ of Lagrangian shape functions in the (physical) cell $T^{k}$. Following an analogous notation as for the cell nodal coordinates, the symbol $\phi^{a,k}$ refers to the shape function associated with the local DOF id $a\in\Lambda(T^{k})$ in cell $T^{k}$. Clearly, $\phi^{a,k}(x^{b,k})=\delta_{ab}$ for any $a,b\in\Lambda(T^{k})$, and $v^{k}=\sum_{a\in\Lambda(T^{k})}v^{a,k}\phi^{a,k}$, with $v^{a,k}\doteq v^{k}(x^{a,k})$ being the cell nodal / DOF values for a given function $v^{k}\in\mathcal{V}(T^{k})$.

On the other hand, we denote by $\mathcal{X}\doteq\{x^{i}\}_{i\in\mathcal{I}}$ and $\{\phi^{i}\}_{i\in\mathcal{I}}$ the set of global nodal coordinates and the set of global shape functions of space $\mathcal{V}^{\rm std}$, resp., where $x^{i}$ and $\phi^{i}$ are the coordinate vector and shape function of the mesh node with global id $i$ and $\mathcal{I}\doteq\{1,\ldots,|\mathcal{V}^{\rm std}|\}$ is the set of global node ids / DOF ids associated with the standard FE space $\mathcal{V}^{\rm std}$. Clearly $\phi^{i}(x^{j})=\delta_{ij}$, for any $i,j\in\mathcal{I}$, and $v=\sum_{i\in\mathcal{I}}v^{i}\phi^{i}$, with $v^{i}\doteq v(x^{i})$ being the global nodal / DOF values for given function $v\in\mathcal{V}^{\rm std}$ .

We also need to introduce, for each cell $T^{k}\in\mathcal{T}$, the usual transformation between the cell local DOF ids and the corresponding global DOF ids. To this end, we denote by $g^{a,k}\in\mathcal{I}$ the global DOF id associated with the local DOF id $a$ in cell $T^{k}$. Using the local-to-global map described by the values $g^{a,k}$, we define the set of global DOF ids touched by a cell $T^{k}$, i.e., $\mathcal{G}(k)\doteq\{g^{a,k}\}_{a\in\Lambda(T^{k})}$. In addition, we introduce the notation $[\phi^{i}]_{k}$ representing the local shape function in cell $T^{k}$ associated with the global one $\phi^{i}$, namely $[\phi^{i}]_{k}\doteq\phi^{a,k}$ with $i=g^{a,k}$. Note that $[\phi^{i}]_{k}$ and $\phi^{i}$ coincide inside the cell $T^{k}$, but they are different elsewhere.

Finally, we introduce the following classification of the global DOF ids into interior and exterior DOF ids. The set of interior DOF ids is defined as $\mathcal{I}^{\rm in}\doteq\{i\in\mathcal{I}:\ i\in\mathcal{G}(k)\text{ for some }k\in\mathcal{K}^{\rm in}\}$, i.e., all global DOF ids touched by the interior cells. On the other hand, the set of exterior DOF ids is defined as $\mathcal{I}^{\rm out}=\mathcal{I}\setminus\mathcal{I}^{\rm in}$. Clearly, the ids in $\mathcal{I}^{\rm out}$ are only touched by cut cells. In Fig. 3(b), interior and exterior DOF ids are represented with blue dots and red crosses, resp. Let us also define a map $K^{\rm own}:\mathcal{I}^{\rm out}\rightarrow\mathcal{K}$ that takes a DOF id $i\in\mathcal{I}^{\rm out}$ and returns the cell id of one of the (possibly) several cells that contain the DOF id $i$ (see Fig.  3(c)). In our implementation, we define $K^{\rm own}(i)$ as the smaller cell id $k\in\mathcal{K}$ such that $i\in\mathcal{G}(k)$. By composing the root cell map $R$ with $K^{\rm own}$, we obtain another map $O:\mathcal{I}^{\rm out}\rightarrow\mathcal{K}^{\rm in}$ that, for each exterior DOF id $i\in\mathcal{I}^{\rm out}$, returns an interior cell id. This latter map is defined as $O(i)\doteq R(K^{\rm own}(i))$. We denote cell $T^{O(i)}$ as the owner cell of the DOF id $i$, or the root cell associated with DOF id $i$, see Fig.  3(c).

With these notations, the agFE space reads

[TABLE]

where $C_{ij}\doteq[\phi^{j}]_{O(i)}(x^{i})$ and $\mathcal{M}(i)\doteq\mathcal{G}(O(i))$. Note that the space $\mathcal{V}^{\rm agg}$ is defined taking as starting point $\mathcal{V}^{\rm std}$, and adding some judiciously defined constraints to all exterior DOFs (see Fig. 3(a) and 3(b)). Only interior DOFs are free in the space $\mathcal{V}^{\rm agg}$. Thus, the dimension of $\mathcal{V}^{\rm agg}$ is equal to $|\mathcal{I}^{\rm in}|$. Moreover, $\mathcal{V}^{\rm agg}\subseteq\mathcal{V}^{\rm std}$ by construction. Each exterior DOF $i\in\mathcal{I}^{\rm out}$ is constrained as a linear combination of the DOF ids in the set $\mathcal{M}(i)$. Consequently, $\mathcal{M}(i)$ denotes the master DOF ids of a given exterior DOF id $i\in\mathcal{I}^{\rm out}$. Note that $\mathcal{M}(i)$ contains all global DOF ids in the owner cell of $i$. Thus, constraints in (2) state that the value of $v$ at an outer node with id $i\in\mathcal{I}^{\rm out}$ is computed as the extrapolation of the values of $v$ at the nodes of the owner cell of $i$ (e.g., in Fig. 3(c), the value at the red cross is computed as an extrapolation of the values at the blue dots).

2.4 Shape functions

Before starting the discussion on the distributed-memory implementation, we detail how to assemble the system of linear algebraic equations associated with the agFE space $\mathcal{V}^{\rm agg}$ in the serial case. To this end, we first need to introduce a shape function basis of $\mathcal{V}^{\rm agg}$, which is the purpose of this section. Let us define a restriction operator that, given a function $v\in\mathcal{V}^{\rm std}$, returns a function $\tilde{v}\in\mathcal{V}^{\rm agg}$ defined as:

[TABLE]

Note that function $\tilde{v}$ is computed in terms of the interior DOF ids only and it fulfills the constraints in (2) by construction. Moreover, $v$ and $\tilde{v}$ are different only in cut cells. The restriction does not affect the value of $v$ in interior cells. A shape function basis $\{\tilde{\phi}^{g}\}_{g\in\mathcal{I}^{\rm in}}$ of $\mathcal{V}^{\rm agg}$ is obtained by applying restriction (3) to the shape functions of $\mathcal{V}^{\rm std}$ associated with interior DOFs. More precisely, the shape function of $\mathcal{V}^{\rm agg}$ associated with the DOF id $g\in\mathcal{I}^{\rm in}$ is defined as

[TABLE]

Since the global shape functions of $\mathcal{V}^{\rm std}$ are such that $\phi^{i}(x^{j})=\delta_{ij}$ for any $i,j\in\mathcal{I}$, previous formula simplifies to

[TABLE]

where we have introduced the set $\mathcal{Z}(g)\doteq\{i\in\mathcal{I}^{\rm out}:\ g\in\mathcal{M}(i)\}$ referred to as the * slave DOF ids* of an interior DOF id $g\in\mathcal{I}^{\rm in}$. Note that the functions $\tilde{\phi}^{i}$ satisfy $\tilde{\phi}^{i}(x^{j})=\delta_{ij}$ for any $i,j\in\mathcal{I}^{\rm in}$. Thus, $\{\tilde{\phi}^{i}\}_{i\in\mathcal{I}^{\rm in}}$ is the shape function basis of $\mathcal{V}^{\rm agg}$ associated with the set of nodes $\{x^{j}\}_{j\in\mathcal{I}^{\rm in}}$. Note that the shape functions of the agFE are associated with the interior nodes of the mesh.

2.5 Finite element assembly

We present now the assembly operations needed to build the matrix $\mathbf{A}$ and vector $\mathbf{b}$ representing the discrete counterparts of a given bi-linear form $\mathrm{a}(\cdot,\cdot)$ and a linear form $\mathrm{b}(\cdot)$ with respect to the agFE space $\mathcal{V}^{\rm agg}$. The strategy we have followed is related with the way linear constraints are imposed in the assembly process in reference [44].

Using the previously defined shape functions of $\mathcal{V}^{\rm agg}$, the entries of the matrix $\mathbf{A}$ and the vector $\mathbf{b}$ are defined as $\left[\mathbf{A}\right]_{ij}\doteq\mathrm{a}(\tilde{\phi}^{i},\tilde{\phi}^{j})$ and $\left[\mathbf{b}\right]_{i}\doteq\mathrm{b}(\tilde{\phi}^{i})$ for each $i,j\in\mathcal{I}^{\rm in}$. As usual, the global matrix and vector are computed as the assembly of elemental matrices and vectors $[\mathbf{A}_{k}]_{ab}\doteq\mathrm{a}_{k}(\phi^{a,k},\phi^{b,k})$ and $[\mathbf{b}_{k}]_{a}\doteq\mathrm{b}_{k}(\phi^{a,k})$, $a,b\in\Lambda(T^{k})$, associated with a generic cell $T^{k}\in\mathcal{T}$. Here, $\mathrm{a}_{k}(\cdot,\cdot)$ (resp. $\mathrm{b}_{k}(\cdot)$) is the restriction of $\mathrm{a}(\cdot,\cdot)$ (resp. $\mathrm{b}(\cdot)$) to cell $T^{k}\in\mathcal{T}$.

For local DOF ids $a\in\Lambda(T^{k})$ such that its corresponding global DOF id $g^{a,k}\in\mathcal{I}^{\rm in}$ is interior, the assembly is standard, i.e., the value $[\mathbf{b}_{k}]_{a}$ is assembled into the global entry $[\mathbf{b}]_{i}$ such that $i=g^{a,k}$. However, for local ids such that $g^{a,k}\notin\mathcal{I}^{\rm in}$, the assembly is slightly different. To clarify this, we use the definition of $\tilde{\phi}^{g}$ (see Eq. (5)) in order to get a more detailed expression for $\left[\mathbf{A}\right]_{ij}$ and $\left[\mathbf{b}\right]_{i}$, namely

[TABLE]

[TABLE]

for $i,j\in\mathcal{I}^{\rm in}$. In addition, we rewrite formula (6) in terms of cell-local quantities:

[TABLE]

with $i=g^{a,k}$, $n=g^{b,k}$, and $\mathcal{B}(i)\doteq\{k\in\mathcal{K}:i=g^{c,k}\text{ for some }c\in\Lambda(T^{k})\}$ being the set containing the ids of cells sharing DOF id $i$. It is clear by observing formula (8) that the assembly of the entry $[\mathbf{b}_{k}]_{a}$ associated with a local DOF id $a\in\Lambda(T^{k})$ such that $g^{a,k}\in\mathcal{I}^{\rm in}$ is done as usual, i.e., $[\mathbf{b}_{k}]_{a}$ is added to the corresponding global entry $[\mathbf{b}]_{i}$ with $i=g^{a,k}$. However, the assembly is different for a local DOF id $a\in\Lambda(T^{k})$ such that $g^{a,k}\notin\mathcal{I}^{\rm in}$. By observing the last term in Eq. (8), the local entry $[\mathbf{b}_{k}]_{b}$ contributes to the global entry $\left[\mathbf{b}\right]_{i}$ for all global DOF ids $i\in\mathcal{I}^{\rm in}$ that satisfy $g^{b,k}\in\mathcal{Z}(i)$, which is equivalent to say that $i\in\mathcal{M}(g^{b,k})$. Thus, the value $[\mathbf{b}_{k}]_{b}$ scaled by a factor $C_{ni}$, with $n=g^{b,k}$, is added to all global entries $\left[\mathbf{b}\right]_{i}$ that satisfy $i\in\mathcal{M}(g^{b,k})$. That is, the local entries $[\mathbf{b}_{k}]_{a}$ associated with constrained DOF ids $g^{a,k}\notin\mathcal{I}^{\rm in}$ contribute to the global entries corresponding to the master DOF ids in $\mathcal{M}(g^{a,k})$. This is a convenient way of imposing constraints, since it does not lead to saddle point problems as it would be the case if we would use the method of Lagrange multipliers.

The same rationale applies for the entries of the system matrix. In this case, the assembly operation to be performed depends on the two local DOF ids associated with $\left[\mathbf{A}_{k}\right]_{ab}$, leading to four different cases. A summary of all different assembly operations is given in Table 1 for all possible cases. The symbol += stands for an in-place sum, i.e., the value of the right is added into the value on the left. Note that the procedures needed to assemble in $\mathcal{V}^{\rm agg}$ are closely related to the ones required in octree-based adaptive mesh refinement codes with non-conforming meshes. The imposition of continuity constraints at so called hanging nodes (see, e.g., [45]) leads to assembly operations equivalent to the ones presented in Table 1. Thus, agFE spaces can be easily incorporated in existing octree-based adaptive codes.

3 Distributed-memory implementation

3.1 Domain decomposition setup

After the presentation of the AgFEM in a serial context, we discuss its implementation in a distributed-memory framework. We start by introducing some notation associated with the partition of the mesh into sub-domains, which is the basis of the distributed-memory implementation. Let $\mathcal{D}$ be a partition of $\Omega^{\rm art}$ into sub-domains obtained by the union of cells in the background mesh $\mathcal{T}^{\rm art}$, i.e., for each cell $T\in\mathcal{T}^{\rm art}$, there is a sub-domain $D\in\mathcal{D}$ such that $T\subset D$ (see Fig. 4(a)). The generation of the sub-domain partition $\mathcal{D}$ is specially simple in the context of embedded boundary methods since Cartesian grids are allowed. Here, we consider the load-balancing mechanism available in p4est, which is based on space-filling curves [46] and has shown excellent scalability up to hundreds of thousands of cores. Note that our approach does not require graph-based partitioning techniques such as metis or parmetis [47].

For convenience, we consider a numbering of the sub-domains in $\mathcal{D}$, where $D_{s}$ is the sub-domain with id $s$ and $\mathcal{S}\doteq\{1,\ldots,|\mathcal{D}|\}$ is the index set containing all sub-domain ids. With this notation the set of sub-domains $\mathcal{D}$ is expressed as $\mathcal{D}=\{D_{s}\}_{s\in\mathcal{S}}$. For a given $s\in\mathcal{S}$, we define the set of local cells $\mathcal{T}^{\rm L}_{s}$ as the cells included in $D_{s}$, i.e., $\mathcal{T}^{\rm L}_{s}\doteq\{T\in\mathcal{T}^{\rm art}:\ T\subset D_{s}\}$. On the other hand, we define the set of ghost cells $\mathcal{T}^{\rm G}_{s}$ of $D_{s}$ as the off-subdomain cells that share at least one vertex, edge or face with $D_{s}$, namely $\mathcal{T}^{\rm G}_{s}\doteq\{T\in(\mathcal{T}^{\rm art}\setminus\mathcal{T}^{\rm L}_{s}):\ {}\mkern 1.0mu\overline{\mkern-1.0muT}\cap{}\mkern 1.0mu\overline{\mkern-1.0muD_{s}}\neq\emptyset\}$ (see Fig. 4(b), and 4(c)). The union of local and ghost cells, $\mathcal{T}_{s}\doteq\mathcal{T}^{\rm L}_{s}\cup\mathcal{T}^{\rm G}_{s}$, is called the set of locally relevant cells in sub-domain $D_{s}$.

For the following developments, we need a local numbering for the locally relevant cells $\mathcal{T}_{s}$. That is, we write $\mathcal{T}_{s}=\{T^{\ell}_{s}\}_{\ell\in\mathcal{L}_{s}}$, where $T^{\ell}_{s}$ denotes the cell with local id $\ell$ in the scope of sub-domain $D_{s}$ and $\mathcal{L}_{s}\doteq\{1,\ldots,|\mathcal{T}_{s}|\}$ is the set of all local cell ids associated with $D_{s}$. In addition, we define $\mathcal{L}^{\rm L}_{s}\doteq\{\ell\in\mathcal{L}_{s}:\ T^{\ell}_{s}\in\mathcal{T}^{\rm L}_{s}\}$ and $\mathcal{L}^{\rm G}_{s}\doteq\{\ell\in\mathcal{L}_{s}:\ T^{\ell}_{s}\in\mathcal{T}^{\rm G}_{s}\}$, which contain the local cell ids of local and ghost cells resp. We also need to identify which is the sub-domain of each cell in terms of local cell ids. To this end, we denote as $S_{s}(\ell)$ the id of the sub-domain, where the cell with id $\ell\in\mathcal{L}_{s}$ is a local cell. Clearly, $S_{s}(\ell)=s$ for all local cell ids $\ell\in\mathcal{L}^{\rm L}_{s}$ and $S_{s}(\ell)\neq s$ for ghost cell ids $\ell\in\mathcal{L}^{\rm G}_{s}$.

For simplicity and without any loss of generality, we will assume (as we have done in the serial case) that exterior cells have been removed since they do not play any role. In particular, it implies that $\mathcal{T}^{\rm L}_{s}$, $\mathcal{T}^{\rm G}_{s}$, and $\mathcal{T}_{s}$ contain only interior or cut cells. As already stated in Sect.  2.1, each cell in the active mesh $\mathcal{T}$ receives a global id $k\in\mathcal{K}$, which is unique across all cells in $\mathcal{T}$. Since $\mathcal{T}_{s}$ is a subset of $\mathcal{T}$, all cells in $\mathcal{T}_{s}$ have a global cell id. We denote by $K_{s}(\ell)$ the global cell id associated with a local cell id $\ell\in\mathcal{L}_{s}$. Conversely, we denote as $L_{s}(k)$ the local cell id associated with the global cell id $k\in\mathcal{K}_{s}$. Here, $\mathcal{K}_{s}\doteq\{k\in\mathcal{K}:\ T^{k}\in\mathcal{T}_{s}\}$ is the set of all global ids assigned to cells in $\mathcal{T}_{s}$. With this notation, we can write $\mathcal{T}_{s}=\{T^{k}\}_{k\in\mathcal{K}_{s}}$. Finally, we need to introduce $\mathcal{S}^{\rm nei}_{s}\doteq\{s^{\prime}\in\mathcal{S}\setminus\{s\}:\ {}\mkern 1.0mu\overline{\mkern-1.0muD_{s^{\prime}}}\cap{}\mkern 1.0mu\overline{\mkern-1.0muD_{s}}\neq\emptyset\}$, i.e., the set containing the sub-domain ids of the nearest neighbors of $D_{s}$.

3.2 Distributed finite element spaces

The goal of this section is to express the agFE space $\mathcal{V}^{\rm agg}$ in terms of sub-domain local quantities in the framework of the domain decomposition setup. By doing so, most of the notation needed in the subsequent text will be introduced in a natural way. For the sake of simplicity, we assume from now on (if not otherwise indicated) that we are located in the local scope of a generic sub-domain $D_{s}$ with id $s\in\mathcal{S}$. Also for simplicity (and without any loss of generality), we assume that we deal with scalar-valued Lagrangian FE interpolations.

First, we define $\mathcal{V}^{\rm std}_{s}\doteq\mathcal{V}^{\rm std}|_{D_{s}}$ the restriction of the standard FE space $\mathcal{V}^{\rm std}$ to the sub-domain $D_{s}$. Given a FE function $v\in\mathcal{V}^{\rm std}$, we also introduce its sub-domain restriction: $v_{s}\doteq v|_{D_{s}}$. Note that $v_{s}$ is characterized by a subset of the DOFs of $v$, which are the ones associated with the nodes of the mesh $\mathcal{T}^{\rm L}_{s}$, i.e., the nodes in the set $\mathcal{X}_{s}\doteq\{x\in\mathcal{X}:\ x\in{}\mkern 1.0mu\overline{\mkern-1.0muT}\text{ for some }T\in\mathcal{T}^{\rm L}_{s}\}$. As we have done previously for cells, we introduce a sub-domain local numbering for nodes in $\mathcal{X}_{s}$, i.e., $\mathcal{X}_{s}=\{x^{j}_{s}\}_{j\in\mathcal{J}_{s}}$, where $x^{j}_{s}$ denotes the nodal coordinate of the mesh node with local id $j$ in sub-domain $D_{s}$ and $\mathcal{J}_{s}\doteq\{1,\ldots,|\mathcal{X}_{s}|\}$ is the set of all sub-domain local node ids (also referred as set of local DOF ids since there is a one-to-one correspondence between nodes and DOFs). Thus, each node in $\mathcal{X}_{s}$ has two different ids, a local one in the local mesh $\mathcal{T}_{s}$ and a global one in the entire mesh $\mathcal{T}$. Let $\mathcal{I}_{s}\doteq\{i\in\mathcal{I}:\ x^{i}\in\mathcal{X}_{s}\}$ be the set of the global node ids (or global DOF ids) in sub-domain $D_{s}$ and let $I_{s}:\mathcal{J}_{s}\rightarrow\mathcal{I}_{s}$ be the map that, for a given local node id $j$, returns the corresponding global id $i=I_{s}(j)$ such that $x^{i}=x^{j}_{s}$. With these notations, we can write a generic function $v_{s}\in\mathcal{V}^{\rm std}_{s}$ as a linear combination of sub-domain local quantities, i.e., $v_{s}(x)=\sum_{j\in\mathcal{J}_{s}}v^{j}_{s}\ \phi^{j}_{s}(x)$, where $v^{j}_{s}\doteq v_{s}(x^{j}_{s})=v^{i}=v(x^{i})$ and $\phi^{j}_{s}\doteq\phi^{i}|_{D_{s}}$ with $i=I_{s}(j)$ for all $j\in\mathcal{J}_{s}$. Clearly, $v^{j}_{s}$ and $\phi^{j}_{s}$ are the DOF and the shape function associated with the sub-domain local id $j$.

On the other hand, we introduce $\mathcal{V}^{\rm agg}_{s}\doteq\mathcal{V}^{\rm agg}|_{D_{s}}$ the restriction of the agFE space $\mathcal{V}^{\rm agg}$ to sub-domain $D_{s}$. In order to characterize $\mathcal{V}^{\rm agg}_{s}$ in terms of sub-domain local quantities, we need to introduce some further notation. First, let us define $\phi^{a}_{s,\ell}$ the shape function associated with a generic DOF id $a\in\Lambda(T^{\ell}_{s})$ in cell $T^{\ell}_{s}$, namely $\phi^{a}_{s,\ell}\doteq\phi^{a,k}$ with $k=K_{s}(\ell)$. Analogously, we introduce $x^{a}_{s,\ell}\doteq x^{a,k}$ and $g^{a}_{s,\ell}\doteq g^{a,k}$, with $k=K_{s}(\ell)$, which represent the nodal coordinates and the global DOF id associated with a generic DOF id $a\in\Lambda(T^{\ell}_{s})$ in cell $T^{\ell}_{s}$. Note that the symbols $g^{a}_{s,\ell}$ and $g^{a,k}$ express the same quantities, but the former expresses them in terms of local cell ids $\ell\in\mathcal{L}_{s}$, and the latter in terms of global cell ids $k\in\mathcal{K}_{s}$. This is true also for $x^{a}_{s,\ell}$ and $x^{a,k}$, and $\phi^{a}_{s,\ell}$ and $\phi^{a,k}$, resp. We also define $j^{a}_{s,\ell}$ as the sub-domain local DOF id in $\mathcal{J}_{s}$ associated with the cell local DOF id $a\in\Lambda(T^{\ell}_{s})$ in cell $T^{\ell}_{s}$. Formally, $j^{a}_{s,\ell}$ is expressed as $j^{a}_{s,\ell}\doteq I^{-1}_{s}(g^{a}_{s,\ell})$ in terms of the inverse application $I^{-1}_{s}$ (which is not needed to be generated in practice). The process that generates $g^{a}_{s,\ell}$ and $j^{a}_{s,\ell}$ for each $a$, $\ell$, and $s$ is standard in FE packages. The sub-domain local cell-wise DOF ids $j^{a}_{s,\ell}$ are simply generated by applying independently in each sub-domain the standard procedure for generating DOF ids in the serial case (see, e.g., [35] for details). On the other hand, the generation of the global DOF ids $g^{a}_{s,\ell}$ is also done by means of standard tools (see, e.g., [48, 45]) as it is detailed later at the end of Sect. 3.6.

We also need to introduce a classification of the locally relevant cells into interior, exterior, and cut. To this end, we define the index sets $\mathcal{L}_{s}^{\rm in}\doteq\{\ell\in\mathcal{L}_{s}:\ T^{\ell}_{s}\in\mathcal{T}^{\rm in}\}$ (resp. for $\mathcal{L}^{\rm out}_{s}$ and $\mathcal{L}^{\rm cut}_{s}$) containing the subset of $\mathcal{L}_{s}$ corresponding to interior cells (resp. exterior and cut). In our case, the set $\mathcal{L}^{\rm out}_{s}$ is empty since we have assumed that exterior cells have been removed from the mesh. The classification of cells as interior or cut induces the following partition of sub-domain local DOF ids in $\mathcal{J}_{s}$. Let $\mathcal{J}^{\rm in}_{s}\doteq\{j\in\mathcal{J}_{s}:\ j=j^{a}_{s,\ell}\text{ for some }a\in\Lambda(T^{\ell}_{s})\text{ and }\ell\in\mathcal{L}^{\rm in}_{s}\}$ be the set of interior sub-domain local DOF ids, and let $\mathcal{J}^{\rm out}_{s}\doteq\mathcal{J}_{s}\setminus\mathcal{J}^{\rm in}_{s}$ be the set of exterior sub-domain local DOF ids. Note that the DOF ids in $\mathcal{J}^{\rm in}_{s}$ are the ones touched by interior cells with ids in $\mathcal{L}^{\rm in}_{s}$, whereas the ids in $\mathcal{J}^{\rm out}_{s}$ are only touched by cut cells.

The last ingredient to characterize $\mathcal{V}^{\rm agg}_{s}$ in terms of sub-domain quantities is a distributed version of the root cell map $R$ and the map $K^{\mathrm{\rm own}}$. On the one hand, the distributed root cell map $R_{s}:\mathcal{L}_{s}\rightarrow\mathcal{K}^{\rm in}$ is a transformation that takes a local cell id $\ell\in\mathcal{L}_{s}$ and returns the global cell id $R_{s}(\ell)\in\mathcal{K}^{\rm in}$ of the root cell of $T^{\ell}_{s}$, namely $R_{s}(\ell)\doteq R(K_{s}(\ell))$. We encode the input value of $R_{s}(\cdot)$ with a local cell id in $\mathcal{L}_{s}$ since we only need to define $R_{s}$ for locally relevant cells in $\mathcal{T}_{s}$. However, the returned value of $R_{s}$ has to be encoded with a global cell id in $\mathcal{K}^{\rm in}$, and not with a local cell id in $\mathcal{L}_{s}$, since the root cell of $T^{\ell}_{s}$ can be outside the set of locally relevant cells $\mathcal{T}_{s}$ (which are the only ones that have a local id in sub-domain $D_{s}$). Note also that, even though we have formally defined $R_{s}$ using $R$, the map $R$ is never built in practice. Constructing and storing $R$ so that each processor can retrieve the value $R(k)$ for any global cell id $k\in\mathcal{K}$ in the entire mesh would kill parallel (memory) scalability. Instead, we only build the distributed map $R_{s}$ as shown later in Sect. 3.4. In addition to $R_{s}$, we also construct another map $S^{\rm root}_{s}:\mathcal{L}_{s}\rightarrow\mathcal{S}$ that for each local cell id $\ell\in\mathcal{L}_{s}$ returns the id $S^{\rm root}_{s}(\ell)\in\mathcal{S}$ of the sub-domain in which the corresponding root cell is a local cell, namely $s^{\prime}=S^{\rm root}_{s}(\ell)$ is such that $T^{\ell^{\prime}}_{s^{\prime}}\in\mathcal{T}^{\rm L}_{s^{\prime}}$ with $\ell^{\prime}=L_{s^{\prime}}(R_{s}(\ell))$. Finally, we need to build a transformation $L^{\rm own}_{s}:\mathcal{J}^{\rm out}_{s}\rightarrow\mathcal{L}_{s}$, which is a distributed version of $K^{\mathrm{\rm own}}$. For a local DOF id $j\in\mathcal{J}^{\rm out}_{s}$, the value $L^{\rm own}_{s}(j)$ is the local cell id $\ell\in\mathcal{L}_{s}$ of one of the (possibly several) cells $T^{\ell}_{s}$ such that $j=j^{a}_{s,\ell}$ for some $a\in\Lambda(T^{\ell}_{s})$. That is, $L^{\rm own}_{s}(j)$ is the local id of one of the (possibly) several cells that contain the sub-domain DOF id $j$. When more than one cell can be selected, we choose one arbitrarily, e.g., the cell that has smaller global id $k\in\mathcal{K}_{s}$. Note that one can choose one cell arbitrarily, but the chosen criterion has to be consistent across all sub-domains. E.g., choosing the cell with smaller global cell id is consistent, but choosing the cell with smaller local id is not.

With these notations, we can finally write $\mathcal{V}^{\rm agg}_{s}$ in terms of sub-domain quantities:

[TABLE]

with $C^{ja}_{s}\doteq\phi^{a}_{s^{\prime},\ell^{\prime}}(x^{j}_{s})$, $\ell^{\prime}=L_{s^{\prime}}(R_{s}(L^{\rm own}_{s}(j)))$ and $s^{\prime}=S^{\rm root}_{s}(L^{\rm own}_{s}(j))$. By observing Eq. (9), the main difficulty of the implementation of the method becomes apparent. Note that in the scope of the current sub-domain $D_{s}$, not all quantities are available in order to define the constraints in Eq. (9). In order to evaluate the coefficient $C^{ja}_{s}$, one needs access to the shape function in physical space $\phi^{a}_{s^{\prime},\ell^{\prime}}$, which is one of the local shape basis functions of the cell with local id $\ell^{\prime}=L_{s^{\prime}}(R_{s}(L^{\rm own}_{s}(j)))$ in sub-domain $s^{\prime}=S^{\rm root}_{s}(L^{\rm own}_{s}(j))$, which can be interpreted as the root cell associated with the DOF id $j$. Building up $\phi^{a}_{s^{\prime},\ell^{\prime}}$ requires the physical coordinates of the nodes of the root cell, namely $x^{b}_{s^{\prime},\ell^{\prime}}$ for $b\in\Lambda(T^{\ell^{\prime}}_{s^{\prime}})$. Since, in general, the sub-domain where the root cell is located is different from the current sub-domain, i.e., $s^{\prime}\neq s$, the quantities $x^{b}_{s^{\prime},\ell^{\prime}}$ are not available in the scope of $D_{s}$ (see Fig. 5). Moreover, one cannot directly use the communication pattern used to exchange cell-based data between nearest neighbor sub-domains because id $s^{\prime}$ is not even necessarily the id of a nearest neighbor sub-domain, namely $s^{\prime}\notin\mathcal{S}^{\rm nei}_{s}$ in general. Thus, an ad-hoc communication strategy has to be built in order to import quantities like $x^{b}_{s^{\prime},\ell^{\prime}}$ into the scope of the current processor.

3.3 Implementation overview

In the following sections (Sects. 3.4–3.7), we discuss the implementation of the main algorithmic phases of the distributed-memory AgFEM method. We start by presenting the distributed version of the cell aggregation method presented in Alg. 1 (see Sect. 3.4), where we show the actual computation of the distributed root cell map $R_{s}$ and the auxiliary maps $K^{\rm own}_{s}$ and $S^{\rm root}_{s}$. Then, in Sect. 3.5, we detail how to build the communication strategy in order to import data from root cells that are located in remote sub-domains. Finally, we detail the data distribution model used for representing the linear system of equations in the distributed-memory implementation (Sect. 3.6) and the parallel FE assembly associated with the agFE space (Sect. 3.7).

3.4 Parallel cell aggregation

The goal of this section is to parallelize the cell aggregation method presented in Alg. 1. To this end, the most naive approach would be to run Alg. 1 on the local cells $\mathcal{T}^{\rm L}_{s}$ in each processor independently, trying to greedily associate each cut cell with a root cell in the same sub-domain. However, this approach does not work when a sub-domain has cut cells that cannot be aggregated to any local cell (i.e., the cut cell does not share a face with any local internal cell), and therefore, it is not general enough. Moreover, even when the method is algorithmically applicable (e.g., for domain decomposition setups where each sub-domain has at least one interior cell), it does not necessarily lead to the same aggregates as the serial method in Alg. 1. The mathematical analysis of the AgFEM method in [10] states that the approximation properties of the agFE spaces are negatively affected as the aggregate size increases. Thus, this naive parallel strategy that can potentially lead to large aggregates has to be avoided.

Instead, we consider a more sophisticated parallel algorithm, which leads to the same aggregates as the ones obtained with the serial aggregation method. This enhanced distributed-memory cell aggregation method is detailed in Alg. 2 and illustrated in Fig. 6 for a simple two dimensional case. The notation $\mathcal{L}^{\rm nei}_{s}(\ell)$, in line 2 of the algorithm represents the set of local cell ids of the face neighbors of cell $T^{\ell}_{s}$. That is, the notations $\mathcal{K}^{\rm neig}(k)$ (see Sect. 2.2) and $\mathcal{L}^{\rm nei}_{s}(\ell)$ both represent the face neighbors of a given cell, but the former in terms of global cell ids and the latter in terms of local cell ids. The parallel cell aggregation method in Alg. 2 follows essentially the same steps as the serial counterpart plus some communications between neighbor sub-domains, which allow one to aggregate a cell into a root belonging to a different sub-domain. The method consists in three main steps (see Fig. 6). (a) We initialize the value of $R_{s}$ for all interior local cells and perform a nearest neighbor communication in order to determine the value of $R_{s}$ at ghost cells. At this point, we mark all interior (local and ghost) cells as "touched". (b) As in the serial case, we perform a loop in (local) cut cells. For a given cut cell, we define its root cell taking the root cell of one of its face neighbors that have been already marked as "touched" (if such a neighbor exists). (c) After the loop, we perform a nearest neighbor communication in order to update the value of $R_{s}$ at ghost cells. At this point, we mark as "touched" all (local and ghost) cells for with $R_{s}$ is already set with a valid root cell id. Finally, we repeat steps (b) and (c) in all processors as many times as needed until all the cells are assigned with a root cell id. In the numerical experiments, 3 repetitions of steps (b) and (c) were needed to aggregate all cut cells for all cases.

Note that, during the execution of Alg. 2, we also build the map $S^{\rm root}_{s}:\mathcal{L}_{s}\rightarrow\mathcal{S}$, which was formally presented in Sect. 3.2. In addition, we build an auxiliary map $N_{s}:\mathcal{L}_{s}\rightarrow\mathcal{K}$, which is the distributed version of the map $N:\mathcal{K}\rightarrow\mathcal{K}$ presented in Sect. 2.2. Precisely, $N_{s}(\ell)$ returns the global id of the cell that was selected in line 2 of the algorithm, i.e., the next cell in the path towards the root cell. This information will be required later in the next section (Sect. 3.5) in order to build the communication strategy for importing data from remote root cells.

3.5 Importing data from root cells

As we have discussed in Sect. 3.2, we need to import data associated with remote root cells in order to build the AgFEM constraints in Eq. (9). Formally, the set containing the global cell ids of those root cells is defined as $\mathcal{K}^{\rm R}_{s}\doteq\{k\in\mathcal{K}\setminus\mathcal{K}_{s}:k=R_{s}(\ell)\text{ for some }\ell\in\mathcal{L}_{s}\}$. For assembling the system matrix and right hand side vector, we need to import the cell-wise nodal coordinates and global DOF ids for all cells with global ids in $\mathcal{K}^{\rm R}_{s}$ as it will be justified in Sect. 3.7. This information is fetched from the corresponding remote sub-domains (as we discuss later) and stored in the local scope of the processor in a separate buffer allocation. In order to address data in this buffer, each cell with global id $k$ in the set $\mathcal{K}^{\rm R}_{s}$ is labeled with a new sub-domain local id $z$, defined as the position of $k$ in the list that represents the index set $\mathcal{K}^{\rm R}_{s}$. This transformation is represented here by a map $Z_{s}:\mathcal{K}^{\rm R}_{s}\longrightarrow\mathcal{L}^{\rm R}_{s}\doteq\{1,\ldots,|\mathcal{K}^{\rm R}_{s}|\}$ such that, given a cell global id $k\in\mathcal{K}^{\rm R}_{s}$, it returns the corresponding $z$ id. We denote by $\tilde{T}^{z}_{s}$ the cell associated with label $z\in\mathcal{L}^{\rm R}_{s}$, namely $\tilde{T}^{z}_{s}\doteq T^{k}$ with $z=Z_{s}(k)$. In addition, we denote as $\tilde{x}^{a}_{s,z}\doteq x^{a,k}$, $z=Z_{s}(k)$, the coordinates of the $a$-th node of cell $\tilde{T}^{z}_{s}$ for each $a\in\Lambda(\tilde{T}^{z}_{s})$. We use an analogous notation for $\tilde{g}^{a}_{s,z}\doteq g^{a,k}$, $\tilde{\phi}^{a}_{s,z}\doteq\phi^{a,k}$, and $\tilde{v}^{a}_{s,z}\doteq v^{a,k}$ (i.e., cell-wise global DOF ids, cell shape functions, and cell DOF values of a given FE function $v\in\mathcal{V}^{\rm agg}$).

Note that the cells with global id in $\mathcal{K}^{\rm R}_{s}$ are local cells in the sub-domains with id in the set $\mathcal{Q}^{\rm rcv}_{s}\doteq\{s^{\prime}\in\mathcal{S}:\ s^{\prime}=S^{\rm root}_{s}(\ell)\text{ for some }\ell\in\mathcal{L}^{\rm cut}_{s}\}$. That is, we import data from all sub-domains that contain the root cells of the cut cells in the current sub-domain. We note that the root cells of ghost cut cells also need to be imported as $L_{s}^{\rm own}(j)$ may be in $\mathcal{L}_{s}^{\rm G}$ for some $j\in\mathcal{J}^{\rm out}_{s}$. We also note that we have constructed the map $S^{\rm root}_{s}$ during the parallel cell aggregation algorithm in order to be able to construct the set $\mathcal{Q}^{\rm rcv}_{s}$. It is also easy to see that the set $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}\doteq\{k\in\mathcal{K}_{s^{\prime}}:\ k=R_{s}(\ell)\text{ for some }\ell\in\mathcal{L}^{\rm cut}_{s}\ \}$ contains, for a given $s^{\prime}\in\mathcal{Q}^{\rm rcv}_{s}$, the global ids of root cells in sub-domain $D_{s^{\prime}}$ from which we (potentially) need to recover data.

In order to grasp the sets $\mathcal{Q}^{\rm rcv}_{s}$ and $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$, we can consider the cell path that goes from each cut cell to its corresponding root cell, jumping along the face neighbors that are selected during the cell aggregation algorithm (see the blue arrows in Fig 6(d)). The set $\mathcal{Q}^{\rm rcv}_{s}$ contains the ids of the sub-domains, where the cell paths starting in the current sub-domain $D_{s}$ end. On the other hand, the set $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$ contains the ids of the last cell in each path that starts at $D_{s}$ and ends at $D_{s^{\prime}}$. The practical construction of the sets $\mathcal{Q}^{\rm rcv}_{s}$ and $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$ from the maps $R_{s}$ and $S^{\rm root}_{s}$ is indicated in Alg. 3. Note that this procedure is done at each sub-domain independently and does not require communication at all.

Once the sets $\mathcal{Q}^{\rm rcv}_{s}$ and $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$ are built, we know from which sub-domains, and from which cells in those sub-domains, we need to recover data. However, we still need to know to which sub-domains we need to send data from the current sub-domain, and which is the data to be sent to each of those sub-domains. The "send" counterparts of the sets $\mathcal{Q}^{\rm rcv}_{s}$ and $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$ are formally defined as $\mathcal{Q}^{\rm snd}_{s}\doteq\{s^{\prime}\in\mathcal{S}:\ s=S^{\rm root}_{s^{\prime}}(\ell^{\prime})\text{ for some }\ell^{\prime}\in\mathcal{L}^{\rm cut}_{s^{\prime}}\}$ and, for a given $s^{\prime}\in\mathcal{Q}^{\rm snd}_{s}$, $\mathcal{R}^{\rm snd}_{s,s^{\prime}}\doteq\{k\in\mathcal{K}_{s}:\ k=R_{s^{\prime}}(\ell^{\prime})\text{ for some }\ell^{\prime}\in\mathcal{L}^{\rm cut}_{s^{\prime}}\}$. We can interpret the set $\mathcal{Q}^{\rm snd}_{s}$ as the ids of the sub-domains, where all the cell paths ending at $D_{s}$ start. On the other hand, the set $\mathcal{R}^{\rm snd}_{s,s^{\prime}}$ contains the ids of the last cell in each path that ends at $D_{s}$ and starts at $D_{s^{\prime}}$. In contrast to the sets $\mathcal{Q}^{\rm rcv}_{s}$ and $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$, it is clear that we cannot construct the sets $\mathcal{Q}^{\rm snd}_{s}$ and $\mathcal{R}^{\rm snd}_{s,s^{\prime}}$ without communications since their definitions are in terms of the maps $R_{s^{\prime}}$ and $S^{\rm root}_{s^{\prime}}$ seen from remote sub-domains $D_{s^{\prime}}$ and thus not available at the current one $D_{s}$. However, thanks to the construction of the auxiliary map $N_{s}:\mathcal{L}_{s}\rightarrow\mathcal{K}$ during the parallel cell aggregation algorithm, we can build the sets $\mathcal{Q}^{\rm snd}_{s}$ and $\mathcal{R}^{\rm snd}_{s,s^{\prime}}$ by using only standard communications between nearest neighbor sub-domains.

The detailed construction of $\mathcal{Q}^{\rm snd}_{s}$ and $\mathcal{R}^{\rm snd}_{s,s^{\prime}}$ is given in Alg. 4. In the algorithm, we create a tuple for each cell path that starts at the current sub-domain (i.e., one tuple for each local cut cell). These tuples have three items each (see line 4): the global id of the first cell in the path, the global id of the next cell in the path, and the current sub-domain id. Conceptually, we gradually move forward these tuples along their cell paths until they arrive to the final cell in the path, which can be possibly in another sub-domain. At the end of the process each sub-domain collects all the tuples corresponding to the cell paths that end in this sub-domain. With this information each sub-domain is able to reconstruct the sets $\mathcal{Q}^{\rm snd}_{s}$ and $\mathcal{R}^{\rm snd}_{s,s^{\prime}}$. Note that Alg. 4 only involves communications with nearest neighbor sub-domains.

Finally, we detail in Alg. 5 how we import the required data from remote sub-domains by using the sets $\mathcal{Q}^{\rm snd}_{s}$ and $\mathcal{R}^{\rm snd}_{s,s^{\prime}}$ and $\mathcal{Q}^{\rm rcv}_{s}$ and $\mathcal{R}^{\rm rcv}_{s,s^{\prime}}$.

3.6 Linear algebra data distribution layout

Before detailing the FE assembly, we need to specify the data model used for representing the matrices and vectors in a distributed-memory context. Here, we consider row-wise partitioned globally addressable matrices and vectors (as we detail in the next paragraph). This choice is motivated by the fact that the vast majority of parallel linear algebra packages (e.g., PETSc [32, 33], or TRILINOS [31]) use this data distribution model. Other data layouts are also possible, such as sub-assembled distributed matrices for non-overlapping domain decomposition solvers [49, 23], but they are out of the scope of the current work.

In a row-wise partitioning model, the global matrix and the global vector are partitioned by rows, i.e., each processor owns a consecutive, non-overlapping, range of global rows of $\mathbf{A}$ and $\mathbf{b}$, respectively. Since each row of $\mathbf{A}$ and $\mathbf{b}$ is associated with a global DOF id in $\mathcal{I}^{\rm in}$, the consecutive row partition implies that each sub-domain $D_{s}$ owns a consecutive subset of the global DOF ids in $\mathcal{I}^{\rm in}$, denoted as $\mathcal{I}^{\rm in,O}_{s}$, where $\{\mathcal{I}^{\rm in,O}_{s}\}_{s\in\mathcal{S}}$ is a partition of $\mathcal{I}^{\rm in}$. Note that we do not need to partition the set of constrained DOF ids $\mathcal{I}^{\rm out}$ since, in our approach, these DOFs are removed from the global linear system. On the other hand, we consider linear algebra data structures that are globally addressable. That is, a given processor with sub-domain id $s\in\mathcal{S}$ may contribute to a global matrix row $[\mathbf{A}]_{gi}$, $i\in\mathcal{I}^{\rm in}$ (or a global vector entry $[\mathbf{b}]_{g}$), which is not necessarily owned by it, namely $g\in\mathcal{I}^{\rm in}\setminus\mathcal{I}^{\rm in,O}_{s}$. This is possible since, once the corresponding entries of the matrix and vector are set in the assembly process, a final communication stage transfers all (touched) values $[\mathbf{A}]_{gi}$, $[\mathbf{b}_{g}]$ such that $g\in\mathcal{I}^{\rm in}\setminus\mathcal{I}^{\rm in,O}_{s}$ to the corresponding remote processor that owns the global id $g$. This communication phase is typically resolved by the linear algebra library and, thus, we do not cover it here for simplicity.

Even though globally addressable linear algebra objects allow one to assign entries that are not owned by the current processor, it is desirable to minimize the number of such assignments in order to reduce inter processor communication. To this end, we define the set of owned global DOF ids $\mathcal{I}^{\rm in,O}_{s}$ following the usual approach in FE applications. Let $\mathcal{I}^{\rm in}_{s}\doteq\{i\in\mathcal{I}^{\rm in}:\ i=g^{a}_{s,\ell}\text{ for some }a\in\Lambda(T^{\ell}_{s})\text{ and }\ell\in\mathcal{L}^{\rm L}_{s}\}$ be the set of interior global DOF ids touched only by local cells of sub-domain $D_{s}$ and let $\mathcal{B}_{s}(i)\doteq\{\ell\in\mathcal{L}_{s}:\ i=g^{a}_{s,\ell}\text{ for some }a\in\Lambda(T^{\ell}_{s})\}$ be the ids of the cells around a given id $i\in\mathcal{I}^{\rm in}_{s}$. With this notation, we define the set of owned global DOF ids as $\mathcal{I}^{\rm in,O}_{s}\doteq\{i\in\mathcal{I}^{\rm in}_{s}:\ s\leq S_{s}(\ell)\text{ for all }\ell\in\mathcal{B}_{s}(i)\}$. Note that each processor owns all DOF ids which are surrounded solely by cells in $\mathcal{T}_{s}^{\rm L}$ plus some of the DOF ids that are on the interface between cells in $\mathcal{T}_{s}^{\rm L}$ and $\mathcal{T}_{s}^{\rm G}$. With this model, most of the cells in $\mathcal{T}^{\rm L}_{s}$ contribute to locally owned DOF ids in $\mathcal{I}^{\rm in,O}_{s}$ during the assembly process (i.e., most of the values set from sub-domain $D_{s}$ will correspond to values owned by $D_{s}$ if we assume that the number of DOFs in the interior of the sub-domain $D_{s}$ is larger than number of DOFs on its boundary). On the other hand, for DOF ids $i\in\mathcal{I}^{\rm in}_{s}$ that are located at the interface between cells in $\mathcal{T}_{s}^{\rm L}$ and $\mathcal{T}_{s}^{\rm G}$ (i.e., at sub-domain boundaries), we define the owner sub-domain as the sub-domain with smaller id among all sub-domain that share the DOF id $i$. We have used this tie-breaker strategy in order to decide which sub-domain is the owner for interface DOFs, but other criteria can be also considered as long as they are consistent across all processors and require no communication.

Even though it is not strictly required, it is convenient, in order to facilitate the interaction with the parallel linear algebra package, that the sets $\{\mathcal{I}^{\rm in,O}_{s}\}_{s\in\mathcal{S}}$ contain contiguous, sequentially increasing DOF ids, namely $\mathcal{I}^{\rm in,O}_{1}=\{1,\ldots,m^{\rm in}_{1}\}$, $\mathcal{I}^{\rm in,O}_{2}=\{m^{\rm in}_{1}+1,\ldots,m^{\rm in}_{1}+m^{\rm in}_{2}\}$, etc., where $m^{\rm in}_{s}\doteq|\mathcal{I}^{\rm in,O}_{s}|$ for $s\in\mathcal{S}$. Methods for the generation of cell-wise global DOF ids $g^{a}_{s,\ell}$ that result in sets $\{\mathcal{I}^{\rm in,O}_{s}\}_{s\in\mathcal{S}}$ fulfilling this property are well known in distributed-memory FE analysis (see, e.g., [48] for the algorithm implemented in deal.II and [45] for the algorithm implemented in FEMPAR). These techniques can readily be used to generate the cell-wise global DOF ids $g^{a}_{s,\ell}$ in the context of the AgFEM method by considering the following steps. First, the values $g^{a}_{s,\ell}$ are distributed over interior cell ids $\ell\in\mathcal{L}^{\rm in}_{s}$ by applying the conventional algorithms [48, 45] restricted to the interior portion of the mesh (i.e., as though cut cells were removed from the mesh). In a second phase, global DOF ids $g^{a}_{s,\ell}$ are assigned to the vertices, edges and faces of cut cells that are in contact with an interior cell (if any) by taking the corresponding values of the neighbor interior cell. Note that, at the end of this process, the cell-wise global DOF ids $g^{a}_{s,\ell}$ are not defined in the interior of cut cells nor on the parts of their boundaries that are not in contact with interior cells. This is not a problem since this information is not required (see Sect. 3.7).

3.7 Distributed-memory finite element assembly

Finally, we discuss the assembly of the global matrix $\mathbf{A}$ and vector $\mathbf{b}$ in a distributed-memory context. The assembly at the current sub-domain $D_{s}$ is performed with a loop in local cells in $\mathcal{T}^{\rm L}_{s}$. At each cell $T^{\ell}_{s}\in\mathcal{T}^{\rm L}_{s}$, we compute the corresponding elemental matrix and vector, namely $\mathbf{A}^{\ell}_{s}$ and $\mathbf{b}^{\ell}_{s}$, whose entries are defined as $[\mathbf{A}^{\ell}_{s}]_{ab}\doteq\mathrm{a}^{\ell}_{s}(\phi^{a}_{s,\ell},\phi^{b}_{s,\ell})$ and $[\mathbf{b}^{\ell}_{s}]_{a}\doteq\mathrm{b}^{\ell}_{s}(\phi^{a}_{s,\ell})$, where $\mathrm{a}^{\ell}_{s}(\cdot,\cdot)$ and $\mathrm{b}^{\ell}_{s}(\cdot)$ represent the restriction of the bi-linear form $\mathrm{a}(\cdot,\cdot)$ and the linear form $\mathrm{b}(\cdot)$ to cell $T^{\ell}_{s}$, resp. Then, the contributions of the elemental matrix and vector $\mathbf{A}^{\ell}_{s}$ and $\mathbf{b}^{\ell}_{s}$ are added to the corresponding entries of the global matrix and vector $\mathbf{A}$ and $\mathbf{b}$.

The assembly operations to be performed are closely related to the ones detailed in Sect. 2.5 for the serial case. In order to grasp the main differences between the serial and the distributed memory implementations, we rewrite the definition of $\mathcal{V}^{\rm agg}_{s}$ in Eq. (9) using the notation introduced in Sect. 3.5, namely

[TABLE]

where $\mathcal{M}_{s}(j)\subset\mathcal{I}^{\rm in}$ denotes the set containing the global DOF ids that constrain the value associated with a given sub-domain local id $j\in\mathcal{J}^{\rm out}_{s}$. The values $C^{jg}_{s}$, for $g\in\mathcal{M}_{s}(j)$ and $j\in\mathcal{J}^{\rm out}_{s}$, are the coefficients of the cell aggregation constraints that are needed to be computed in sub-domain $D_{s}$, and the value $[v^{g}]_{s}$ represents a variable that stores (in the local scope of sub-domain $D_{s}$) the value of $v^{g}$ for each master global DOF id $g\in\mathcal{M}_{s}(j)$ and each $j\in\mathcal{J}^{\rm out}_{s}$. These quantities are formally defined as

[TABLE]

with $\ell=L_{s}(k)$, $z=Z_{s}(k)$, $k=R_{s}(L^{\rm own}_{s}(j))$, and being $b$ such that $g=g^{b}_{s,\ell}$ or $g=\tilde{g}^{b}_{s,z}$. Note that the definition of these quantities takes one form or another depending on whether the global cell id $k=R_{s}(L^{\rm own}_{s}(j))$ belongs to $\mathcal{K}_{s}$ or not (i.e., if the root cell $T^{k}$ associated with $j$ is a locally relevant cell in $D_{s}$ or not). If $k\in\mathcal{K}_{s}$, the quantities can be computed from data stored in the local scope of the current processor. Otherwise, if $k\notin\mathcal{K}_{s}$ (which implies $k\in\mathcal{K}^{\rm R}_{s}$), the quantities need to be computed from the specially allocated variables $\tilde{g}^{b}_{s,z}$, $\tilde{x}^{b}_{s,z}$, and $\tilde{v}^{b}_{s,z}$ presented in Sect. 3.5.

The main differences between the serial and distributed-memory assembly are readily grasped by comparing formulas (2) and (10). The first main difference is the following. Note that, in both the serial and distributed implementations, the master DOFs are represented with global DOF ids. However, the constrained DOFs are represented with global ids in the serial implementation and with sub-domain local ids in the distributed-memory version. For this reason, it is not needed to generate a global numbering for the constrained DOFs. The second major difference is that, in the distributed memory implementation, one has to import data from remote sub-domains in order to compute the quantities $\mathcal{M}_{s}(j)$, $C^{jg}_{s}$, and $[v^{g}]_{s}$, for $g\in\mathcal{M}_{s}(j)$, if the root cell associated with the constrained DOF id $j\in\mathcal{J}^{\rm out}_{s}$ is not a locally relevant cell. The assembly operations to be performed in the distributed-memory case are obtained by modifying the operations of the serial case (Table 1) in accordance with these two main differences and taking into account that cells are represented by sub-domain local ids in the parallel implementation. The result is detailed in Table 2.

4 Numerical examples

4.1 Experimental environment

The numerical experiments are run at the Marenostrum-IV supercomputer [50], hosted by the Barcelona Supercomputing Center. The Marenostrum-IV is a petascale machine equipped with 3,456 compute nodes interconnected with the Intel OPA HPC network. Each node has 2x Intel Xeon Platinum 8160 multi-core CPUs, with 24 cores each (i.e. 48 cores per node) and 96 GBytes of RAM.

With respect to the software, we used the MPI-parallel implementation of the AgFEM method available at FEMPAR [35] linked against p4est v2.0 [41] as the Cartesian grid manipulation engine, and PETSc v3.9.0 [33] for distributed-memory linear algebra data structures and solvers. These software were compiled with Intel v18.0.5 compilers using system recommended optimization flags and linked against the Intel MPI Library (v2018.4.057) for message-passing and the BLAS/LAPACK available on the Intel MKL library for optimized dense linear algebra kernels. All floating-point operations were performed in IEEE double precision.

4.2 Model problem

As a model problem for the numerical examples, let us consider the Poisson equation with Dirichlet boundary conditions. After scaling with the diffusion term, the equation reads: find $u\in H^{1}(\Omega)$ such that

[TABLE]

where $f\in H^{-1}(\Omega)$ is the source term and $g\in H^{1/2}(\Gamma)$ is the prescribed value on the Dirichlet boundary. We denote by $\mathcal{V}^{\rm x}$ any of the FE spaces $\mathcal{V}^{\mathrm{std}}$, or $\mathcal{V}^{\mathrm{agg}}$, when it is not necessary to distinguish between them. We approximate problem  (12) in $\mathcal{V}^{\rm x}$ with the following variational equation: find $u_{h}\in\mathcal{V}^{\rm x}$ such that $\mathrm{a}(v_{h},u_{h})=\mathrm{b}(v_{h})$ for all $v_{h}\in\mathcal{V}^{\rm x}$, where

[TABLE]

with $n$ being the outward unit normal on $\Gamma$. Note that the forms $\mathrm{a}(\cdot,\cdot)$ and $\mathrm{b}(\cdot)$ include the usual terms resulting from the integration by parts of (12) plus additional terms associated with the weak imposition of Dirichlet boundary conditions with Nitsche’s method [51, 52]. The coefficient $\tau>0$ is a mesh-dependent parameter that has to be large enough to ensure the coercivity of $\mathrm{a}(\cdot,\cdot)$. We consider a weak imposition of boundary conditions since it is not straightforward to include prescribed values in the approximation space in a strong manner, when dealing with an unfitted grid. Nitsche’s method is commonly used in the embedded boundary community for circumventing this problem (see, e.g., [7, 9]) since it provides a consistent numerical scheme with optimal converge rates (also for high-order FEs).

When using the agFE space (i.e., $\mathcal{V}^{\rm x}=\mathcal{V}^{\rm agg}$), we compute the mesh-dependent parameter $\tau$ as

[TABLE]

where $h$ is the cell size in the computational grid and $\beta$ is a user defined constant parameter. In the numerical examples, we take two different values of $\beta$, namely $\beta=10.0$ and $\beta=100.0$. The choice of $\beta=10.0$ is common in the literature of unfitted methods (see, e.g., reference [53]). We also have considered $\beta=100.0$ in order to study the sensibility of the method with respect to this parameter. We note that using the space $\mathcal{V}^{\rm agg}$ and formula (14) for the computation of $\tau$, the bilinear form in (13) is coercive independently of the position of the cuts. We refer to [10] for the mathematical proof of this statement. However, when using the standard FE space $\mathcal{V}^{\rm std}$, formula (14) does not necessarily lead to a coercive bi-linear form. For the simulations with the standard space $\mathcal{V}^{\rm std}$, we compute $\tau$ following an alternative approach in order recover a coercive problem independently of the position of the cuts (see [15] for details). For the sake of completeness, we briefly present this methodology here. It consists in building the penalty parameter $\tau$ as an element-wise constant function, where the value in a generic cut cell, namely $\tau_{k}\doteq\tau|_{T^{k}}$, $k\in\mathcal{K}^{\rm cut}$, is computed as $\tau_{k}=\beta\lambda^{\rm max}_{k}$. In this context, $\beta$ is also a user-defined constant (with values $\beta=10.0$ and $\beta=100.0$ in the numerical examples). On the other hand, $\lambda^{\rm max}_{k}$ is the maximum eigenvalue of the following generalized eigenvalue problem: find $\eta_{k}\in\mathcal{V}(T^{k})$ and $\lambda_{k}\in\mathbb{R}$ such that

[TABLE]

Note that this problem is local at each cut cell, and its definition only depends on how the physical domain intersects the cut cell and the functional space used within the cell. We refer to [15] for the specific details on the solution of this problem. As indicated in [15], one of the main drawbacks of using the standard FE space $\mathcal{V}^{\rm std}$ is that the optimal value of $\tau$ needed for ensuring coercivity of the bilinear form $\mathrm{a}(\cdot,\cdot)$ (i.e., $\lambda^{\rm max}_{k}$) might be arbitrarily large if the active part of a cell $T\in\mathcal{T}$, namely $\Omega\cap T$, has a volume that tends to zero. This fact leads to strongly ill-conditioned systems of linear equations when considering the standard space $\mathcal{V}^{\rm std}$. This issue is solved when considering the agFE space $\mathcal{V}^{\rm agg}$ as the optimal value of $\tau$, which leads to a coercive bi-linear form, is independent of the position of the cuts (see [10] for further details).

4.3 Linear solver setup

The FE discretization of the weak problem represented by the bi-linear and linear forms in (13) leads to a system of linear algebraic equations, namely $\mathbf{A}\mathbf{x}=\mathbf{b}$, which, at large scales, requires efficient parallel linear solvers. Our goal is to show that the usage of the agFE spaces allows one to efficiently solve the linear systems using solvers that are well established for conventional FE analysis based on body-fitted grids. To this end, we consider the widely available linear solvers available in the PETSc library [33]. In particular, we use a conjugate gradient method from the KSP module of PETSc, preconditioned with a smoothed-aggregation AMG preconditioner called GAMG (see [54] for specific details).

We set up the linear solver by relying as much as possible on the default configuration given by GAMG to effectively show that the agFE spaces lead to linear systems that are efficiently solved using standard multigrid tools. The modifications that we have introduced in the default parameters are, basically, in order to build a preconditioner tailored for symmetric definite positive matrices (since the default configuration of GAMG is designed for the non-symmetric case). We report the configuration file used to set up the PETSc solver and preconditioner in Listing 1. The file is written using a domain specific language provided by the PETSc library (see the PETSc users’ manual [33] for further details on the file syntax).

In line 1 of Listing 1, we select a conjugate gradient method to compute an approximation $\mathbf{x}^{\rm cg}\approx\mathbf{x}$ of the solution of the symmetric definite positive system $\mathbf{A}\mathbf{x}=\mathbf{b}$. We declare convergence when $\|\mathbf{r}\|_{2}/\|\mathbf{b}\|_{2}<10^{-6}$ within the first $500$ iterations (see lines  2-5), where $\mathbf{r}\doteq\mathbf{b}-\mathbf{A}\mathbf{x}^{\rm cg}$ is the un-preconditioned residual and $\|\cdot\|_{2}$ denotes the Euclidean norm. In lines 6-7, we select the smoothed aggregation AMG preconditioner called GAMG and in line 8 we select a Cholesky factorization in order to solve the problem on the coarsest level of the AMG hierarchy. The GAMG preconditioner uses Chebyshev smoothers in order to approximate the solution at the finest and intermediate levels of the AMG hierarchy. In order to set up a Chebyshev smoother one needs to estimate the maximum eigenvalue of the corresponding system matrix (see, e.g., [55] for details). To compute this eigenvalue estimate, GAMG uses internally an iterative Krylov sub-space method. In line 9, we instruct GAMG to use a conjugate gradient iteration for this purpose. In particular, we are using the default number of iterations provided by PETSc (i.e., 10 iterations) in the conjugate gradient method used to estimate the eigenvalue. In the process of generating the prolongation and restriction operators associated with the AMG hierarchy, GAMG allows one to use the graph of the matrix $\mathbf{A}^{T}\mathbf{A}$ instead of the graph of the original matrix $\mathbf{A}$. By default, GAMG uses the graph of $\mathbf{A}^{T}\mathbf{A}$ in the finest level and the graph of (the coarse representation of) $\mathbf{A}$ at the other levels. In line 10, we instruct GAMG to use the graph of $\mathbf{A}$ also in the finest level. We have observed that this setting leads to a faster solver (since the product $\mathbf{A}^{T}\mathbf{A}$ is not computed) even though the complexity of the AMG preconditioner is slightly increased. Finally, we have included the option in line 11 in order to bypass a software BUG in the Intel MPI implementation used in the experiments. Please refer to the GAMG package documentation [39] for the default methods and parameters that are not covered in Listing 1.

4.4 Problem setup

In this section, we detail the rest of the parameter values and methods used in the numerical examples. In order to show that the methods are robust enough to deal with different geometrical data, we have performed all the experiments for three different complex geometries. On the one hand, we have considered a bulky 3D body whose shape reminds the one of a popcorn flake (see Fig. 7(a)). This "popcorn-flake" geometry is often used in the literature to study the performance of unfitted FE methods (see, e.g., [7]). The popcorn flake geometry considered here is obtained by taking the one defined in [7], scaling it by a factor of 0.5 and translating it by a value of 0.5 in each direction such that the body fits in the unit cube $[0,1]^{3}$. The second geometry is a massive spiral pipe (see Fig. 7(b)). The radius of its tubular cross section is 0.1, whereas the radius of the spiral central axis is 0.875. The third geometry, called "Swiss cheese", is a complex body often used is the literature of unfitted FE methods. The precise definition of this geometry is found in [7]. In all cases, the bounding box used to define the background mesh is the unit cube $[0,1]^{3}$, which is also depicted in Fig. 7.

The source term function $f$ and the Dirichlet condition function $g$ of the Poisson problem (12) are defined using the method of the manufactured analytic solutions so that problem (12) has the following exact solution: $u(x,y,z)=x+y+z$. This exact solution is used to check the quality of the computed FE approximations. Since the exact solution $u$ belongs to the considered FE spaces, the computed FE approximations $u_{h}$ have to coincide with $u$ up to the tolerance used to declare convergence in the linear solver. In order to check that this is true in our experiments, we have computed the (relative) error norms $|u-u_{h}|_{H^{1}}/|u|_{H^{1}}$ and $\|u-u_{h}\|_{L^{2}}/\ \|u\|_{L^{2}}$, where

[TABLE]

are the $H^{1}$ semi-norm and $L^{2}$ norm, resp.

As indicated in Sect. 2.1, we consider background meshes of hexahedral cells generated and partitioned in parallel with the p4est library. We generate the meshes by recursive sub-division of the unit cube $[0,1]^{3}$. That is, in the first refinement level, the cube is decomposed into $(2^{1})^{3}=8$ cells, in the second refinement level into $(2^{2})^{3}=64$ cells, etc. For the weak scaling test shown below, we generate three families of partitions such that the average number of cells per sub-domain is constant as the mesh is refined (see Table 3). The first one (referred as load 1) has approximately $1,000$ cells per sub-domain, and the second and the third ones (load 2 and load 3) have approximately $8,000$ and $64,000$ cells per sub-domain, resp. The Swiss cheese geometry is not simulated with the coarsest mesh since it is unable to capture the main geometrical features. The finest mesh considered in the examples has $1,073,741,824$ cells (corresponding to refinement level 10), and it is partitioned into $16,777$ sub-domains (which are mapped to the same number of MPI tasks). We used the default MPI task placement policy of Intel MPI (v2018.4.057) with partially filled nodes. The number of nodes $N$ is selected as $N=\texttt{ceil}(|\mathcal{S}|/48)$, where $\texttt{ceil}(x)$ rounds $x$ to the nearest integer that is greater or equal to $x$. If $|\mathcal{S}|$ is not multiple of $48$, the placement policy fully populates the first $N-1$ nodes with $48$ MPI tasks per node. The remaining $\texttt{mod}(|\mathcal{S}|,48)$ MPI tasks are placed on the last node, with mod denoting integer division remainder.

Each mesh is decomposed using the load-balancing mechanism of the p4est library. By default, the method splits the mesh such that the number of cells in each of the resulting sub-domains is approximately equally balanced. p4est also includes the option to provide a user-defined weight to each of the cells. In this case, the algorithm leads to partitions, where the sum of the weights in each sub-domain is equally balanced across the entire partition. We use this functionality in order to give different weights to active and exterior cells. In this context, one has to find a trade-off between two conflicting goals. On the one hand, it is desirable to load balance the storage of the background mesh and the cost of the geometrical computations that run on top of it. In this case, one has to assign the same weight to active and exterior cells. On the other hand, it is also desirable to load balance the assembly, storage, and solution of the underlying systems of linear algebraic equations. In that case, exterior cells do not play any role, and one has to balance the number of active cells per sub-domain. This is achieved in practice by assigning a much smaller weight to exterior cells than to active ones. In the experiments, we have used a ratio $w=10/1$ (i.e., weight $10$ for active cells and $1$ for exterior ones). We observed experimentally (see Figs. 8 and 9) that this value of $w$ leads to the best trade-off among the two aforementioned factors for the considered geometries. The wall clock time of the main phases of the AgFEM method (Fig. 8) and the wall clock time of the linear solver (Fig. 9) decreases significantly when the weight factor $w$ moves from $w=1$ to $w=10$. Beyond this point, the timings tend to stagnate or slightly increase for $w\geq 10$. An important exception is the assembly of the algebraic linear system, whose wall clock time significantly increases again for $w\geq 10$ (see the orange curve with pentagonal markers in Fig. 8). For this reason, we have chosen a value of $w=10$ for the subsequent numerical experiments.

We build the FE spaces $\mathcal{V}^{\rm std}$ and $\mathcal{V}^{\rm agg}$ using hexahedral elements with continuous piece-wise trilinear shape functions. We consider both agFE and conventional spaces in order to evaluate the benefits of using the AgFEM method. A summary of the main parameters and methods considered in the numerical examples is given in Table 4.

4.5 Weak scaling test

We start this section devoted to the weak scaling test by checking that the problems were accurately solved in the experiments. To this end, Fig. 10 shows the relative $L^{2}$ norm of the error associated with the computed numerical solutions versus problem size for the three loads per processor detailed in Table 3. The results in Figs. 10(a), 10(b), and 10(c) are for the "popcorn flake", the spiral and the "Swiss cheese" geometries. Note that the results are qualitatively equivalent for the three geometries. In Fig. 10, some points of the red lines (the ones obtained with standard FE spaces) are missing since the linear solver was not able to provide a solution. Either the setup of the preconditioner failed or the computed preconditioner was not positive definite (which makes it unusable in a conjugate gradients iteration). These problems are related with the strong ill-conditioning that results from the usage of standard FE spaces. Note, however, that the linear solvers properly worked in all cases, when considering the agFE spaces. Note that (as expected) the relative errors of the computed numerical solution are of the same order of magnitude as the tolerance imposed to declare convergence of the linear solver. Moreover, the computed errors in the $L^{2}$ norm are reduced about three orders of magnitude, if the problem is solved by setting a linear solver tolerance three orders of magnitude smaller, i.e., $10^{-9}$ instead of $10^{-6}$ (see Fig. 11). This clearly demonstrates that the solution quality is dictated by the linear solver tolerance as one would expect. Results for the $H^{1}$ norm of the error are qualitatively equivalent as the ones showed in Fig. 10 and they are not included for the sake of brevity.

We start the reporting of the weak scaling test, by discussing the impact of the different parameters on the algorithmic performance of the considered AMG preconditioners. First, we consider the impact of using agFE or standard FE spaces in the number of solver iterations. As the results in Fig. 12 show, it is clear that the usage of AgFEM is beneficial (and essential) for the good performance of the linear solver. As it was already discussed for Fig. 10, some results in Fig. 12 (the ones associated with the standard FE spaces) are missing since the linear solver was not able to provide a converged solution. This shows that standard FE spaces are unusable in practice since they do not lead to a robust discretization method. In contrast, the usage of agFE spaces allows one to effectively solve the underlying linear systems. In particular, note that the results obtained with agFE spaces (blue lines) are very close to the expected optimal performance of multigrid methods (i.e., number of linear solver iterations is asymptotically independent of the problem size). This optimal behavior is observed for all considered geometries, which further demonstrates the usability and robustness on the computational strategy based on agFE spaces.

Next, we study the influence the user-defined coefficient $\beta$ associated with Nitsche’s method (see Sect. 4.2) has in the performance of the solver. Since it is clear from Figs. 10 and 12 that the method based on standard FE spaces is unusable in practice, from now on, we only provide results associated with the agFE spaces. The results in Fig. 13 show the influence of $\beta$ in the solver performance for the two values of $\beta$ previously detailed in Table 4 and for the three loads per processor detailed in Table 3. As it is seen in the figure, the number of linear solver iterations are weakly scalable for the two values of $\beta$ (i.e., iterations asymptotically independent of the problem size). However, the absolute number of iterations is significantly larger for $\beta=100$ than for $\beta=10$, which suggests that $\beta$ has to be chosen large enough to ensure coercivity of the underlying bi-linear form, but cannot be chosen arbitrary large since it has a negative impact on the solver performance. This requirement is well known in the literature of Nitsche-based FE methods (see, e.g., [7]).

We conclude the numerical examples by reporting the wall clock time consumed in the main phases of the AgFEM method (see Fig. 14) and the one in the linear solver step (see Fig. 15). As it is observed in Fig. 14, all phases have remarkable weak scaling. For the linear solver (see Fig. 15), we report the wall clock times of the solver setup (i.e., the setup of the AMG preconditioner) and the wall clock time of the solver run (i.e., the preconditioned conjugate gradient iterations). Note that the scaling of the linear solver is not as optimal as for the different phases reported in Fig. 14. However, it can still be considered excellent since the solver wall clock times increase only moderately as the problem is scaled, and also if one takes into account that we have customized the linear solver as little as possible. As an illustrative example, for the popcorn flake body and for load 3, the total solver wall clock time (setup plus iterations) scales from $0.36$ to $4.33$ seconds, while the problem size scales from $67,192$ to $298,911,096$ DOFs. That is to say, the total solver time increases by a factor of $12x$, while the problem size increases by a factor of $4,449x$, which are excellent results.

4.6 Weak scaling test for the Poisson equation on the unit cube

In order to demonstrate that the (moderate) degradation of the linear solver time observed in Fig. 15 is not caused by the AgFEM method, we have performed (as a reference) a weak scaling test with a Poisson problem defined on the unit cube, with a body-fitted mesh optimally partitioned into sub-domains by means of a Cartesian partition. The test is done for three different local loads (referred to as "load 1", "load 2" and "load 3") with $20^{3}$, $30^{3}$, and $40^{3}$ cells per sub-domain respectively. The largest mesh considered has $1,053,696,000$ cells which are distributed over $16,464$ sub-domains (see Table 5). The linear solver considered here is the same as the one used before in Section 4.5 (i.e., we have used the configuration file displayed in Listing 1). The results of this weak scaling test are reported in Fig. 16. Note that the number of iterations of the conjugate gradient solver exhibits the optimal scaling behavior expected for this problem type (see Fig. 16(a)). However, even though this is the most favorable setup one can consider for an AMG preconditioner (i.e., Poisson problem, Cartesian body fitted meshes, optimal load balance, no AgFEM constraints), the obtained weak scaling of the solver time is not the expected optimal one (constant time as the problem size increases in the same proportion as the number of processors). This can be clearly seen in Figs. 16(b) and 16(c). These results show that the moderate loss of scalability previously observed also for the computations based on the AgFEM method in Fig. 15 cannot be directly attributed to the AgFEM method itself. The degradation of the linear solver time could be caused by other factors such as the underlying parallel system (i.e., the HPC cluster architecture plus the software stack that runs on top of it), but confirming this is out of the scope of the current work. In any case, by performing the weak scaling test on the unit cube with a body-fitted mesh, we have discarded that the loss of efficiency is caused by the AgFEM method.

5 Conclusions

In this work, we have presented the distributed-memory implementation and performance of the so-called AgFEM method for the solution of large-scale problems with unfitted FE techniques. The method is based on removal of basis functions associated with badly cut cells by introducing carefully designed constraints, which result in well-conditioned systems of linear algebraic equations independently of the position of the cuts. As the main contributions of the paper, (a) we have proposed a parallel implementation of AgFEM, (b) shown that the proposed parallel setup is scalable, and (c) shown that, by using AgFEM, the underlying system of linear equations can be efficiently solved using preconditioned iterative linear solvers designed for conventional FE analysis. In this case, we have considered well-known AMG preconditioners available in PETSc. This is in contrast to previous works, which considered highly customized preconditioners in order to handle the severely ill-conditioned operators resulting from the discretization based on unfitted grids. Also in contrast to previous works, which mainly consider non-scalable serial algorithms, we have studied the performance and scalability of the method with weak scaling tests up to 16K processors and up to 300M DOFs. The proposed parallel AgFEM method can be easily incorporated in existing large-scale codes since (a) it allows one to use standard parallel linear solvers, and (b) the algorithmic phases of the method can be implemented using standard functionality available in distributed-memory FE codes. In this case, the parallel implementation of the AgFEM method has been done in the large-scale FE package FEMPAR.

Acknowledgments

Financial support from the European Commission under the FET-HPC ExaQUte project (Grant agreement ID: 800898) within the Horizon 2020 Framework Programme is gratefully acknowledged. This work has been partially funded by the project MTM2014-60713-P from the “Ministerio de Economía, industria y Competitividad” of Spain. S. Badia gratefully acknowledges the support received from the Catalan Government through the ICREA Acadèmia Research Program. F. Verdugo gratefully acknowledges the support received from the Secretaria d’Universitats i Recerca of the Catalan Government in the framework of the Beatriu Pinós Program (Grant Id.: 2016 BP 00145). The authors thankfully acknowledge the computer resources at Marenostrum-IV and the technical support provided by the Barcelona Supercomputing Center (RES-ActivityID: FI-2018-1-0014, FI-2018-2-0009, FI-2018-3-0029, FI-2019-1-0007).