BDDC: Balancing Domain Decomposition by Constraints

• BDDC is a domain decomposition preconditioner that partitions a computational domain into nonoverlapping subdomains to efficiently solve large sparse systems from PDE discretizations.
• It employs a dual-primal substructuring method by combining local harmonic extensions with coarse corrections through selected primal constraints and weighted averaging.
• Its adaptive, multilevel, and surrogate variants enhance scalability and robustness across diverse discretization methods and challenging computational settings.

Balancing Domain Decomposition by Constraints (BDDC) is a nonoverlapping domain decomposition preconditioner for large sparse linear systems arising from PDE discretizations. In its standard form, the domain is partitioned into subdomains, interior unknowns are eliminated by iterative substructuring, and the remaining interface problem is preconditioned by combining local subdomain solves with a coarse correction defined by selected primal constraints. A characteristic feature of BDDC is the use of a partially assembled space, in which continuity is enforced only on a small set of coarse degrees of freedom, followed by a weighted averaging operator that reconstructs a globally conforming interface field. Across the literature, this structure appears in classical conforming finite elements, mixed and hybrid formulations, DG and HDG methods, virtual elements, and saddle-point systems, with polylogarithmic condition-number estimates in many cases and a broad range of parallel implementations (Šístek et al., 2011, Brenner et al., 2018, Bertoluzza et al., 2017).

1. Foundational construction and interface algebra

The algebraic core of BDDC is formulated after decomposition of the computational domain into nonoverlapping subdomains $\Omega_i$ and elimination of subdomain-interior degrees of freedom. In the iterative-substructuring formulation used for Stokes flow, the reordered discrete system yields the interface Schur complement problem

$$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$

The subdomain-local discrete harmonic condition is

$A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$

so the interface unknown determines an energy-minimizing extension into each subdomain interior (Šístek et al., 2011).

The relevant interface spaces are typically denoted $W_i$ for the subdomain harmonic traces, $W=W_1\times\cdots\times W_N$ for the fully broken interface space, and $U_\Gamma\subset W$ for the globally assembled interface space. BDDC introduces an intermediate partially assembled space

$U_\Gamma\subset \widetilde W\subset W,$

whose defining property is that functions are generally discontinuous on the interface except at the selected primal variables. In the global-matrix formulation used in the parallel Stokes implementation, the preconditioner acts as

$M_{BDDC}:r\mapsto v=Ew,$

where $w\in\widetilde W$ solves

$$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$

equivalently

$$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$0

Here $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$1 is the operator on the partially assembled space and $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$2 is the averaging or projection operator that maps back to the assembled interface space (Šístek et al., 2011).

The same architecture reappears in other formulations with different notation. For 3D divergence-free VEM Stokes, the interface space is written as $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$3, and the BDDC operator takes the canonical form

$$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$4

with $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$5 the partially assembled Schur complement and $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$6 a scaled restriction operator (Bevilacqua et al., 2023). For DG cell-by-cell cardiac models, the same formula is used on a broken interface $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$7,

$$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$8

but continuity is enforced only with respect to the BDDC splitting, not across the physical cell membranes, because the underlying global solution is intentionally discontinuous (Huynh et al., 2022).

These formulations make explicit that BDDC is not merely a coarse-grid correction. It is a dual-primal substructuring method in which the coarse problem, the local harmonic extensions, and the averaging operator are defined simultaneously through the partially assembled space.

2. Primal constraints, scaling operators, and spectral estimates

The defining design choice in BDDC is the selection of primal constraints. The minimal choice in several settings consists of corners or subdomain vertices. Enriched coarse spaces add edge averages, face averages, or related flux-preserving constraints. For the $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$9 interior penalty discretization of the biharmonic problem, the coarse variables are associated with subdomain corners, and the global harmonic space is split into functions vanishing at all corners and an energy-orthogonal coarse subspace (Brenner et al., 2018). In 3D VEM Stokes, all vertices are primal and additional primal constraints are imposed so that the averaging operator preserves normal fluxes across faces and edges; in some geometries, additional tangential edge averages may be needed (Bevilacqua et al., 2023).

