---
title: 'BDDC: Balancing Domain Decomposition by Constraints'
url: https://www.emergentmind.com/topics/balancing-domain-decomposition-by-constraints-bddc
type: topic
---

# BDDC: Balancing Domain Decomposition by Constraints

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 [1101.1775], [1811.06941], [1708.03599].

## 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 [1101.1775].

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
\[
v=E\widetilde A^{-1}E^Tr.
\]
Here \(\widetilde A\) is the operator on the partially assembled space and \(E\) is the averaging or projection operator that maps back to the assembled interface space [1101.1775].

The same architecture reappears in other formulations with different notation. For 3D divergence-free VEM Stokes, the interface space is written as \(\widetilde{\mathbf V}_\Gamma=\widehat{\mathbf V}_\Pi\oplus \mathbf V_\Delta\), and the BDDC operator takes the canonical form
\[
M^{-1}=\widetilde R_D^T\,\widetilde S^{-1}\,\widetilde R_D,
\]
with \(\widetilde S\) the partially assembled Schur complement and \(\widetilde R_D\) a scaled restriction operator [2304.09770]. For DG cell-by-cell cardiac models, the same formula is used on a broken interface \(\Gamma'\),
\[
M^{-1}_{\mathrm{BDDC}}=\widetilde R_{D,\Gamma'}^T(\widetilde S_{\Gamma'})^{-1}\widetilde R_{D,\Gamma'},
\]
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 [2212.12295].

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 \(C^0\) 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 [1811.06941]. 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 [2304.09770].

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 [1101.1775]. More generally, the literature uses multiplicity scaling, \(\rho\)-scaling, deluxe scaling, and problem-specific weighted averages. In VEM for elliptic problems, coefficient-aware multiplicity scaling is written
\[
d^{\ell,i}=\frac{\rho_\ell^\gamma}{\sum_{j\in\mathcal N_i}\rho_j^\gamma},\qquad \gamma\in[1/2,\infty),
\]
which yields robustness under jumps in the diffusion coefficient [1708.03599]. In cardiac EMI and related DG settings, \(\rho\)-scaling is defined pointwise by conductivity ratios, for example
\[
\delta_i^\dagger(x)=\frac{\sigma_i}{\sum_{j\in N_x}\sigma_j},
\]
with the corresponding jump and averaging operators \(P_D=I-E_D\) and \(E_D\) [2502.07722].

In many SPD settings, the spectral estimate has the classical polylogarithmic form
\[
\kappa \le C(1+\ln(H/h))^2.
\]
For the \(C^0\) interior penalty biharmonic problem, the paper proves
\[
\kappa(S_h)\le C(1+\ln(H/h))^2
\]
with \(C\) independent of \(h\), \(H\), and the number of subdomains [1811.06941]. For conforming VEM discretizations of scalar diffusion, the estimate becomes
\[
\omega_{\mathrm{BDDC}}\lesssim (1+\log(Hk^2/h))^2,
\]
independent of the number of subdomains, the fine mesh size, and jumps in the coefficient \(\rho\) [1708.03599]. For 3D divergence-free VEM Stokes, the condition number bound is
\[
\kappa(M^{-1}\widehat S)\le C\,\frac{1}{\beta_h^2}\left(1+\log\frac{H}{h}\right)^2,
\]
where \(\beta_h\) is the discrete inf-sup constant [2304.09770].

The effect of coarse-space enrichment is repeatedly observed in computation. In the Bidomain solver, using only vertices gives condition numbers growing from \(8.4\) to \(110.1\) as \(H/h\) increases from \(4\) to \(24\), while \(V+E\) and \(V+E+F\) remain around \(1.9\)–\(4.5\) and \(1.7\)–\(4.5\), respectively [2101.02959]. 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 [1101.1775]. 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
\[
\kappa(\hat G)\le C\Theta,
\]
where \(\Theta\ge 1\) is a user-defined tolerance and \(C\) depends only on the maximum number of interfaces per subdomain [1704.07674]. For weighted plane wave least-squares Helmholtz discretizations, the adaptive two-level algorithm also yields
\[
\kappa(\hat G)\le C\Theta,
\]
but the interface eigenproblems are complex-valued and primal constraints are selected using the complex modulus of eigenvalues [1801.08800]. 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 \(\tau\); on the SPE 10 benchmark this sharply reduces condition numbers and CG iterations [1901.02090].

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 \(L\) levels, the cited estimate is
\[
\kappa(M_{BDDC}\widehat S)\le \prod_{\ell=1}^{L-1} C_\ell\left(1+\log\frac{H_\ell}{H_{\ell-1}}\right)^2,
\]
with \(H_0\equiv h\) [1201.5991]. 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
\[
\kappa \le C(1+\log(L/h))^2,
\]
where \(L\) is the characteristic size of the subobjects rather than the original subdomain diameter \(H\). Since \(L\) can be chosen almost freely, the paper states that the condition number can theoretically be as small as \(O(1)\), at the price of a larger coarse problem [2001.07289].

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 [2104.09162]. 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 [2510.05993].

| Variant | Distinguishing mechanism | Representative statement |
|---|---|---|
| Adaptive BDDC | local generalized eigenproblems | \(\kappa(\hat G)\le C\Theta\) [1704.07674] |
| Multilevel BDDC | recursive coarse-problem replacement | coarse bottleneck reduced, more iterations [1201.5991] |
| BDDC-SO | constraints on subobjects | \(\kappa \le C(1+\log(L/h))^2\) [2001.07289] |
| Learned/stochastic BDDC | DNN or PC/KL surrogate coarse data | cheaper online stage [2104.09162], [2510.05993] |

## 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 \(LDL^T\), 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 [1101.1775]. 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 \(\widehat S_\Gamma\) 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 [2101.02959]. 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 \(C\), but iterative substructuring reduces them to an SPD interface problem on Lagrange multipliers, and PCG with BDDC becomes applicable [1504.07085].

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 [2212.12295], [2502.07722]. 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 [2508.13997]; 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 \(\beta\) [2504.02072].

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 [2504.04859].

## 5. Discretization ecosystems and computational performance

The range of discretizations treated with BDDC is unusually wide. The cited literature includes conforming \(Q_1\), \(Q_2\), and Taylor–Hood finite elements; \(C^0\) 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 [1811.06941], [1704.07674], [1901.02090], [1801.08800], [2506.11956], [2304.09770], [2212.12295].

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 [1101.1775]. 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 [1201.5991]. For mixed-hybrid porous-media flow, the implementation combines Flow123d and BDDCML [1504.07085]. For 3D VEM Stokes, the parallel solver uses PETSc/PCBDDC [2304.09770]. For the 3D EMI model, BDDC is implemented in the Ginkgo library and used as a preconditioner for CG [2502.07722]. For non-conforming polytopal hybrid methods, the implementation is in GridapSolvers.jl with Gridap.jl and is run on the GADI supercomputer [2506.11956].

The numerical record is similarly diverse. In the \(C^0\) interior penalty biharmonic problem on 16 subdomains, the unpreconditioned system requires \(95\), \(235\), \(434\), \(704\), and \(1026\) iterations as \(h\) is refined, while the BDDC-preconditioned system remains at about \(23\)–\(27\) iterations; the preconditioned full-system condition number is around \(13\)–\(19\), whereas the unpreconditioned \(A_h\) grows from about \(10^3\) to over \(4\times 10^5\) [1811.06941]. 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 [1101.1775].

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 \(22\) at 32 processors to \(78\) at 2048, while on the ellipsoid case they stay around \(30\)–\(37\) over the same range. In strong scaling on the Indaco cluster, BDDC time per time step drops from \(348.2\) s at 32 processors to \(8.9\) s at 256 processors for the slab, and from \(571.8\) s to \(15.1\) s for the ellipsoid [2101.02959]. 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 [2001.07289].

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 [1101.1775], [2101.02959], [2508.13997]. 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 [2101.02959]. Subobject constraints can reduce the theoretical condition number to \(C(1+\log(L/h))^2\) with \(L\) close to \(h\), even theoretically to \(O(1)\), but the paper emphasizes the corresponding increase in coarse-space size [2001.07289]. In highly heterogeneous porous media, lowering the adaptive threshold \(\tau\) sharply improves convergence but increases the number of adaptive constraints, especially in 3D [1901.02090].

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 \(\rho\)-scaling and can increase CPU time slightly because the weighted averaging is more expensive [2101.02959]. By contrast, in adaptive mortar and adaptive PWLS Helmholtz formulations, deluxe-type scaling significantly reduces the number of primal unknowns or improves robustness [1704.07674], [1801.08800].

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 \(50\) to \(79\) to \(568\) and the iterations from \(46\) to \(56\) to \(131\), yet the best total solution time at the largest processor count is achieved with four levels [1201.5991]. 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 [2104.09162]. 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 [2510.05993].

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.

Source: https://www.emergentmind.com/topics/balancing-domain-decomposition-by-constraints-bddc