SLIDE 1
✫ ✬ ✪ ✩
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Outline
- The coupled scattering problem
- Discretization by Finite Elements
- Multigrid Method
– Multigrid algorithms – Computational results in 2D
- Subsctructuring Method
Iterative Solvers for Coupled Fluid-Solid Scattering Jan Mandel - - PowerPoint PPT Presentation
Iterative Solvers for Coupled Fluid-Solid Scattering Jan Mandel - - PowerPoint PPT Presentation
Iterative Solvers for Coupled Fluid-Solid Scattering Jan Mandel Center for Computational Mathematics Department of Mathematics University of Colorado at Denver Center for Aerospace Structures Department of Aerospace Engineering
SLIDE 2
SLIDE 3
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Main results
- treating a coupled problem with vastly different scales
⇒ algorigthms are invariant to scaling of physical units ⇒ physical units must match ⇒ define residual with care
- Multigrid on a properly scaled coupled problem with GMRES as
smoother converges well if coarse problem fine enough
- extended
FETI-H (De La Bourdonnaye, Farhat, Macedo, Magoul` es, Roux, 1998, 2000; Farhat, Macedo, Tezaur, DD11, 1999) to coupled problem ⇒ coupled algorithm reduces to FETI- H in the limit for a very stiff obstacle
SLIDE 4
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Coupled problem Ωf Γn Γn Γd Γa Γ Ωe
✲ ν
- solid obstacle in fluid
- time harmonic fluid pressure p(x, t) = p(x)eiωt
= ⇒ Helmholtz equation in the fluid
- time harmonic solid displacement u(x, t) = u(x)eiωt
= ⇒ elastodynamic Lam´ e equation in the solid
- continuity of displacement & balance of normal forces
= ⇒ wet interface conditions
SLIDE 5
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Coupled Problem Fluid Medium -the Helmholtz Equation
- isotropic, homogeneous, inviscid
- compressible, irrotational
- time-harmonic solution of the wave equation
- complex amplitude p(x): the fluid pressure at location x and time t is
p(x)eiωt −∆p − k2p = f in Ω p = p0
- n Γd
∂p ∂n
=
- n Γn
∂p ∂n + ikp
=
- n Γa (radiation b.c.)
- k = ω
cf, ω is wave frequency, cf is speed of sound in acoustic medium
- radiation boundary condition - Sommerfeld
– does not reflect waves in the normal direction – first-order approximation for radiation to infinity
SLIDE 6
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Coupled Problem Elastic Medium
- isotropic, homogeneous
- small time-harmonic displacement u(x)eiωt of an elastic body
- complex displacement amplitude u(x) satisfies the Lam´
e equations ∇ · τ + ω2ρe u = 0 in Ωe + boundary conditions, where : ρe − density of the elastic medium
τ
= λI(∇ · u) + 2µe(u) stress tensor e(u) = 1 2(∇u + (∇u)T) strain tensor λ, µ − Lam´ e coefficients of elastic medium
SLIDE 7
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Coupled Problem Fluid-Solid Interface
- 1. continuity across interface
n · u =
1 ρfω2 ∂p ∂n
- 2. balance of normal forces:
n · τ · n = −p
- 3. nonviscous fluid → zero tangential tension:
n × τ · n = 0
SLIDE 8
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Coupled Problem Summary of the Pressure-Displacement Formulation Helmholtz equation in the fluid region: ∆p + k2p = 0 in Ωf, p = p0 on Γd, ∂p ∂n = 0 on Γn, ∂p ∂n + ikp = 0 on Γa. Lam´ e equation in the elastic region: ∇ · τ + ω2ρeu = 0 in Ωe,
τ = λI(∇ · u) + 2µe(u),
eij(u) = 1 2(∂ui ∂uj + ∂uj ∂xi ), Wet interface conditions:
n · u =
1 ρfω2 ∂p ∂n
n · τ · n = −p n × τ · n = 0 on Γ
On Γ, the value of u provides load for the Helmholtz problem for p and the value of p provides load for the elastodynamic problem for u.
SLIDE 9
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Existence of solution and non-radiating modes Solution exists (p, u) ∈ H1(Ωf) × H1(Ωe)3 and is unique up to non-radiating modes in the solid
- the boundary value problem
∇ · τ + ω2ρeu = in Ωe
τ · n
=
- n
Γ
u · n
=
- n
Γ has nonzero solution for certain frequencies and geometries
- called non-radiating modes.
- for this to happen, ω2ρe needs to be an eigenvalue of the pure traction
problem and in addition u · n = 0.
- bodies
with certain symmetries have non-radiating modes (e.q., sphere can sustain torsional oscillations with zero radial component displacement)
- almost all elastic bodies do not have non-radiating modes [Harg´
e]
SLIDE 10
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Variational Formulation Find (p, pΓ), (p − p0, pΓ) ∈ Vf and (u, n · uΓ) ∈ Ve such that −
- Ωf
∇p∇¯
q + k2
- Ωf
p¯ q − ik
- Γa
pΓ¯ qΓ −
- Γ
ρfω2(n · uΓ)¯ qΓ = 0 −
- Ωe
(λ(∇ · u)(∇ · ¯
v) + 2µe(u) : e(¯ v)) + ω2
- Ωe
ρeu · ¯
v −
- Γ
pΓ(n · ¯
vΓ) = 0
∀ ¯ q, ¯ qΓ ∈ Vf and ∀ ¯
v, n · ¯ vΓ ∈ Ve,
where Vf = { (p, pΓ) | p ∈ H1(Ωf), pΓ ∈ H
1 2(Γ) |
p = 0
- n
Γd} Ve = { (u, n · uΓ) |
u ∈ (H1(Ωe))3, uΓ ∈ (H
1 2(Γ))3},
Same wet interface term in both equations = ⇒ symmetric system Industry standard (Morand and Ohayon 1995)
SLIDE 11
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Discretized system Finite elements = ⇒ 2 × 2 symmetric block system
- −Kf + k2Mf − ikGf
−ρfω2T −T∗ −Ke + ω2Me
p u
- =
- r
- T is the coupling matrix: p∗Tv =
Γ
p(ν · v)
T is transfer of load =
⇒ lower order, weak coupling
SLIDE 12
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Scaling Since λ and µ are large ⇒ use a scaling of the form u = su′ and v = sv′, where s is a scalar, to control leading terms Discretization and Error Bound Iterative Soution leading terms differ by factor of k2 comparable leading terms consistent with definition of norm scale to O(1) diagonal s = 1 c
- ρf max{λ, 2µ}
. s = 1 ck
- ρf max{λ, 2µ}
- −Sf + k2Mf + ikGf
−ρfω2sT −ρfω2sTt −ρfω2s2Se + ρeρfω4s2Me p
u
- = R
SLIDE 13
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Finite element error estimates If the solution (p, u) ∈ H1+α(Ωf) × H1+α(Ωe)3 and h is small enough then the discretization error satisfies p−ph2
H1(Ωf) +k2p−ph2 L2(Ωf) +k2(u−uh2 H1(Ωe)3 +k2u−uh2 L2(Ωe)3) ≤ C(k)h2α
Proof: G˚ arding inequality, scales of spaces, regularity, approximation. PhD thesis of M. Popa, 2002
SLIDE 14
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Model Problem Ωf Γn Γn Γd Γa Γ Ωe
✲ n
- 2D channel filled with water, unit square
- excitation on the left (Dirichet boundary condition) on Γd
- sound-hard sides(Neumann boundary condition) on Γn
- outgoing radiation on the right (absorbing boundary condition) on Γa
- aluminum scatterer in the middle
- discretization by uniform mesh, P1 elements
SLIDE 15
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Numerical Solution Square scatterer size 0.2 m in the middle
30 60 30 60 −1.5 1.5
y−axis x−axis Γd Γn Γn Γa
30 60 30 60
x−axis y−axis Γn Γn Γa Γd
Fluid pressure Solid displacement
SLIDE 16
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Numerical Solution Obstacle with a gap
30 60 30 60 −3 3
y−axis x−axis Γd Γn Γn Γ a
20 40 60 20 40 60 −2 −1 1 2
y−axis x−axis Γd Γn Γn Γ a
20 40 60 20 40 60
x−axis y−axis Γn Γn Γa Γd
20 40 60 20 40 60
x−axis y−axis Γn Γn Γa Γd
gap on x axis gap on y axis
SLIDE 17
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Multigrid for the Coupled Problem
- coarsening: standard variational bilinear interpolation, mesh ratio 2
- smoothing
– GMRES, BICG-STAB – Preconditioners ∗ inverse of lower triangular part of A ∗ inverse of nodal block diagonal – Gauss-Seidel
- for stable smoothers can increase number of steps without causing
divergence
- all smoothers work well for fluid or solid part separately
- wet interface aligned with coarsest grid, otherwise nothing special
SLIDE 18
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Multigrid convergence, scatterer size 0.2 Decreasing h, increasing k, k3h2 constant Adding coarse meshes
2 3 4 5 2e−82e−7 2e−6 2e−5 10
−0.3
10
−0.2
10
−0.1
multigrid levels GMRES as smoother h residual reduction
2 3 4 5 2e−5 2e−6 2e−7 2e−8 10
−0.17
10
−0.13
10
−0.09
multigrid levels BICG−STAB as smoother h residual reduction
2 3 4 5 2e−8 2e−7 2e−6 2e−5 10
−5
10 10
5
10
10
multigrid levels GAUSS−SEIDEL as smoother h residual reduction
2 3 4 5 2e−5 2e−6 2e−7 2e−8 10
−0.4
10
−0.3
10
−0.2
multigrid levels GMRES preconditioned by inv lower triangular part of A as smoother h residual reduction
SLIDE 19
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Multigrid convergence, scatterer size 0.2 Decreasing h, increasing k, k3h2 constant Adding smoothing steps
3 10 20 2e−6 2e−7 2e−8 2e−9 10e−3 10e−1
smoothing steps GMRES as smoother h residual reduction
3 10 20 2e−6 2e−7 2e−8 2e−9 10
−3
10
−2
10
−1
10
smoothing steps BICG−STAB as smoother h residual reduction
3 10 20 2e−5 2e−6 2e−7 2e−8 10
−4
10
−3
10
−2
10
−1
GMRES preconditioned by inv D as smoother smoothing steps h residual reduction
3 10 20 2e−6 2e−7 2e−8 2e−9 10
−4
10
−3
10
−2
10
−1
GMRES precondition by inv lower triangular part of A as smoother smoothing steps h residual reduction
SLIDE 20
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Multigrid convergence, scatterer size 0.2 GMRES preconditioned by Multigrid Multigrid with different smoothers
1/256 1/128 1/64 10
−3
10
−2
10
−1
h GMRES preconditioned by multigrid residual reduction BICG−STAB GMRES prec by inv lower triang A GMRES
1/256 1/128 1/64 10
−4
10
−2
10
Multigrid residual reduction BICG−STAB GMRES prec by inv lower triang A GMRES
Decreasing h, k3h2 constant
3 10 20 10
−3
10
−2
10
−1
soothing steps GMRES preconditioned by multigrid residual reduction BICG−STAB GMRES prec by inv lower triang A GMRES
3 10 20 10
−3
10
−2
10
−1
10
smoothing steps Multigrid residual reduction BICG−STAB GMRES prec by inv lower triang A GMRES
h = 1/128, k = 25
SLIDE 21
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Multigrid convergence, scatterer with a gap GMRES preconditioned by Multigrid and Multigrid with smoother GMRES
1/512 1/256 1/128 1/64 1/32 10
−0.8
10
−0.1
h Obstacle 0.2 in x−dir 0.4 in y −dir gap on x−axis of size 40% on x 50% on y residual reduction GMRES prec Multigrid Multigrid smoother GMRES
Decreasing h, k3h2 constant
SLIDE 22
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Summary of multigrid performance
- efficient multigrid method for the discrete coupled system
- Krylov subspace smoothers, in particular GMRES, make robust algorithms when used
with multigrid.
- large number of smoothing steps - GMRES preconditioned by multigrid gives comparable
results to multigrid with smoother GMRES
- more smoothing steps have a stabilizing effect on the solution.
- adding coarse meshes - Gauss Seidel as smoother diverges, Krylov smoothers are robust
- multigrid methods are found to diverge, or give very poor convergence for high wave
numbers, totally different methods are required for this purpose.
- in some cases, multigrid with GMRES as smoother gives slightly better convergence
rates than GMRES preconditioned by one multigrid cycle, otherwise about same.
- multigrid with smoother GMRES and multigrid with smoother GMRES preconditioned
by inverse of lower triangular part of A were most robust and gave best results, however they require the most memory to store all fine grid vectors.
SLIDE 23
✫ ✬ ✪ ✩
University of Colorado Jan Mandel GMRES with 2 × 2 block diagonal ILU precondioning scatterer size 0.2 in the middle
- in fluid let Lf · Rf be approximate factorization of −Kf + k2Mf − ikGf
- in solid let Le · Re be approximate factorization of −ρfω2s2Ke + ρfρeω4s2Me
GMRES preconditioned by the block diagonal matrix
- (Lf · Rf)−1
(Le · Re)−1
- For some frequencies, the matrix −Ke + ω2Me will be singular - scatterer is at
resonance.
1/64 1/128 1/200 10
−6
10
−4
10
−2
10
h GMRES preconditioned by standard ILU to 2 × 2 block diagonal relative residual natural frequency k3h2=const
SLIDE 24
✫ ✬ ✪ ✩
University of Colorado Jan Mandel FETI-H for coupled problem
- FETI-H for Helmholtz equation (De La Bourdonnaye, Farhat, Macedo,
Magoul` es, F. Roux 1998; Farhat, Macedo Lesoinne 2000)
- Extension to solid
- Coupled system
SLIDE 25
✫ ✬ ✪ ✩
University of Colorado Jan Mandel FETI-H for two subdomains (continuous form) Ω1 Ω2 Γ12 −∆p1 − k2p1 = in Ω1 + b.c. −∆p2 − k2p2 = in Ω2 + b.c. p1 = p2
- n Γ12
∂p1 ∂ν1
= −∂p2
∂ν2
- n Γ12
FETI idea: enforce p1 = p2 by Lagrange multiplier, eliminate p1, p2, iterate on the system for the Lagrange multiplier. But here p1 p2 may not be determined due to (near) resonance. Hence ∂p1
∂ν1 = −∂p2 ∂ν2 replaced by a linear combination ∂p1 ∂ν1 − ikp1
= −∂p2
∂ν2 − ikp2
- n Γ12
SLIDE 26
✫ ✬ ✪ ✩
University of Colorado Jan Mandel FETI-H for two subdomains (discrete form) Discrete equations: Kp = f, where K = S − k2M + ikMS Decomposed, continuity p1 = p2
- n Γ12
to be enforced by Lagrange multipliers: subassembly in Ω1 :
K1p1
+
B1Tλ
=
f1
subassembly in Ω2 :
K2p2
+
B2Tλ
=
f2
continuity across Γ12 :
B1p1
+
B2p2
= After replacement of interface condition by a linear combination becomes (K1 + R1)p1 =
f1 − B1Tλ
(K2 + R2)p2 =
f2 − B2Tλ B1p1
+
B2p2
= Upon assembly the contributions of R1 and R2 to the global system would cancel = ⇒ equivalent system The artificial radiation condition gives R1 = ikM12, R2 = −ikM12 (the interface mass matrix
- r
its diagonal approximation, with natural dof mapping)
SLIDE 27
✫ ✬ ✪ ✩
University of Colorado Jan Mandel FETI-H: multiple subdomains
- for each subdomain decide sign of the artificial radiation condition
- for each interface between subdomains, add artificial radiation matrix to
subdomain matrix if not in conflict between neighbors
- subdomain problems guaranteed solvable (Farhat 1998)
- eliminate original dofs, get system for Lagrange multipliers
FIλ = d
- reduced problems has the size of the interface
- symmetric, but not Hermitian
− → GCR(GMRES)
- Numerical experiments indicate scalability with respect to h
SLIDE 28
✫ ✬ ✪ ✩
University of Colorado Jan Mandel FETI-H: Coarse space Enforce an optional constraint QTrk = 0 for the residuals rk = d − FIλk by using the splitting
λ = P¯ λ + λ0
The original interface problem is transformed into the modified problem
PTFIPλ = FIPλ = PTd
where
P = I − Q(QTFIQ)−1QTFI
and
λ0 = Q(QTFIQ)−1QTd
SLIDE 29
✫ ✬ ✪ ✩
University of Colorado Jan Mandel FETI-H method: Construction of the coarse space The coarse space basis is
Q =
- B1Q1
· · ·
BsQs
· · ·
BNsQNs
where Qs
j is the discrete representation of the plane wave eikdT
j x, with the
directions dj = [cos θj, sin θj], where θj = (j − 1)2π
Nθ,
j = 1, . . . , Nθ :
SLIDE 30
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Extension of FETI-H to elastic scattering, two subdomains Ω1 Ω2 Γ12 ∇ · τ(u1) + ω2ρeu1 = in Ω1 + b.c. ∇ · τ(u2) + ω2ρeu2 = in Ω2 + b.c. u1 = u2
- n Γ12
continuity of displacement τ(u1) · ν1 = −τ(u2) · ν2
- n Γ12
continuity of traction FETI idea: enforce u1 = u2 by Lagrange multipliers, eliminate u1, u2, iterate
- n the system for the Lagrange multiplier.
But here u1 u2 may not be determined due to (near) resonance. Hence τ(u1) · ν1 = −τ(u2) · ν2 replaced by linear combination: τ(u1) · ν1 + iωaν1(ν1 · u1) = −τ(u2) · ν2 + iωaν2(ν2 · u2) on Γ12.
SLIDE 31
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Elastic scattering, two subdomains, discrete form Discrete equations: Ku = f where K = S − ω2M Decomposed and u1 = u2 on Γ12 enforced by Lagrange multipliers: subassembly in Ω1 :
K1u1
+
B1Tλ
=
f1
subassembly in Ω2 :
K2u2
+
B2Tλ
=
f2
continuity across Γ12 :
B1u1
+
B2p2
= Equivalent to: (K1 + R1)u1 =
f1 − B1Tλ
(K2 + R2)u2 =
f2 − B2Tλ B1u1
+
B2u2
= Upon assembly the contributions of R1 and R2 to the global system would cancel = ⇒ equivalent system
SLIDE 32
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Substructuring for the coupled problem: Overview
- FETI-H:
– decomposition of solid and fluid domains – artificial radiation condition to ensure solvability on subdomains – Lagrange multipliers on subdomain interfaces – eliminate primary variables = ⇒ uncoupled local problemss
- here: local problems coupled across wet interface
- solution:
DUPLICATE variables on wet interface, eliminate original primary variables, leave copies in the system iterated on
- =
⇒ numerically almost triangular system containing FETI-H for fluid and solid as diagonal blocks
Jan Mandel, Iterative Substructuring with Lagrange Multipliers for Coupled Fluid-Solid Scattering, DD14, 2002 Jan Mandel, An Iterative Substructuring Method for Coupled Fluid-Solid Acoustic Problems, JCP 2002 http://www-math.cudenver.edu/~jmandel/papers
SLIDE 33
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Model coupled problem Ωf Γn Γn Γd Γa Γ Ωe
✲ ν
- 2D channel filled with water
- excitation on the left (Dirichet boundary condition) on Γd
- sound-hard sides (Neumann boundary condition) on Γn
- outgoing radiation on the right (absorbing boundary condition) on Γa
- aluminum scatterer in the middle
SLIDE 34
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Model 2D Problem Decomposed in 5x5 Fluid and 2x2 Solid Subdomains
SLIDE 35
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Decomposition and subassembly Non-overlapping subdomains: Ωf = Nf
s=1 Ω s e,
Ωe = Ne
s=1 Ω s e.
Subdomain local vectors and matrices:
ˆ p = p1
. . .
pNf , ˆ u = u1
. . .
uNe , ˆ Kf = K1
f
. . . . . . ... . . . . . .
KNf
f
, ps∗Ks
fq =
- Ωs
f
∇p∇q, ( ˆ
Ke ˆ Mf ˆ Me defined similarly) ˆ T = T11
. . .
T1,Ne
. . . ... . . .
TNf,1
. . .
TNf,Ne , pr∗Trsvs =
- ∂Ωr
f∩∂Ωs e
p(ν · v)
SLIDE 36
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Intersubdomain continuity Local to global maps Nf and Ne:
Kf = N∗
f ˆ
KfNf, Ke = N∗
e ˆ
KeNe ˆ p = Nfp, ˆ u = Neu.
To enforce same values between subdomains:
Bf = [B1
f, . . . , BNf f ],
Be = [B1
e, . . . , BNe e ]
such that
Bfˆ p = 0
⇐ ⇒
ˆ p = Nfp
for some
p Beˆ u = 0
⇐ ⇒
ˆ u = Neu
for some
u
SLIDE 37
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Decomposed system Introduce Lagrange multipliers to enforce intersubdomain continuity:
− ˆ
Kf + k2 ˆ Mf
−ω2ρf ˆ
T B∗
f
−ˆ
T∗
− ˆ
Ke + ω2 ˆ Me B∗
e
Bf Be ˆ p ˆ u
λf λe
= ˆ r
The system for (ˆ
p, ˆ u, λf, λe) is equivalent to the original system via the local
to global maps: ˆ
p = Nfp and ˆ u = Neu.
But − ˆ
Kf + k2 ˆ Mf and − ˆ Ke + ω2 ˆ Me are typically (close to) singular due to
(near) resonance . . .
SLIDE 38
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Regularized system
ˆ Af
−ω2ρf ˆ
T B∗
f
−ˆ
T∗ ˆ Ae B∗
e
Bf Be ˆ p ˆ u
λf λe
= ˆ r
where with σst = ±1, σst = −σst,
ˆ Af
= − ˆ
Kf + k2 ˆ Mf + ˆ Rf ˆ Ae
= − ˆ
Ke + ω2 ˆ Me + ˆ Re ˆ Rf
= (Rrs
f )rs,
ps∗Rfqs = iα0k
- t=s
σst
- ∂Ωs
f∩∂Ωt f
pq,
ˆ Re
= (Rrs
e )rs,
us∗Revs = iα0ω
- (λ + 2µ)ρe
- t=s
σst
- ∂Ωs
e∩∂Ωt e
(n · u)(n · v),
ˆ N∗
fRf ˆ
Nf = 0, ˆ N∗
eRe ˆ
Ne = 0 ⇒ equivalent to the decomposed system
SLIDE 39
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Choice of the artificial radiation matrices If σst does not change sign on Ωs
f, then − ˆ
Ks
f + k2 ˆ
Ms
f + ˆ
Rr
f is regular (Farhat,
Macedo, Lesoinne 2000). Similary for solid subdomains. In computational tests, we assigned σst by counting, did not try to avoid change of sign. If α0 = ±1, the artifical intersubdomain conditions are satisfied by wave in the normal direction (pressure wave in solid). But we will see numerically that α0 = 0 is OK except when very close to resonance.
SLIDE 40
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Augmented system Key idea: ˆ
Tˆ u, ˆ T∗ˆ p depend on the values of ˆ u, ˆ p on the wet interface Γ only.
Define ˆ
Jf, ˆ Je as expanding vector on Γ by zero entries, then ˆ Tˆ u = ˆ TJeˆ uΓ, ˆ uΓ = J∗
eˆ
u, ˆ T∗ˆ p = ˆ T∗Jfˆ pΓ, ˆ pΓ = J∗
fˆ
p.
Get the augmented system
ˆ Af B∗
f
−ω2ρf ˆ
TJe ˆ Ae B∗
e
−ˆ
T∗Jf Bf Be
−J∗
f
I
−J∗
e
I ˆ p ˆ u
λf λe
ˆ pΓ ˆ uΓ
=
ˆ r
SLIDE 41
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Primal-dual method for the coupled problem
- scale the augmented system to symmetrize the coupling terms and so
that the fluid and the elastic equations have the same order of magnitude
- eliminate the subdomain degrees of freedom ˆ
p and ˆ u from the augmented
system
- solve the resulting reduced system for the intersubdomain Lagrange
multipliers λf, λe and the wet interface degrees of freedom ˆ
pΓ, ˆ uΓ by
Generalized Conjugate Residuals
- precondition by projection on a coarse space
SLIDE 42
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Scaling
- multiply the second equation by ω2ρf
- symmetric diagonal scaling
- Af
- B∗
f
−
TJe
- Ae
- B∗
e
−
T∗Jf
- Bf
- Be
−J∗
f
I
−J∗
e
I
- p
- u
- λf
- λe
- pΓ
- uΓ
=
- r
, where
- Af = Df ˆ
AfDf,
- Ae = ω2ρfDe ˆ
AeDe,
- T = ω2ρfDf ˆ
TDe,
- Bf = EfBfDf,
- Be = EeBeDe,
- r = Dfˆ
r, ˆ p = Df p, ˆ u = De u,
λf = Df λf, λe = ω2ρfDf λf.
SLIDE 43
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Reduced system Eliminating
- p
=
- A−1
f (
r − B∗
f
λf +
TJe uΓ)
- u
=
- A−1
e (−
B∗
e
λe +
T∗Jf pΓ)
gives the final reduced system
Fx = b,
where
F =
- Bf
A−1
f
- B∗
f
−
Bf A−1
f
TJe
- Be
A−1
e
- Be
−
Be A−1
e
- T∗Jf
J∗
f
A−1
f
- B∗
f
I J∗
f
A−1
f
TJe Je A−1
e
- Be
−Je
A−1
e
- T∗Jf
I ,
and
T is small for small kh or stiff scatterer. Here, x =
λf λe
- pΓ
- uΓ
, b =
- Bf
A−1
f
r
−J∗
f
A−1
f
r .
SLIDE 44
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Decoupling of fluid and elastic equations
- the first diagonal block
Bf A−1
f
- B∗
f is the FETI-H operator
- the second diagonal block
Be A−1
e
- Be is the same for elasticity
- as the material is more rigid λ, µ → ∞, the scaling matrix De → 0
⇒ the coupling term
T = ω2ρfDf ˆ TDe → 0
⇒ the fluid and elastic equations are uncoupled in the limit ⇒ iterations “should” behave much like for fluid or elasticity alone (can anything be said in general about the convergence of GMRES for a block diagonal matrix from convergence on the diagonal blocks... let alone block triangular...)
SLIDE 45
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Estimate of the coupling matrix Since De ≈ (ω2ρfµhn−2)−1/2, it follows that
- T
= ω2ρfDf ˆ
TDe ≤ ω2ρfDf ˆ T De
≈ ω2ρf(hn−2)−1/2hn−1(ω2ρfµhn−2)−1/2 = ωhρ1/2
f
µ−1/2 = khcfρ1/2
f
µ−1/2. Consequently, the system will become numerically decoupled if khcfρ1/2
f
µ−1/2 << 1
SLIDE 46
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Iterative solution Enforce the residual condition
Q∗(Fx − b) = 0
throughout the iterations. For given v use the initial approximation
x(0) = v + Qw, w obtained by solving the residual correction equation, Q∗(F(v + Qw) − b) = 0. Fx = b solved by GCR with left preconditioning by the projection P = I − Q(Q∗FQ)−1Q∗F
and initial iterate x(0). Equivalently, GCR applied to
PFx = Pb.
Iterations run in a subspace:
Q∗F(x(n) − x(0)) = 0
SLIDE 47
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Selection of coarse space
Q = DfBfYf DeBeYe DfJ∗
fZf
DfJ∗
fZe
Coarse space selection for multipliers (same/similar to FETI-H):
Yf = diag(Ys
f), columns of Ys f are discrete representations of plane waves in
a small number of equally distributed directions, or discrete representation of the constant function.
Ye = diag(Ys
e), columns of Ys e are discrete representations of elastic plane
waves (both pressure and shear) in a small number of equally distributed directions, or discrete representation of the rigid body motions Coarse space selection for wet interface (New): The matrices Zs
f and Zs e are chosen in the same way as Ys f and Ys e, with possibly
different selection of the number of directions and selection of constant or rigid body modes. Some of the matrices Ys
f, Ys e, Zs e, or Ys e may be void.
SLIDE 48
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Artificial radiation condition on wet interface
- −Kf + k2Mf − ikGf
−ρfω2T −T∗ −Ke + ω2Me
p u
- =
- r
- But if there is only one solid subdomain Ke − ω2Me may be (nearly) singular
⇒ no convergence of iterations on reduced system with u eliminated Solution: simulate radiation of the solid into the fluid (physics: solid damped by the radiation condition in the fluid via the fluid): add the first equation multiplied by iαT′ to the second equation The diagonal block becomes −Ke + ω2Me + iβρfω2T′T choose β = β0ω
- ρe(λ + 2µ)/T1 so that the added term approximates the
radiation condition satisfied by pressure waves in normal direction and physical units match ⇒ proceed as before, but now iterations converge Generalization to multiple fluid subdomains straightforward: add decomposed equations for fluid subdomains, multiplied by the appropriate submatrix of T, to the equation for the elastic domain Numerical results: loss of convergence when ω2 very close to eigenvalue, but all is fine with multiple solid subdomains
SLIDE 49
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Computational complexity Convergence improves with inreasing coarse space size but complexity increases Assume: complexity of LU decomposition of a block sparse matrix of order nm, where m is the order of dense blocks, is about nαm3, complexity of one solution is nβm2 Denote: N the total number of degrees of freedom Ns the number of subdomains Ni the total number of degrees of freedom on subdomain interfaces, Nc the number of coarse space basis vectors per subdomain Nt interfaces between the subdomains Then the dominant cost of the method is Tsetup ≈ Ns(N/Ns)α + NcNt(N/Ns)β + Nα
s N3 c + N2 c Ni,
Titeration ≈ Ns(N/Ns)β + Nβ
s N2 c .
SLIDE 50
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Numerical results: scalability with number of solid subdomains
Problem description Coarse directions Number of Subdomains Multipliers Wet interf. degrees of freedom
h k f e f e f e Orig Red Crs It 1/250 20 4x4 rigid 4 60600 1416 63 32 1/250 20 4x4 1x1 4 65802 1824 63 37 1/250 20 4x4 1x1 4 4 65802 1824 86 32 1/250 20 4x4 rigid 4 60600 1416 63 32 1/250 20 4x4 2x2 4 65802 2034 72 48 1/250 20 4x4 2x2 4 4 65802 2034 100 42 1/250 20 4x4 rigid 4 60600 1416 63 32 1/250 20 4x4 3x3 4 65802 1824 63 37 1/250 20 4x4 3x3 4 4 65802 1824 86 32 1/250 20 4x4 rigid 8 60600 1416 102 26 1/250 20 4x4 1x1 8 4 65802 1824 125 27 1/250 20 4x4 1x1 8 4 4 4 65802 1824 126 27 1/250 20 4x4 rigid 8 60600 1416 102 26 1/250 20 4x4 2x2 8 4 65802 2034 139 39 1/250 20 4x4 2x2 8 4 4 4 65802 2034 148 36 1/250 20 4x4 rigid 8 60600 1416 102 26 1/250 20 4x4 3x3 8 4 65802 1824 125 27 1/250 20 4x4 3x3 8 4 4 4 65802 1824 126 27
SLIDE 51
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Numerical results: scalability with number of fluid and solid subdomains
Problem description Coarse directions Number of Subdomains Multipliers Wet interf. degrees of freedom
h k f e f e f e Orig Red Crs It 1/315 32 3x3 rigid 4 96012 1268 32 45 1/315 32 3x3 1x1 4 104204 1776 32 56 1/315 32 3x3 1x1 4 4 104204 1776 40 51 1/315 32 5x5 rigid 4 96012 2292 96 45 1/315 32 5x5 3x3 4 104204 2804 96 57 1/315 32 5x5 3x3 4 4 104204 2804 111 46 1/315 32 3x3 rigid 8 96012 1268 52 27 1/315 32 3x3 1x1 8 4 104204 1776 60 34 1/315 32 3x3 1x1 8 4 4 4 104204 1776 61 34 1/315 32 5x5 rigid 8 96012 2292 156 28 1/315 32 5x5 3x3 8 4 104204 2804 171 31 1/315 32 5x5 3x3 8 4 4 4 104204 2804 172 31
SLIDE 52
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Numerical results: Decreasing h, k3h2 constant, constant number of elements per subdomain
Problem description Coarse directions Number of Subdomains Multipliers Wet interf. degrees of freedom
h k f e f e f e Orig Red Crs It 1/40 10 2x2 rigid 4 1632 68 11 10 1/40 10 2x2 1x1 4 1794 140 18 14 1/40 10 2x2 1x1 4 4 1794 140 34 11 1/80 16 4x4 rigid 4 6336 464 63 15 1/80 16 4x4 2x2 4 6914 674 84 25 1/80 16 4x4 2x2 4 4 6914 674 100 24 1/80 16 4x4 2x2 4 4 4 6914 674 102 24 1/80 16 4x4 2x2 4 4 4 4 6914 674 111 21 1/80 16 4x4 rigid 8 6336 464 102 11 1/80 16 4x4 2x2 8 4 6914 674 139 21 1/80 16 4x4 2x2 8 4 4 4 6914 674 150 18 1/80 16 4x4 2x2 8 4 8 8 6914 674 172 18 1/80 16 4x4 2x2 8 8 8 8 6914 674 183 13 1/160 25 8x8 rigid 8 24960 2240 462 16 1/160 25 8x8 4x4 8 4 27138 2930 543 33 1/160 25 8x8 4x4 8 4 4 4 27138 2930 601 31 1/160 25 8x8 4x4 8 4 8 8 27138 2930 635 30
SLIDE 53
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Numerical results: Decreasing h, k3h2 constant, same subdomains
Problem description Coarse directions Number of Subdomains Multipliers Wet interf. degrees of freedom
h k f e f e f e Orig Red Crs It 1/40 10 4x4 rigid 4 1632 240 63 11 1/40 10 4x4 2x2 4 1794 354 84 18 1/40 10 4x4 2x2 4 4 1794 354 100 16 1/80 16 4x4 rigid 4 6336 464 63 15 1/80 16 4x4 2x2 4 6914 674 84 25 1/80 16 4x4 2x2 4 4 4 4 6914 674 111 21 1/80 16 4x4 rigid 8 6336 464 102 11 1/80 16 4x4 2x2 8 4 6914 674 139 21 1/80 16 4x4 2x2 8 8 8 8 6914 674 183 13 1/160 25 4x4 rigid 8 24960 912 102 17 1/160 25 4x4 2x2 8 4 27138 1314 139 35 1/160 25 4x4 2x2 8 8 8 8 27138 1314 181 24 1/320 40 4x4 rigid 8 99072 1808 102 24 1/320 40 4x4 2x2 8 4 107522 2594 139 48 1/320 40 4x4 2x2 8 8 8 8 107522 2594 183 34
SLIDE 54
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Varying strength of artificial radiation conditions
10
−7
10
−6
10
−5
10
−4
10
−3
10
−2
10
−1
10 20 30 40 50 20 25 30 35 40 45 50 alpha0 k iter
h=1/40, 4x4 subdomains in fluid and 4x4 in solid, 8 coarse directions in fluid, 4 pressure waves and 4 shear waves in solid, scatterer size 0.4x0.4
SLIDE 55
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Varying strength of artificial radiation conditions
10
−7
10
−6
10
−5
10
−4
10
−3
10
−2
10
−1
10 20 30 40 50 20 25 30 35 40 45 50 alpha0 k iter
h=1/40, 4x4 subdomains in fluid and 4x4 in solid, 8 coarse directions in fluid, 4 pressure waves and 4 shear waves in solid, scatterer size 0.4x0.4, selection of basis instead of orthogonalization of the rows of B and the columns of Q
SLIDE 56
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Varying strength of artificial radiation conditions
10
−7
10
−6
10
−5
10
−4
10
−3
10
−2
10
−1
10 20 25 30 35 40 45 50 5 10 15 20 25 30 35 40 beta0 k iter
h=1/40, 4x4 subdomains in fluid and 1x1 in solid, 8 coarse directions in fluid, 4 pressure waves and 4 shear waves in solid, scatterer size 0.2x0.2
SLIDE 57
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Spectrum of preconditioned operator for FETI-H Fast convergence of FETI-H is due to clustering of the spectrum of the preconditioned operator at a point different from zero.
−5 5 10 15 20 25 30 −2 2 4 6 8 10 12 14 Re Im Complete spectrum
h=1/40, k=10, 4x4 subdomains, 4 coarse directions
SLIDE 58
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Spectrum of preconditioned operator for coupled problem The spectrum of the method for the coupled problem is very similar to the spectrum of FETI-H. Fast convergence is again due to clustering of the spectrum of the preconditioned operator at a point different from zero.
−5 5 10 15 20 25 30 −2 2 4 6 8 10 12 14 Re Im Complete spectrum
h=1/40, k=10, 4x4 and 2x2 fluid subdomains, 4 fluid coarse directions
SLIDE 59
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Computational test for unstructured 3D grids
- Numbers reported here from Matlab runs on imported data
- Coarse problem from eigenvectors of subdomain matrices
- Artificial radiation strength α0 = 0.1
- Cube scatterer 0.2x0.2x0.2 in cube waveguide 1x1x1, k=10
3362 fluid dofs and 1203 solid dofs on subdomain interfaces, reduced system (augmented by wet dofs) size 7114 unstructured grid, 16 fluid and 4 solid subdomains Number of iterations coarse/subdomain fluid only solid only decoupled coupled 12 47 54 93 126 20 35 44 90 90 40 28 32 88 88 Most of the slowdown in a coupled problem comes from solving the solid and fluid at the same time... not from the coupling
SLIDE 60
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Convergence for unstructured meshes 2D, coarse by directions
2D, 1x1 square discretized regularly by 400x400 Q1 elements. For the coupled case, there is a solid square 0.4x0.4 in the middle fluid solid coupled k subd dirs NC IT NC IT NC IN 25.12 4 + 4 61 88 5 44 26 100 104 34 84 + 16 219 456 1e-1 5 540 16 1260 30 736 33 62.8 4 + 4 111 6e-2 111 5 44 71 100 64 104 50 9 72 43 183 37 176 31 84 + 16 1.5e-6 703 8e-3 5 540 42 1260 41 736 59 9 900 13 1701 29 1148 22 125.6 4 + 4 203 2e-5 165 5 44 163 100 5e-5 104 105 9 72 131 196 47 188 67 17 128 67 296 25 294 43 33 228 12 84 + 16 967 8e-3 3e-1 5 540 483 1260 45 736 424 9 900 122 2160 25 1240 142 13 1260 12 1576 22
SLIDE 61
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Convergence for unstructured meshes 2D, coarse by directions
2D, 1x1 square discretized regularly by 400x400 Q1 elements. For the coupled case, there is a solid square 0.4x0.4 in the middle fluid solid coupled k subd dirs NC IT NC IT NC IN 125.6 4 + 4 203 2e-5 165 5 44 163 100 5e-5 104 105 9 72 131 196 47 188 67 17 128 67 296 25 294 43 33 228 12 84 + 16 967 8e-3 3e-1 5 540 483 1260 45 736 424 9 900 122 2160 25 1240 142 13 1260 12 1576 22 36 + 6 13 1182 52
SLIDE 62
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Convergence for unstructured meshes 3D, coarse by directions
3D, 1x1x1 cube discretized by P1 elements, for the coupled case, there is a solid cube 0.4x0.4x0.4 in the middle, 22396 nodes, 27548 dofs fluid solid coupled k subd dirs NC IT NC IT NC IN 12.5 4+8 0 366 349 7 318 54 1051 68 441 125 13 578 38 1644 58 767 91 21 870 38 2044 68 1088 85 15 + 32 383 787 644 1 255 143 1524 122 561 274 7 1678 37 5223 108 2454 153 13 2933 31 3827 101
SLIDE 63
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Convergence for unstructured meshes 3D, coarse by directions
3D, 1x1x1 cube discretized by P1 elements, for the coupled case, there is a solid cube 0.4x0.4x0.4 in the middle, 158514 nodes, 188073 dofs fluid solid coupled k subd dirs NC IT NC IT NC IN 25.6 16 + 32 1140 1787 7 1753 175 2758 390 13 3185 60 3243 182
SLIDE 64
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Spectrum of cube, fluid block only, no coarse
−15 −10 −5 5 10 15 20 −20 −15 −10 −5 5 10 3dcube reduced operator fluid block
SLIDE 65
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Spectrum of solid block operator, no coarse, 3D example
−500 −400 −300 −200 −100 100 −200 −150 −100 −50 50 100 150 200 250 3dcube reduced operator elastic block
SLIDE 66
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Spectrum for cube preconditioned by coarse
−20 −10 10 20 30 −30 −25 −20 −15 −10 −5 5 10 15 3dcube prec. red. blue=fluid, green=solid, red=all
- Spectrum of coupled problem is close to the union of the spectra of the
fluid and the solid blocks
- Coarse space tends to focus the spectrum along the positive real half-axis
SLIDE 67
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Spectrum of coupled operator preconditioned by coarse, 3D example
−20 −15 −10 −5 5 10 15 20 25 30 −25 −20 −15 −10 −5 5 10 15 3dcube prec. red. blue=fluid, green=solid, red=all
- For good convergence of GMRES or other Krylov space methods,
spectrum should be concentrated away from the origin
SLIDE 68
✫ ✬ ✪ ✩
University of Colorado Jan Mandel Conclusion
- FETI-H: good performance, numerically/computationally scalable
- coupled method reduces in the limit of rigid scatterer to FETI-H for solid
and fluid at the same time (triangular system)
- coupled method takes often about the same number of iterations as
FETI-H for rigid scatterer, sometimes the sum of numbers of iterations for fluid and solid separately
- coupled method prototype in Matlab, 3D parallel implementation Radek
Tezaur based on FETI code by Charbel Farhat and Michel Lesoinne
- MG: 2D Matlab prototype, very promising
SLIDE 69
✫ ✬ ✪ ✩