Non-Iterative Domain Decomposition Integrator

• Non-iterative domain decomposition integrators are time-marching schemes that decompose spatial operators and update via fixed local solves without iterative coupling.
• They employ operator splitting, overlapping subdomains, and partition-of-unity weights to incorporate interface corrections and guarantee unconditional stability.
• These integrators have been applied to parabolic, fluid–structure, and acoustic wave problems, achieving robust convergence and, in some cases, second-order time accuracy.

In the cited literature, a non-iterative domain decomposition integrator denotes a domain decomposition scheme for an evolutionary problem in which the solution is advanced from one time level to the next by a fixed number of algebraic or subdomain-local substeps, rather than by inner Schwarz, fixed-point, or Krylov iterations on every time slab. The common construction begins from an evolution equation such as $\frac{du}{dt}+Au=f(t)$ or its finite element form $B\frac{dy}{dt}+Ay=\varphi(t)$, decomposes the spatial operator into subdomain contributions, and embeds the coupling into an operator-splitting, factorized, projection, or Robin-transmission update. This viewpoint is explicit in the regionally additive schemes for parabolic and first-order evolutionary equations, in additive Stokes solvers in primitive variables, in partitioned fluid–structure interaction algorithms, and in recent acoustic-wave integrators based on local prediction plus subdomain Crank–Nicolson solves (Vabishchevich et al., 2014, Vabishchevich, 2011, Vabishchevich, 2011, Seboldt et al., 2020, Buchholz et al., 25 Jul 2025).

1. Definition and distinguishing criteria

A central distinction in this literature is between domain decomposition as an iterative solver and domain decomposition as a time integrator. For time-dependent problems, the iteration-free approach constructs a splitting so that “a transition to a new time level is performed via solving problems in particular subdomains,” with no convergence loop at fixed time level (Vabishchevich, 2011). In the parabolic setting, this is described as an “iteration-free domain decomposition algorithm” in which the advance from $t^n$ to $t^{n+1}$ is realized by a fixed number of substeps after an additive decomposition of the operator and the right-hand side (Vabishchevich et al., 2014). The same formulation appears in the two-component overlapping scheme for parabolic equations, where each time step requires only a fixed sequence of local solves and “there is no inner Schwarz iteration or fixed-point iteration to convergence” (Vabishchevich, 2017).

The classification of such methods is tied to three choices: the method of domain decomposition, the choice of decomposition operators, and the splitting scheme employed (Vabishchevich et al., 2014). Within this framework, overlapping subdomain methods are preferred when one seeks homogeneous numerical algorithms, because the overlap and the associated weighting functions encode interface exchange directly inside the decomposed operators. Substructuring variants reorganize the same idea by separating interior subdomain nodes from common boundary or interface nodes, still within a fixed-step, non-iterative update (Vabishchevich, 2011).

This definition excludes several nearby but distinct classes. BDD and FETI, for example, are iterative Krylov substructuring solvers for interface problems, even when their iterates admit useful local error decompositions (Rey et al., 2013). Likewise, some nonlocal substructuring frameworks establish exact multi-domain equivalence without yet providing an iteration-free solver (Capodaglio et al., 2020). The term therefore has a stricter meaning than “domain decomposition without a global monolithic matrix”; it refers to the absence of inner domain-decomposition iterations in the update itself.

2. Abstract operator formulation

The basic analytical template is an evolution equation with a nonnegative operator,

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$

or, after finite element semidiscretization of a parabolic problem,

$B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$

with $B=B^*>0$ and $A=A^*>0$ (Vabishchevich et al., 2014). A standard two-level weighted scheme is then

$$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$

and, in the non-self-adjoint first-order setting, unconditional stability holds for $\sigma\ge \frac12$ with transition operator

$B\frac{dy}{dt}+Ay=\varphi(t)$0

(Vabishchevich, 2011). The self-adjoint parabolic theory uses the analogous weighted estimate and then embeds domain decomposition into the operator itself (Vabishchevich et al., 2014).

