PaperlayerView exact revision on arXiv

Source-linked AI summary

An efficient coupled block preconditioner for the Marker-and-Cell discretization of the Stokes-Darcy equations

Zhen Chao, Tianwen Li

arXiv:2610.34692v1math.NAcs.NA

TL;DR

The Stokes–Darcy system yields a poorly conditioned double saddle-point problem, while existing preconditioners rely on approximating a difficult second Schur complement. The paper proposes a coupled, regularized preconditioner that avoids that approximation and analyzes its spectrum. The analysis supplies a scaling rule for the regularization parameter and characterizes spectral clustering for exact and inexact versions.

  • Problem

    Existing block triangular preconditioners require approximating two Schur complements, and the second is difficult to approximate uniformly as physical parameters vary.

    Source: §1 · p. 2

  • Method

    The proposed coupled preconditioner retains interface coupling and regularizes the constraint block with αI, avoiding any approximation of S2.

    Source: §1 · p. 2, §4.2 · p. 6

  • Results

    For the ideal preconditioner, the spectrum includes θ_i/(θ_i + α), which cluster near 1 at rate O(α); for the inexact version, eigenvalues lie within O(∥E∥^1/2) of the exact-spectrum set.

    Source: §1 · p. 2

  • Takeaways & Limitations

    Choosing α according to α/θ_min, with α = α0 min{1, κ}, keeps spectral clustering tight across the tested parameter range.

    Source: §1 · p. 2

    Replacing the direct solve for Wα with an inexact multigrid or augmented-Lagrangian solve remains open, and would remove the setup bottleneck.

    Source: §8 · p. 26

Abstract

from arXiv · show

We consider the iterative solution of the large, sparse and nonsymmetric double saddle-point systems arising from the stationary Stokes-Darcy equations in two dimensions, discretized by the Marker-and- Cell (MAC) finite difference method. Existing block triangular preconditioners rest on approximations of two Schur complements, the second of which is hard to approximate uniformly in the physical parameters. We propose instead a coupled preconditioner that keeps the interface coupling and regularizes the constraint block by a small multiple of the identity, so that the second Schur complement never has to be approximated. We analyze the spectrum of the preconditioned matrix, both in the ideal case and for the inexact version used in computations, and derive from it a rule for scaling the regularization parameter with the physical parameters. Numerical results show mesh-independent convergence of preconditioned GMRES and robustness with respect to the viscosity and the permeability.

1 Introduction

The MAC-discretized Stokes–Darcy equations produce a nonsymmetric, indefinite double saddle-point system that becomes poorly conditioned as the mesh is refined and viscosity and permeability separate. The paper adapts a coupled preconditioning idea to this system, avoiding approximation of the difficult second Schur complement, and analyzes the resulting spectrum and parameter scaling.

1 Introduction

  • The system is nonsymmetric, indefinite, and badly conditioned as h →0 and as viscosity ν and permeability κ separate.The Darcy potential, Stokes velocity, and Stokes pressure are coupled through interface conditions, giving the system its 3 × 3 double saddle-point structure.Source: §1 · p. 2
  • Existing practical block triangular preconditioners require separate approximations of S1 and S2, with the dense S2 spectrum moving with ν and κ.The prior analysis found the block triangular choice performed best, but its approximation of S2 is obtained through sparsification arguments.Source: §1 · p. 2
  • The paper adapts a coupled preconditioning idea by retaining interface coupling and regularizing the constraint block, so S2 need not be approximated.The approach adapts the idea to the MAC Stokes–Darcy system, whose nonzero, nonsymmetric (2, 2) block was outside the earlier setting.Source: §1 · p. 2
  • For the ideal preconditioner, the spectrum contains eigenvalues 1 and −1 plus values θ_i/(θ_i + α), with the third group clustering near 1 at rate O(α).When the symmetric part of S1 is positive definite, all eigenvalues lie in the disk |z − 1| < 1/2.Source: §1 · p. 2
  • For the inexact preconditioner, eigenvalues are controlled by a scalar quadratic and lie within O(∥E∥^1/2) of the exact-spectrum set, where E is supported only on interface unknowns.The regularization rule uses α/θ_min; α = α0 min{1, κ}, with α0 = 10−2, keeps clustering tight over the tested parameter range.Source: §1 · pp. 2–3