The averaging operator is equally central. In the global-matrix Stokes implementation, the weights are “the inverse to the number of subdomains, in which the degree of freedom is contained,” and this scaling is used for both velocity and pressure unknowns (Šístek et al., 2011). More generally, the literature uses multiplicity scaling, $A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$0-scaling, deluxe scaling, and problem-specific weighted averages. In VEM for elliptic problems, coefficient-aware multiplicity scaling is written

$A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$1

which yields robustness under jumps in the diffusion coefficient (Bertoluzza et al., 2017). In cardiac EMI and related DG settings, $A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$2-scaling is defined pointwise by conductivity ratios, for example

$A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$3

with the corresponding jump and averaging operators $A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$4 and $A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$5 (Goebel et al., 11 Feb 2025).

In many SPD settings, the spectral estimate has the classical polylogarithmic form

$A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$6

For the $A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$7 interior penalty biharmonic problem, the paper proves

$A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$8

with $A_{11}u_{\Gamma 1}+A_{12}u_{\Gamma 2}=0,$9 independent of $W_i$0, $W_i$1, and the number of subdomains (Brenner et al., 2018). For conforming VEM discretizations of scalar diffusion, the estimate becomes

$W_i$2

independent of the number of subdomains, the fine mesh size, and jumps in the coefficient $W_i$3 (Bertoluzza et al., 2017). For 3D divergence-free VEM Stokes, the condition number bound is

$W_i$4

where $W_i$5 is the discrete inf-sup constant (Bevilacqua et al., 2023).

The effect of coarse-space enrichment is repeatedly observed in computation. In the Bidomain solver, using only vertices gives condition numbers growing from $W_i$6 to $W_i$7 as $W_i$8 increases from $W_i$9 to $W=W_1\times\cdots\times W_N$0, while $W=W_1\times\cdots\times W_N$1 and $W=W_1\times\cdots\times W_N$2 remain around $W=W_1\times\cdots\times W_N$3–$W=W_1\times\cdots\times W_N$4 and $W=W_1\times\cdots\times W_N$5–$W=W_1\times\cdots\times W_N$6, respectively (Huynh et al., 2021). In the Stokes experiments with Taylor–Hood elements, adding edge and face averages significantly stabilizes iteration counts, whereas with only corner constraints the counts grow much more strongly (Šístek et al., 2011). The common interpretation is straightforward: enriched primal spaces remove low-energy interface modes more effectively, but enlarge the coarse problem.

3. Adaptive, multilevel, and surrogate variants

A substantial part of the BDDC literature concerns the construction of coarse spaces beyond fixed vertex/edge/face choices. One major line is adaptive BDDC based on local generalized eigenproblems. In high-order mortar discretizations, the adaptive coarse space is defined through interface transformation operators obtained from local generalized eigenvalue problems, and the resulting preconditioned operator satisfies

$W=W_1\times\cdots\times W_N$7

where $W=W_1\times\cdots\times W_N$8 is a user-defined tolerance and $W=W_1\times\cdots\times W_N$9 depends only on the maximum number of interfaces per subdomain (Peng et al., 2017). For weighted plane wave least-squares Helmholtz discretizations, the adaptive two-level algorithm also yields

$U_\Gamma\subset W$0

but the interface eigenproblems are complex-valued and primal constraints are selected using the complex modulus of eigenvalues (Peng et al., 2018). In mixed Darcy flow through heterogeneous porous media, the adaptive algorithm enriches the coarse space by adding face constraints associated with local eigenvectors whose eigenvalues exceed a target threshold $U_\Gamma\subset W$1; on the SPE 10 benchmark this sharply reduces condition numbers and CG iterations (Sousedík, 2019).

A second line is multilevel BDDC. When the number of substructures is large, the exact coarse solve becomes the scalability bottleneck. The multilevel construction replaces the coarse solve by another BDDC step, recursively, so that subdomains at one level become “elements” on the next. For $U_\Gamma\subset W$2 levels, the cited estimate is

$U_\Gamma\subset W$3

with $U_\Gamma\subset W$4 (Šístek et al., 2012). The estimate records the standard trade-off: more levels improve parallel scalability of the coarse solve but worsen convergence constants.

A third line introduces nonstandard coarse objects. In BDDC-SO, the interface is refined into subobjects—subedges, subfaces, and vertices between pairs of subedges—and the condition-number bound becomes

$U_\Gamma\subset W$5