The decomposed form is additive: $B\frac{dy}{dt}+Ay=\varphi(t)$1 with $B\frac{dy}{dt}+Ay=\varphi(t)$2 in the standard parabolic formulation (Vabishchevich et al., 2014). For non-self-adjoint problems, the decomposition is applied after splitting

$B\frac{dy}{dt}+Ay=\varphi(t)$3

so that $B\frac{dy}{dt}+Ay=\varphi(t)$4 and $B\frac{dy}{dt}+Ay=\varphi(t)$5, and then constructing $B\frac{dy}{dt}+Ay=\varphi(t)$6 with $B\frac{dy}{dt}+Ay=\varphi(t)$7 preserved by design (Vabishchevich, 2011). This preservation of nonnegativity is the key algebraic property behind the unconditional stability of the regularized additive schemes.

A further technical device is symmetrization. In the parabolic finite element setting one introduces

$B\frac{dy}{dt}+Ay=\varphi(t)$8

which yields

$B\frac{dy}{dt}+Ay=\varphi(t)$9

(Vabishchevich et al., 2014). This form places the DD integrator inside standard operator-splitting stability theory.

3. Overlap, partitions of unity, and interface corrections

The most characteristic spatial construction is an overlapping covering

$t^n$0

together with partition-of-unity functions $t^n$1 such that

$t^n$2

(Vabishchevich et al., 2014). These weights localize both the operator and the forcing: $t^n$3 with the standard choice

$t^n$4

(Vabishchevich et al., 2014). Equivalent forms,

$t^n$5

show that the overlap is not merely geometric padding; it is part of the interface-transfer mechanism.

A different overlapping construction uses indicator functions rather than smooth weights. For two subdomains with overlap $t^n$6, define

$t^n$7

so that

$t^n$8

This yields the corrected decomposition

$t^n$9

where the overlap contribution is subtracted once so that it is not double counted (Vabishchevich, 2017). Rewriting with

$t^{n+1}$0

reduces the three-component overlap-corrected formulation to a standard positive two-component factorization (Vabishchevich, 2017).

Substructuring schemes isolate interface degrees of freedom even more explicitly. In the two-component interior/interface split, one chooses

$t^{n+1}$1

where $t^{n+1}$2 is the set of internal boundary nodes. The resulting operator split,

$t^{n+1}$3

separates the interior subdomain problem from the boundary/interface problem (Vabishchevich, 2011). This suggests a common structural interpretation across overlapping and substructured methods: interface exchange is encoded algebraically in the decomposition operators rather than enforced by repeated transmission iterations.

4. Factorized updates and stability theory

The classical non-iterative update for two subdomains is the factorized two-stage scheme. In its weighted variational form it reads

$t^{n+1}$4

$t^{n+1}$5

For $t^{n+1}$6 this gives the Peaceman–Rachford scheme, and for $t^{n+1}$7 the Douglas–Rachford scheme (Vabishchevich et al., 2014). In matrix form,

$t^{n+1}$8

which is the paper’s “regionally additive” factorization (Vabishchevich et al., 2014). The corresponding two-stage operator form in the two-component overlapping formulation is also described as the iteration-free analogue of ADI-type methods (Vabishchevich, 2017).

The principal stability statement is that the factorized regionally additive difference scheme is unconditionally stable for $t^{n+1}$9 (Vabishchevich et al., 2014). In symmetrized coordinates this becomes

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$0

with estimate

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$1

For the indicator-function three-component decomposition, unconditional stability again holds for $\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$2, with

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$3

where

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$4

(Vabishchevich, 2017).

For general first-order evolutionary equations with non-self-adjoint $\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$5, the same threshold $\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$6 underlies two regularized non-iterative families. The additive regularized transition operator is

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$7

and the multiplicative analogue is

$\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$8