2 The Stokes–Darcy problem

The model couples free Stokes flow and porous-medium Darcy flow across an interface. Its two subdomains have different governing operators and scales that separate as viscosity and permeability vary.

2 The Stokes–Darcy problem

  • The domain is split into a free-flow region Ωs and a porous region Ωd, joined at interface Γ.The unit normal on Γ points from Ωd into Ωs.Source: §2 · p. 3, Figure 1 · p. 3
  • In the porous region, Darcy’s law and mass balance yield an equation for the Darcy pressure, denoted the Darcy potential ϕ.Source: §2 · p. 3
  • The free-flow region is governed by the Stokes equations, with velocity u, pressure p, and kinematic viscosity ν; the permeability κ is constant, assuming a homogeneous and isotropic porous medium.Source: §2 · p. 3
  • Three interface conditions enforce mass conservation, balance normal forces while allowing a pressure jump, and impose Beavers–Joseph–Saffman slip.Dirichlet conditions on the remaining boundaries make the discrete Darcy operator nonsingular.Source: §2 · pp. 3–4
  • The model’s analytical and computational difficulty comes from coupling subdomains governed by different operator types whose scales separate as ν and κ vary.Source: §2 · p. 4

3 The MAC discretization and the algebraic system

The MAC scheme places velocity and pressure unknowns on a staggered grid and discretizes the interface conditions to produce a structured algebraic system. The resulting block structure includes a low-rank interface coupling and a dense second Schur complement that the proposed preconditioner avoids forming.

3 The MAC discretization and the algebraic system

  • The MAC grid places Stokes velocity components and pressure, and Darcy potential, at distinct staggered locations in the two unit-square subdomains.Source: §3 · p. 4
  • Away from the interface, the equations use the usual five-point discretizations; at Γ, the scheme follows and eliminates ghost values to form interface equations.Eliminating the Darcy-side interface value changes the adjacent stencil’s diagonal from 4κ/h² to 3κ/h² and leaves a coupling term h−1v.Source: §3 · p. 4
  • The resulting algebraic system has structured blocks, including positive-eigenvalue As and a full-row-rank B.Source: §3 · pp. 4–5
  • The interface coupling matrix G has rank n, and its contribution through A_d^-1 is nonzero only in an n × n block.Source: §3 · p. 5
  • S2 = BS1^-1B^T is a dense n² × n² matrix that is never formed, and the preconditioner exploits the interface structure while sidestepping its approximation.Source: §3 · p. 5

4 Block preconditioning

The section contrasts ideal block preconditioners that approximate both Schur complements with a coupled alternative that retains interface blocks and regularizes the constraint block. The coupled design avoids S2, while some spectral guarantees depend on a definiteness condition that is not satisfied for every tested parameter pair.

4.1 Ideal block preconditioners

  • In practice, the best of the three ideal preconditioners loses ground as κ decreases because its implementation approximates both S1 and S2.The ideal M3 requires exact solves with S1 and S2, while the practical version uses an incomplete-Cholesky-based S1 approximation and a diagonal S2 surrogate.Source: §4.1 · p. 6

4.2 A coupled preconditioner with scalar regularization

  • The coupled preconditioner keeps G^T, B, and ±B^T, and replaces the zero constraint block with αI rather than approximating S2.Compared with K, it drops the (2, 1) block, changes the (2, 2) block to S1, reverses the (2, 3) block’s sign, and adds αI.Source: §4.2 · p. 6
  • The coupled construction does not involve S2, and Qα is nonsingular when −α is not in σ(S2); if H(S1) is positive definite, it is nonsingular for every α > 0.Source: §4.2 · p. 6
  • When H(S1) is positive definite, every eigenvalue θ of S2 has positive real part.Source: §4.2 · p. 6
  • The condition H(S1) ≻ 0 fails at κ = ν = 1 on grids with n ≥32, although the clustering result does not require that hypothesis.Source: §4.2 · p. 7