where $U_\Gamma\subset W$6 is the characteristic size of the subobjects rather than the original subdomain diameter $U_\Gamma\subset W$7. Since $U_\Gamma\subset W$8 can be chosen almost freely, the paper states that the condition number can theoretically be as small as $U_\Gamma\subset W$9, at the price of a larger coarse problem (Badia et al., 2020).

A fourth line replaces repeated adaptive setup by learned or stochastic surrogates. For stochastic elliptic problems with oscillatory and high-contrast coefficients, a deep neural network is trained to map truncated Karhunen–Loève variables to the dominant eigenvectors that define the adaptive coarse space; the predicted coarse space is then used in the standard adaptive BDDC algorithm (Chung et al., 2021). In stochastic BDDC for repeated sampling of elliptic problems with random coefficients, the algorithm is split into offline and online stages: local Karhunen–Loève parametrizations and Polynomial Chaos expansions are built offline, while the sample-dependent BDDC preconditioner is assembled online at small cost (Tu et al., 7 Oct 2025).

Variant	Distinguishing mechanism	Representative statement
Adaptive BDDC	local generalized eigenproblems	$U_\Gamma\subset \widetilde W\subset W,$0 (Peng et al., 2017)
Multilevel BDDC	recursive coarse-problem replacement	coarse bottleneck reduced, more iterations (Šístek et al., 2012)
BDDC-SO	constraints on subobjects	$U_\Gamma\subset \widetilde W\subset W,$1 (Badia et al., 2020)
Learned/stochastic BDDC	DNN or PC/KL surrogate coarse data	cheaper online stage (Chung et al., 2021, Tu et al., 7 Oct 2025)

4. Extensions beyond classical conforming SPD elliptic systems

Although BDDC is often associated with SPD elliptic problems, the published formulations extend far beyond that setting. The key structural point is that BDDC is usually applied after static condensation or iterative substructuring to an interface operator whose properties may differ significantly from those of the original system.

For incompressible Stokes flow discretized by Taylor–Hood elements, the original bilinear form is symmetric but indefinite, and the assembled linear system is a saddle-point system. The cited implementation factorizes the global matrices with MUMPS using $U_\Gamma\subset \widetilde W\subset W,$2, applies the BDDC preconditioner in global algebraic form, and reports outer iterations by a parallel PCG solver, while also showing MATLAB experiments with GMRES and BiCGStab (Šístek et al., 2011). The paper explicitly notes that the available theory of Li and Widlund for incompressible Stokes with benign subspaces requires discontinuous pressure approximation and therefore does not cover the Taylor–Hood and serendipity elements used in that work.

For the cardiac Bidomain equations, the Jacobian system at each Newton step is reduced by static condensation to an interface Schur complement $U_\Gamma\subset \widetilde W\subset W,$3 that is symmetric positive semidefinite, so PCG can be used there even though the full problem is nonlinear and coupled to ODEs. The BDDC preconditioner is combined with a Newton–Krylov outer iteration and deluxe scaling of the dual variables (Huynh et al., 2021). In the mixed-hybrid formulation of Darcy flow with combined mesh dimensions, the global and local systems are symmetric saddle-point systems with a nonzero penalty block $U_\Gamma\subset \widetilde W\subset W,$4, but iterative substructuring reduces them to an SPD interface problem on Lagrange multipliers, and PCG with BDDC becomes applicable (Šístek et al., 2015).

In DG and HDG settings, the interface unknowns are naturally duplicated. For cardiac cell-by-cell models and the 3D EMI model, this duplication is not an artifact of the preconditioner but a consequence of the physical model: the global solution is discontinuous across cell boundaries, and the BDDC spaces are designed to transfer information without imposing physical continuity across membranes (Huynh et al., 2022, Goebel et al., 11 Feb 2025). For HDG discretizations of elliptic optimal control, local elimination yields a global interface system for the numerical traces. In one paper this condensed system is nonsymmetric positive definite, and the BDDC-preconditioned solver is GMRES rather than CG (Liu et al., 19 Aug 2025); in another HDG optimal-control paper, edge-average primal constraints and multiplicity scaling are used in a BDDC preconditioner for the condensed interface operator, and numerical evidence is reported for robustness with respect to the regularization parameter $U_\Gamma\subset \widetilde W\subset W,$5 (Liu et al., 2 Apr 2025).