both unconditionally stable for $\frac{du}{dt}+Au=f(t), \qquad u(0)=u^0,$9 (Vabishchevich, 2011). The same paper also develops vector additive schemes in which $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$0 components are advanced by successive inversions of operators $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$1, again without subdomain iterations.

For the unsteady Stokes equations in primitive variables, non-iterativity appears through a different split: a viscous evolution step followed by a pressure projection step. After decomposition into $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$2, the viscous substep is implemented by ordered triangular sweeps,

$B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$3

followed by additive componentwise pressure correction (Vabishchevich, 2011). Theorem 5.1 states unconditional stability of the additive domain decomposition scheme, with final estimate

$B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$4

(Vabishchevich, 2011).

Numerical behavior in the parabolic papers is more nuanced than the unconditional stability theory. Reported experiments show $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$5 behavior for larger $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$6 and $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$7 for smaller $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$8, decomposition error close to $B\frac{dy}{dt}+Ay=\varphi(t), \qquad y(0)=y^0,$9 in the reported mesh-refinement studies, and strong dependence on overlap width (Vabishchevich et al., 2014). In the indicator-function study, reducing the overlap width decreases accuracy, refining the spatial grid also decreases accuracy, and the results indicate conditional convergence of the domain decomposition approximations (Vabishchevich, 2017).

5. Representative equation classes and modern extensions

The non-iterative template has been adapted beyond self-adjoint parabolic diffusion to incompressible flow, fluid–structure interaction, and wave propagation. The common pattern remains a fixed-step coupling: a local or subdomain solve uses previously available interface data, then a complementary local solve uses the newly computed state, and no interface convergence loop is introduced.

Problem class	Characteristic integrator	Reported property
Parabolic equations	Factorized regionally additive Douglas–Rachford / Peaceman–Rachford updates	Unconditionally stable for $B=B^*>0$0
First-order evolutionary equations with non-self-adjoint operators	Additive, multiplicative, and vector regularized schemes	Unconditionally stable for $B=B^*>0$1
Stokes equations	Viscous step plus pressure projection in decomposed product space	Unconditionally stable
Thick-structure FSI	Sequential structure solve, geometry update, fluid solve with generalized Robin conditions	Unconditionally stable; $B=B^*>0$2 in time and optimal convergence in space
Acoustic wave equation	Interface prediction plus local Crank–Nicolson solves and averaging	Second-order in time; global convergence $B=B^*>0$3 under a CFL-type condition

In moving-domain fluid–structure interaction, the non-iterative construction is based on generalized Robin transmission conditions obtained by combining the kinematic and dynamic interface conditions with a parameter $B=B^*>0$4: $B=B^*>0$5 At each step, one performs: a structure solve using the previous fluid data, a geometry update by harmonic extension into the fluid domain, and a fluid solve on the updated domain using the new structure state (Seboldt et al., 2020). The method is non-iterative because each time step requires only one structure solve, one geometry update, and one fluid solve. Its semi-discrete stability theorem provides an energy estimate with total discrete energy

$B=B^*>0$6

and the convergence analysis proves $B=B^*>0$7 convergence in time and optimal spatial convergence under the stated regularity assumptions, while the fully discrete proof uses the CFL-like condition

$B=B^*>0$8

(Seboldt et al., 2020).

For the linear acoustic wave equation, a more recent non-iterative time integrator combines a local explicit prediction step on interfaces with independent implicit Crank–Nicolson solves on overlapping subdomains, followed by nodal averaging (Buchholz et al., 25 Jul 2025). In first-order form,

$B=B^*>0$9

the interface values on artificial boundaries $A=A^*>0$0 are predicted by one explicit leapfrog step, local Crank–Nicolson problems are solved on each $A=A^*>0$1, and the results are averaged back to a globally continuous discrete function (Buchholz et al., 25 Jul 2025). Under the prediction CFL bound

$A=A^*>0$2

the paper proves second-order accuracy in time and the final error estimate

$A=A^*>0$3