5 Spectral analysis

The spectral analysis characterizes the exact preconditioner through a Möbius map of the second Schur complement’s spectrum, then bounds localization and the effects of an inexact interface-block approximation.

5.1 The exact preconditioner

  • The exact preconditioned spectrum contains eigenvalues 1 and −1, while its remaining eigenvalues are θ/(θ + α), the Möbius images of the second Schur complement’s eigenvalues.This characterization assumes the stated nonsingularity conditions, including that −α is not in σ(S2).Source: §5.1 · pp. 7–8

5.2 Localization and clustering

  • For small α, the remaining eigenvalues approach 1 at rate O(α), yielding clusters at −1 and 1; when H(S1) ≻ 0, they also lie in the disk centered at 1/2 with radius 1/2 for every α > 0.The clustering bound assumes S2 is nonsingular and α ≤ θmin/2; the disk result uses H(S1) ≻ 0.Source: §5.2 · pp. 8–9

5.3 The inexact preconditioner

  • For every eigenvalue λ other than 0 and 1, a unit vector y gives a scalar quadratic with Rayleigh-quotient coefficients: (ŝ + γ)λ² + (ŝ + ε)λ − γ = 0.If S2 is nonsingular, zero is excluded as an eigenvalue.Source: §5.3 · p. 10

6 Practical implementation

The implementation approximates only the first Schur complement, applies the coupled preconditioner by block solves, and scales its regularization using the physical-parameter dependence of θmin.

6.1 Approximation of the first Schur complement

  • An incomplete Cholesky factorization in the natural ordering approximates the first Schur complement, avoiding the exact factorization’s O(N^2) work and O(n^3) storage.The approximation is used to construct the interface block; only an n × n triangular factor is inverted at negligible cost.Source: §6.1 · p. 12

6.2 Applying the preconditioner

  • Block elimination applies the preconditioner with solves involving Ad and Wα; the delicate second Schur complement disappears, while cond(Wα) grows like α^-1.Each GMRES iteration uses two solves with stored factorizations and sparse matrix–vector products.Source: §6.2 · pp. 12–13

6.3 Choice of the regularization parameter

  • The clustering scale is α/θmin, with θmin essentially mesh-independent and approximately min{c1(ν), c2κ}, where c2 ≈ 8.A fixed α loses effect as κ approaches zero; reducing α improves clustering but eventually worsens inner-solve accuracy as cond(Wα) grows.Source: §6.3 · p. 13

7 Numerical experiments

The experiments verify the MAC discretization and reproduce reported unpreconditioned spectral behavior across parameter settings.

7 Numerical experiments

  • The experiments use restarted, left-preconditioned GMRES with explicit stopping tests; default runs use GMRES(20), tol = 10−8, and a 500-iteration cap.The reported stopping test and iteration counts are stated to be unambiguous and exactly reproducible.Source: §7 · p. 13

7.1 Test problems

  • Example 18 verifies the discretization at fixed ν = αBJ = 1; Example 19 varies ν and κ, and Table 1 lists system sizes.Example 19 uses αBJ = 1 throughout.Source: §7.1 · p. 14, Table 1 · p. 14

7.2 Verification of the discretization

  • Example 18’s velocity components and Darcy potential converge with order two, while all four Example 19 components converge with order one.For Example 18, the pressure order rises from 1.714 on the coarsest grid pair to 1.851 on the finest; the order is still climbing at n = 512.Source: §7.2 · p. 14, Table 2 · p. 15, Table 3 · p. 15