Block-coupled systems also admit BDDC components. For the three-field mixed finite element discretization of Biot’s consolidation model, elimination of displacement and interior pressure variables produces an SPD interface system for interface total pressure, interface pressure, and a displacement continuity multiplier. The preconditioner is block-structured: BDDC is used for the interface total-pressure block and the interface pressure block, while FETI-DP is used for the Lagrange multiplier block (Chu et al., 7 Apr 2025).

5. Discretization ecosystems and computational performance

The range of discretizations treated with BDDC is unusually wide. The cited literature includes conforming $U_\Gamma\subset \widetilde W\subset W,$6, $U_\Gamma\subset \widetilde W\subset W,$7, and Taylor–Hood finite elements; $U_\Gamma\subset \widetilde W\subset W,$8 interior penalty biharmonic discretizations; mortar finite elements; RT0 mixed and mixed-hybrid flow discretizations; weighted plane wave least-squares methods for Helmholtz; HDG and HHO on polytopal meshes; conforming and divergence-free VEM; and composite DG methods for microscopic cardiac models (Brenner et al., 2018, Peng et al., 2017, Sousedík, 2019, Peng et al., 2018, Badia et al., 13 Jun 2025, Bevilacqua et al., 2023, Huynh et al., 2022).

Several papers provide explicit parallel realizations. The global-matrix Stokes implementation is based on the MUMPS parallel sparse direct solver, with one subdomain assigned to one processor and distributed matrix storage (Šístek et al., 2011). The multilevel implementation is given in the MPI-based Fortran 95 library BDDCML, with serial MUMPS for subdomain problems and parallel MUMPS for the final coarse problem (Šístek et al., 2012). For mixed-hybrid porous-media flow, the implementation combines Flow123d and BDDCML (Šístek et al., 2015). For 3D VEM Stokes, the parallel solver uses PETSc/PCBDDC (Bevilacqua et al., 2023). For the 3D EMI model, BDDC is implemented in the Ginkgo library and used as a preconditioner for CG (Goebel et al., 11 Feb 2025). For non-conforming polytopal hybrid methods, the implementation is in GridapSolvers.jl with Gridap.jl and is run on the GADI supercomputer (Badia et al., 13 Jun 2025).

The numerical record is similarly diverse. In the $U_\Gamma\subset \widetilde W\subset W,$9 interior penalty biharmonic problem on 16 subdomains, the unpreconditioned system requires $M_{BDDC}:r\mapsto v=Ew,$0, $M_{BDDC}:r\mapsto v=Ew,$1, $M_{BDDC}:r\mapsto v=Ew,$2, $M_{BDDC}:r\mapsto v=Ew,$3, and $M_{BDDC}:r\mapsto v=Ew,$4 iterations as $M_{BDDC}:r\mapsto v=Ew,$5 is refined, while the BDDC-preconditioned system remains at about $M_{BDDC}:r\mapsto v=Ew,$6–$M_{BDDC}:r\mapsto v=Ew,$7 iterations; the preconditioned full-system condition number is around $M_{BDDC}:r\mapsto v=Ew,$8–$M_{BDDC}:r\mapsto v=Ew,$9, whereas the unpreconditioned $w\in\widetilde W$0 grows from about $w\in\widetilde W$1 to over $w\in\widetilde W$2 (Brenner et al., 2018). In the Stokes cavity tests, the parallel BDDC approach is reported to be about twice as fast as direct MUMPS on 8 processors in a larger cavity test, and in the arterial-flow example a small problem is solved in parallel in 30 seconds versus 133 seconds for the serial frontal solver (Šístek et al., 2011).

Large-scale scalability studies are equally prominent. For the Bidomain equations, weak scaling on the Galileo cluster is reported up to 2048 processors, with BDDC and FETI-DP described as scalable; on the slab case, BDDC linear iterations rise from $w\in\widetilde W$3 at 32 processors to $w\in\widetilde W$4 at 2048, while on the ellipsoid case they stay around $w\in\widetilde W$5–$w\in\widetilde W$6 over the same range. In strong scaling on the Indaco cluster, BDDC time per time step drops from $w\in\widetilde W$7 s at 32 processors to $w\in\widetilde W$8 s at 256 processors for the slab, and from $w\in\widetilde W$9 s to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$0 s for the ellipsoid (Huynh et al., 2021). For BDDC-SO on MareNostrum IV, weak-scaling tests reach more than 260 million unknowns on up to 8K processors; in the reported large cases, BDDC-SO converges in fewer than about 21 iterations and under 9 seconds, while standard BDDC may require over 1500 iterations and over 300 seconds (Badia et al., 2020).