The same work reports that the allowed $A=A^*>0$4 grows approximately linearly in the number of overlap layers $A=A^*>0$5 and that the method permits much larger time steps than leapfrog while retaining the accuracy of Crank–Nicolson (Buchholz et al., 25 Jul 2025).

A discontinuous Galerkin variant extends this wave-equation strategy to higher-order approximations and heterogeneous material parameters. It combines a local leapfrog prediction on a narrow interface strip with a local Crank–Nicolson solve on each overlapping subdomain and then local exchange in overlap regions, avoiding global averaging and remaining non-iterative in the coupling step (Buchholz et al., 31 Oct 2025). In the reported prism test, the wall time per step is $A=A^*>0$6 s for DS with Cholesky and $A=A^*>0$7 s for CN, with total simulation times $A=A^*>0$8 s and $A=A^*>0$9 s respectively, while the final-time relative $$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$0 errors are essentially the same (Buchholz et al., 31 Oct 2025). This suggests that the non-iterative template is compatible with higher-order DG discretizations, not only with linear finite elements and mass lumping.

6. Limitations, adjacent methods, and terminological boundaries

A recurring misconception is that any domain decomposition method with local solves is “non-iterative.” The cited literature repeatedly rejects that equivalence. BDD and FETI remain iterative Krylov solvers for interface problems; what changes in the error-bound paper is not the solver type but the availability of a strict upper bound separating algebraic and discretization error,

$$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$1

which yields an objective stopping criterion for the iterative solver (Rey et al., 2013). This is a non-iterative error decomposition, not a non-iterative integrator.

The same boundary applies in nonlocal problems. A substructuring-based framework for nonlocal operators proves exact equivalence between the single-domain problem and a multi-domain reformulation, both continuously and after finite element discretization, but explicitly states that the subdomain equations remain coupled and that the actual parallel solution methods are deferred to later work (Capodaglio et al., 2020). A later nonlocal FETI-type solver performs non-iterative elimination of interior variables by Schur complements and non-iterative reconstruction by local back substitution, but the reduced interface problem is solved iteratively by projected preconditioned CG, so the overall method is hybrid rather than fully non-iterative (Klar et al., 2023).

There are also non-iterative domain decomposition methods that are not time integrators. For the Helmholtz equation, the Method of Difference Potentials reduces each subdomain problem to a Calderón boundary equation with projection, assembles a single global linear system for Dirichlet and Neumann traces, solves that system once, and reconstructs the subdomain fields independently (North et al., 2021). The method is explicitly non-iterative in its interface coupling, but it is a stationary boundary solver rather than an evolutionary integrator.

Even within the true time-integration setting, the absence of inner iterations does not imply decomposition-independent accuracy. The parabolic and substructuring papers show that overlap width strongly affects the error, that increasing the number of subdomains can worsen accuracy, and that the schemes can be conditionally convergent with respect to decomposition parameters (Vabishchevich et al., 2014, Vabishchevich, 2011, Vabishchevich, 2017). In thick-structure FSI, the combination parameter $$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$2 must be tuned problem by problem; moderate values such as $$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$3 or $$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$4 behave best in the reported tests, whereas very large $$\frac{y^{n+1}-y^n}{\tau}+A y^{n+\sigma}=\varphi^{n+\sigma},\qquad y^{n+\sigma}=\sigma y^{n+1}+(1-\sigma)y^n,$$5 can degrade convergence rate (Seboldt et al., 2020). In wave propagation, the CFL-type restriction in the convergence proof depends on the overlap width, so larger overlap permits larger stable time steps (Buchholz et al., 25 Jul 2025).

Taken together, these results define a precise technical meaning for the subject. A non-iterative domain decomposition integrator is not merely a decomposed solver, nor merely a substructured formulation. It is a time-marching algorithm in which operator decomposition, overlap or interface correction, and factorized or partitioned updates are organized so that one advances the solution by a prescribed local sequence, with stability and convergence controlled analytically rather than by an inner interface iteration.