7.3 The unpreconditioned spectrum

  • For n = 32, unpreconditioned spectra reproduce ’s qualitative pattern: imaginary parts stay below 2.5, while real parts reach ±104.For ν = 10−4 and κ = 10−8, the spectrum is entirely real and nearly symmetric about the origin.Source: §7.3 · p. 14, Figure 2 · p. 15, Table 4 · p. 16

7.4 Verification of the spectral analysis

The ideal-preconditioner spectra agree with the stated theorems, and the experiments check the assumptions and quantities used in the spectral analysis.

7.4 Verification of the spectral analysis

  • For M1, the observed multiplicities attain the theorem’s bounds, and both eigenvalue-count bounds are attained.The 15 eigenvalues left unlocated by [16, Rem. 4.2] are small perturbations of −1 and (−1 − √5)/2.Source: §7.4 · p. 16, Table 5 · p. 16
  • M2 has four distinct eigenvalues with multiplicities matching Theorem 3, while M3 has only eigenvalue 1 and Qα has 64 further eigenvalues clustered near 1.For M2, the eigenvalues are 1, −1, and (−1 ± √5)/2, with multiplicities 64, 56, 64, and 64; for M3, eigenvalue 1 has multiplicity 248.Source: §7.4 · p. 16, Table 5 · p. 16
  • M3’s sole eigenvalue does not imply one-step GMRES: its preconditioned matrix is nondiagonalizable and GMRES needs three iterations.The matrix is the unit upper triangular factor U, whose minimal polynomial is (z − 1)^3; forming M3 requires exact solves with both Schur complements.Source: §7.4 · p. 16
  • Although H(S1) is not positive definite for κ = ν = 1 on grids with n ≥ 32, Re θi > 0 in every tested case and θmin stays bounded away from zero.For κ < 1, the positive definite interface block restores H(S1) ≻ 0; the clustering result does not require that hypothesis.Source: §7.4 · pp. 16–17, Table 6 · p. 17
  • θmin is the only quantity entering the clustering bound and remains bounded away from zero on the grids considered.The observed parameter dependence of θmin motivates the scaling rule discussed in the paper.Source: §7.4 · pp. 16–17

7.5 The role of α

The regularization parameter α controls spectral clustering, while the coupled preconditioner bQα delivers mesh-independent GMRES counts and the lowest total CPU time among the four practical preconditioners in the reported comparison.

7.5 The role of α

  • The count falls from about 40 at α = 10 to 6–8 at α = 10^-3, while for each fixed α it grows only mildly with mesh size.For α ≤ 10^-2, counts rise by at most two iterations over the full grid range.Source: §7.5 · p. 17, Table 7 · p. 18
  • The sharper eigenvalue bound overestimates the measured maximum distance from 1 by at most a factor of 1.4 and becomes asymptotically exact as α decreases.It is correct to two digits at α = 10^-3 and to three digits from α = 10^-4 on.Source: §7.5 · p. 17, Figure 4 · p. 19, Table 8 · p. 19

7.6 Comparison of the preconditioners

  • bQα needs 8 iterations at n = 16 and 9 at every n from 32 to 1024, versus 16–22 for c M3 and 52–64 for c M2.The largest system has more than 4 × 10^6 unknowns; bQα's total CPU time is also smallest at every n.Source: §7.6 · p. 20, Table 9 · p. 20
  • At n = 64, the bQα residual curve decreases more rapidly than those of the other three practical preconditioners throughout the iteration.The plotted runs use a stopping tolerance of 10^-8.Source: §7.6 · p. 20, Figure 5 · p. 21
  • bQα is the cheapest of the four during iteration, by factors of 1.2–1.5 relative to c M3 and 2.8–4.5 relative to c M1.Its total-time advantage over c M3 is smaller on the finest meshes because setup accounts for 71% of bQα's total time at n = 1024.Source: §7.6 · p. 21
  • At n = 256, c M1 and c M2 take 2.0 and 2.1 times bQα's total time, respectively.The setup times are nearly identical for n ≥ 64, differing by at most 4–7%.Source: §7.6 · pp. 20–21

7.7 Quality of the Schur complement approximation