These results do not identify a single universal performance regime. Instead, they show that BDDC can remain competitive across very different algebraic structures provided that the primal space, scaling, and coarse solver are chosen in a formulation-aware way.

6. Limitations, misconceptions, and current directions

A common misconception is that BDDC is only a solver for conforming SPD second-order elliptic problems. The published evidence does not support that restriction: BDDC has been developed for symmetric indefinite Stokes systems, nonlinear Newton Jacobians in Bidomain models, DG and HDG discretizations with duplicated traces, mixed Darcy systems with pressure and flux variables, VEM saddle-point systems, and nonsymmetric positive definite interface operators solved by GMRES (Šístek et al., 2011, Huynh et al., 2021, Liu et al., 19 Aug 2025). The accurate statement is narrower: the cleanest theory is often available for reduced interface operators with SPD or benign-subspace structure, while more general settings require problem-specific analysis.

A second misconception is that more primal constraints are unconditionally better. The literature repeatedly reports the opposite trade-off. Adding edge and face averages usually improves conditioning and stabilizes iteration counts, but it also enlarges the coarse problem (Huynh et al., 2021). Subobject constraints can reduce the theoretical condition number to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$1 with $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$2 close to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$3, even theoretically to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$4, but the paper emphasizes the corresponding increase in coarse-space size (Badia et al., 2020). In highly heterogeneous porous media, lowering the adaptive threshold $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$5 sharply improves convergence but increases the number of adaptive constraints, especially in 3D (Sousedík, 2019).

A third issue concerns scaling. Deluxe scaling is central in several analyses and can be decisive for robustness in heterogeneous settings, but it is not uniformly dominant in wall-clock time. In the Bidomain study, deluxe scaling is essential for the convergence proof yet “does not dramatically change the iteration counts” relative to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$6-scaling and can increase CPU time slightly because the weighted averaging is more expensive (Huynh et al., 2021). By contrast, in adaptive mortar and adaptive PWLS Helmholtz formulations, deluxe-type scaling significantly reduces the number of primal unknowns or improves robustness (Peng et al., 2017, Peng et al., 2018).

A fourth issue is coarse-problem scalability. The multilevel literature states the point explicitly: adding levels replaces a coarse bottleneck by more iterations. In the geocomposite elasticity example, moving from two to three to four levels changes the condition number from about $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$7 to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$8 to $$\widetilde a(w,z)=(r,Ez)\qquad \forall z\in\widetilde W,$$9 and the iterations from $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$00 to $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$01 to $$Su_{\Gamma 2}=g_2,\qquad S=A_{22}-A_{21}A_{11}^{-1}A_{12},\qquad g_2=f_2-A_{21}A_{11}^{-1}f_1.$$02, yet the best total solution time at the largest processor count is achieved with four levels (Šístek et al., 2012). This is not a contradiction; it is the expected trade-off between spectral quality and parallel coarse-solver cost.

Recent surrogate approaches introduce another limitation: transferability. In the learned adaptive BDDC study, the neural network generalizes well when test coefficients have similar stochastic properties to the training data, but prediction degrades substantially when the mean permeability changes significantly, as in the transition from modified SPE10 Layer 35 to Layer 34 (Chung et al., 2021). The stochastic BDDC framework based on local KL and Polynomial Chaos addresses repeated-solve costs more directly, but its analysis requires sufficiently accurate KL truncation and PC approximation (Tu et al., 7 Oct 2025).

This suggests that contemporary BDDC research is driven primarily by coarse-space design, scaling design, and coarse-solver design rather than by changes to the basic dual-primal architecture. The recurrent themes are robust adaptive constraints, multilevel and subobject coarse spaces, formulation-specific averaging for DG and hybrid methods, and stochastic or learned surrogates that reduce setup cost without destroying the spectral properties on which BDDC depends.