The interface-block approximation is central to the preconditioner’s performance. With scaled α, bQα converges across the tested permeability range, while the computed solution resolves the stated interface behavior.

7.7 Quality of the Schur complement approximation

  • The exact interface block T gives iteration counts that barely depend on κ, whereas the surrogate T ≈ (3κ)^-1 I fails at κ = 10^-6 on the n = 128 grid.Incomplete Cholesky approximations degrade more as the drop tolerance is loosened, especially at small κ, where ∥T∥ = O(κ^-1).Source: §7.7 · p. 22, Table 11 · p. 22

7.8 Robustness with respect to the physical parameters

  • Across 72 parameter cases, bQα with scaled α converges in all cases, while c M3 stagnates in seven, all at small κ on finer grids.For κ ≥ 10^-2, bQα requires 8–19 iterations versus 12–36 for c M3 at every tested mesh size and both tested viscosities.Source: §7.8 · p. 22, Table 12 · p. 23, Table 13 · p. 23
  • At n = 256, ν = 1, κ = 10^-6, scaling α reduces bQα's count from 102 to 27; all table entries stay at or below 34.The largest counts occur at ν = 10^-2, n = 256, κ ≤ 10^-7, where α ≤ 10^-9 and conditioning of Wα begins to affect performance.Source: §7.8 · p. 22, Table 12 · p. 23, Table 13 · p. 23
  • For ν = 10^-2 and n = 256, c M3 converges in 14 iterations for every κ ≤ 10^-4, but at ν = 1 it needs 264 iterations at κ = 10^-5 and fails for κ ≤ 10^-6.The comparison is summarized on the finest grid in Figure 6.Source: §7.8 · p. 22, Figure 6 · p. 24

7.9 A computed solution

  • The computed solution resolves the permitted jump across Γ and shows no visible kink in the normal velocity component or error layer along the interface.The pointwise error is largest in the Stokes region, consistent with the error norms of Table 2.Source: §7.9 · p. 24, Figure 7 · p. 25

8 Concluding remarks

The paper introduces a coupled preconditioner that replaces approximation of the second Schur complement with scalar regularization. Its spectral analysis supports parameter-scaled regularization, and the reported experiments show mesh-independent iteration counts and robustness across tested parameters.

8 Concluding remarks

  • Qα replaces the second Schur complement with scalar regularization, leaving only the structurally known first Schur complement to approximate.This distinguishes it from the block triangular preconditioners discussed in the paper.Source: §8 · p. 26
  • The ideal spectrum has eigenvalues 1 and −1 with multiplicities n^2 and 2n^2 − n, while the remaining n^2 eigenvalues are θ/(θ + α) images of the second Schur-complement spectrum.The analysis also gives an O(α/θmin) clustering rate and a rule for scaling α; for the inexact variant, the spectral displacement is bounded by O(∥E∥^1/2), with E confined to interface unknowns.Source: §8 · p. 26
  • At unit parameters, the preconditioner cuts GMRES iterations by factors of two to seven relative to block triangular preconditioners with no mesh dependence, and scaled α converges in all 72 tested configurations.The exception is ν = 10^-2 on the finest grid, where c M3 uses fewer iterations in cases where it converges.Source: §8 · p. 26
  • Replacing the direct Wα solve with an inexact solve remains open; it would remove the setup bottleneck, but E would no longer be confined to the interface and Lemma 15 would be lost.The paper also identifies extensions to three dimensions, non-uniform grids, and time-dependent Navier–Stokes/Darcy coupling as open questions.Source: §8 · p. 26

Use of Artificial Intelligence

The authors report using Claude to polish the manuscript’s English and check consistency, while retaining responsibility for its scientific interpretation and final form.

Use of Artificial Intelligence

  • Claude (Anthropic) was used to polish the manuscript’s English and check consistency; the authors state that they are responsible for the scientific interpretation and final manuscript.Source: Use of Artificial Intelligence · p. 26
Loading 2610.34692v1…