Adaptive BDDC algorithms for the system arising from plane wave discretization of Helmholtz equations
Jie Peng Junxian Wang Shi Shu
math.NA 65N30, 65F10, 65N55




XTU1]Jie Pengcof
xtu_pengjie@163.com

XTU1,XTU2]Junxian Wangcof
wangjunxian@xtu.edu.cn

XTU1,XTU2]Shi Shucor
shushi@xtu.edu.cn


[cof]These authors contributed equally to this work and should be considered co-first authors
[cor]Corresponding author
[XTU1]School of Mathematics and Computational Science, Xiangtan University, Xiangtan 411105, China
[XTU2]Hunan Key Laboratory for Computation and Simulation in Science and Engineering,
Xiangtan University, Xiangtan 411105, China






 Balancing domain decomposition by constraints (BDDC) algorithms with adaptive primal constraints are developed in a concise variational framework for the weighted plane wave least-squares (PWLS) discritization of Helmholtz equations with high and various wave numbers.
 The unknowns to be solved in this preconditioned system are defined on elements rather than vertices or edges,
 which are different from the well-known discritizations such as the classical finite element method.
 Through choosing suitable “interface" and appropriate primal constraints with complex coefficients and  introducing some local techniques, we developed a two-level adaptive BDDC algorithm for the PWLS discretization,
 and the condition number of the preconditioned system is proved to be bounded above by a user-defined tolerance and a constant which is only dependent on the maximum number of interfaces per subdomain. A multilevel algorithm is also attempted to resolve the bottleneck in large scale coarse problem.
 Numerical results are carried out to confirm the theoretical results and illustrate the efficiency of the proposed algorithms.












Helmholtz equation, high wave number, plane wave discretization, BDDC algorithm, adaptive primal constraints
[2010] 65N30 65F10 65N55


Adaptive BDDC algorithms for the system arising from plane wave discretization of Helmholtz equations
    [
    December 30, 2023
=====================================================================================================







§ INTRODUCTION


Helmholtz equations have many applications in  electromagnetic radiation, acoustics scattering and exploration seismology.







As the oscillatory behavior of the solution of the Helmholtz equation,
the corresponding discrete system is usually huge and highly indefinite, especially for high wave numbers.

The plane wave methods, which fall into the class of Trefftz methods <cit.>,

are popular discretization methods for solving this kind of equations <cit.>.
Compared with the classical finite element method (FEM) <cit.>, the plane wave methods can significantly
reduce the required degree of freedom under the same error precision, and with the increase of the wave number, the superiority is more obvious.
The weighted plane wave least-squares method (PWLS) is a frequently-used plane wave method <cit.>.

One advantage of PWLS over the other plane wave methods is that the stiffness matrix of the PWLS discrete system is Hermitian positive definite, this
lead to solve the resulting system by preconditioned conjugate gradient (PCG) method, and
the preconditioner plays an important role in the iterative process.



The development of an efficient solver or preconditioner for the Helmholtz equation has led to a great interest over the course
of the past decades <cit.>.

Domain decomposition (DD) methods are powerful parallel methods for solving the systems arising from finite element discretization of elliptic problems.
There exist many well known nonoverlapping DD methods for solving indefinite systems of Helmholtz equations, like the Robin-type DD method <cit.>,
the substructuring method <cit.>, the finite element tearing and interconnecting (FETI) method <cit.> and the dual-primal finite element tearing and interconnecting (FETI-DP) method <cit.>.

Alternative advanced nonoverlapping DD method is
the balancing domain decomposition by constraints (BDDC) methods <cit.>.

In the works by Li and Tu <cit.>, the BDDC, which incorporated some plane waves in the coarse problem to accelerate the convergence rate,  were extended to solving the FEM discrete system of Helmholtz equation. Numerical experiments illustrate that the convergence rate depends on a logarithmic pattern of the dimension of the local subdomain problems, improves with the decrease of the subdomain diameters, and depends on the wave number but it can be improved by including more plane wave continuity constraints in the coarse space.
Therefore, to enhance the robustness of the BDDC methods for solving the Helmholtz systems, the selection of good primal constraints should be necessary.

The main objective of this paper is to propose an adaptive BDDC preconditioner.
Adaptive BDDC preconditioner is an advanced BDDC method using a transformation of basis, the
primal unknowns are always selected by solving some generalized eigenvalue problems with respect to the local problems, and adaptively depends on a given tolerance <cit.>.
Since these local problems can indicate the bad behavior of the standard coarse problem,
they can be used to select the primal constraints to enhance the convergence of the iteration <cit.>.




However, there is an undeniable fact that the number of primal unknowns increases as the number of subdomains increases,
the corresponding coarse problem will become too large and hard to solve directly.


This leads to multilevel extension of this algorithm naturally <cit.>.

In this paper, we will develop an adaptive BDDC preconditioner. To be more specific,
we will extend the existing methods in <cit.> to  PWLS discretization of Helmholtz equation with high and various wave numbers,
and present a complete theory.

Contrast to the classical FEM, the dofs in the PWLS method are defined on elements rather than vertices or edges,
we thus introduce a kind of special “interface",
which is different from the existing nonoverlapping DD methods.





Since the PWLS discrete system consists of complex coefficients, we construct our transformation operators by using a series of local generalized eigenvalue problems with respect to the parallel sum and the primal constraints are formed by the eigenvectors with their complex modulus of eigenvalues greater than a given tolerance Θ, which is different from <cit.>.
As the spectral condition number of the PWLS discretizations of the Helmholtz equations with high wave numbers
grows with the increase of the number of plane wave bases in each element and the decrease of the grid size <cit.>,
some local techniques are introduced to overcome this difficulty.
Then, by introducing some other auxiliary spaces and operators, we arrive at our two-level adaptive BDDC algorithm in variational framework for PWLS discretizations.
The condition number bound of the two-level adaptive BDDC preconditioned systems, CΘ, can be derived by using the properties of the auxiliary spaces and involved operators, where C is a constant which depends only on the maximum number of interfaces per subdomain. Compared with the previous work for mortar discretizations in <cit.>, the variational framework in this paper is more concise.


We perform numerical experiments for various model problems.
These results verify the correctness of theoretical results, and show that our two-level adaptive BDDC algorithms are scalability with respect to the angular frequency, the number of subdomains and mesh size.
It is worth pointing out that the algorithm with deluxe scaling matrices has more advantages over the algorithm with multiplicity scaling matrices in the size of coarse problem even for model problem with constant medium, which is different from the adaptive BDDC algorithms for the two-order elliptic problems.



However, since the number of primal unknowns increase as the wave number or the number of subdomains increase,
we attempt to construct a multilevel adaptive BDDC algorithm to resolve the bottleneck in solving large scale coarse problem.

The coarser subdomains are gathered by a certain amount of subdomains at the finer level, and the new “interface" can be obtained naturally.
Numerical results show that the multilevel algorithm can reduce the size of the coarse problem, and it is also robust to solve the Helmholtz equations with high wave number.


The rest of this paper is organized as follows. In Section 2, The PWLS formulation will be presented for a Helmholtz equation with Robin boundary condition.
In Section 3, we will firstly derive the Schur complete variational problem and its corresponding function space, introduce some auxiliary spaces and dual-primal basis functions further, and then illustrate our two-level and multilevel adaptive BDDC algorithms. An estimate of the condition numbers will be analyzed in Section 4, and various numerical experiments are presented to verify the performance of our algorithm in Section 5.
Finally, we will give a conclusion in Section 6.







§ WEIGHTED PLANE WAVE LEAST SQUARES FORMULATION










 §.§ Model problem


Let Ω∈ℝ^2 be a bounded and connected Lipschitz domain with boundary ∂Ω.
Consider the Helmholtz equation with Robin boundary condition (<cit.>)

    {[  -Δ u  - κ^2 u              =              0          in Ω,; (∂_ n + i κ )u              =              g         on ∂Ω, ].

where i = √(-1), ∂_ n and κ are separately the imaginary unit, the outer normal derivative and the wave number, g ∈ L^2(Ω).
The wave number κ = ω/c > 0, where ω and c are separately called the angular frequency and the wave speed.
The above problem is usually seen as an approximation of the acoustic scattering problem, and the wave speed c (and hence κ)  can be a constant or variable function. For g = 0, the second equation of (<ref>) becomes a general representation of an absorbing boundary condition <cit.>.




 §.§ Weighted plane wave least squares discretization




Following Hu and Zhang <cit.>, we first define a quadrilateral mesh 𝒯_h, namely, dividing Ω into

    Ω̅ = ⋃_k=1^N_hΩ̅_k,

where the quadrilateral elements {Ω_k} satisfy that Ω_m ∩Ω_l = ∅, m ≠ l, h_k is the size of Ω_k and h = max_1≤ k ≤ N_h h_k.
Define

    [ γ_kj = ∂Ω_k ∩∂Ω_j,   k,j = 1,⋯,N_h    k≠ j,;                γ_k = ∂Ω_k ∩∂Ω,   k=1,⋯,N_h,;     ℱ_B = ⋃_k=1^N_hγ_k,  ℱ_I = ⋃_k ≠ jγ_kj. ]


Throughout this paper, we assume that each
κ_k:= κ|_Ω_k is a constant. Let V(Ω_k) be the local space whose members satisfy the
homogeneous Helmholtz's equation (<ref>) on Ω_k:

    V(Ω_k) = {v_k ∈ H^1(Ω_k): Δ v_k + κ_k^2 v_k = 0}, k=1,⋯,N_h.

Define the global space

    V(𝒯_h) = ⋃_k=1^N_h V(Ω_k).






As we all know,
problem (<ref>) to be solved is equivalent to find the local solution
u_k: = u|_Ω_k∈{v ∈ H^1(Ω_k): ∇ v ∈ H(div;Ω_k)}
such that

    {[ -Δ u_k  -κ_k^2 u_k                  =                  0           in  Ω_k,                   ;  (∂_ n + iκ_k )u_k                  =                  g            on γ_k,                    ]. k=1,2,⋯,N_h,

with the continuity conditions for u and its normal derivative on the interfaces between the
elements:

    u_k  - u_j  =  0, ∂_ n_ku_k +  ∂_ n_ju_j =   0,   γ_kj,  k,j = 1,⋯,N_h and k≠ j.




In the weighted plane wave least squares (PWLS) formulation, a finite dimensional subspace of  V(𝒯_h) is introduced,

    V_p(𝒯_h) =span{φ_m,l: 1 ≤ l ≤ p, 1≤ m ≤ N_h},

where

    φ_m,l()=
    {[ y_m,l()   ∈Ω̅_m;       0 ∈Ω\Ω̅_m ]. ,

here y_m,l (l=1,⋯,p) denote the wave shape functions on Ω_m, which satisfy

    {[           y_m,l()                 = e^iκ(·_l), ∈Ω̅_m,;              |_l|                 =                1,;                _l                 ≠       _s,  l ≠ s, ].

and _l (l = 1,⋯,p) are unit wave propagation directions.
In particular, during numerical simulations, we set

    _l = (
      [ cos(2π(l-1)/p); sin(2π(l-1)/p);                ]).


By using the plane wave finite dimensional space V_p(𝒯_h) defined above, the PWLS formulation associated with problem (<ref>) and (<ref>) can be described as follows: find u ∈ V_p(𝒯_h) such that

    a(u,v) =  ℒ(v), ∀ v ∈ V_p(𝒯_h),

where

    a(u,v)    =   ∑_j≠ k (α_kj∫_γ_kj (u_k-u_j)·(v_k-v_j)ds
    +β_kj∫_γ_kj
    (∂_ n_ku_k+∂_ n_ju_j)·(∂_ n_kv_k+∂_ n_jv_j)ds)
      
           + ∑_k=1^N_hθ_k ∫_γ_k ((∂_ n+iκ_k )u_k) ·(∂_ n+i κ_k)v_kds,
    ℒ(v)    =   ∑_k=1^N_hθ_k ∫_γ_k g·(∂_ n+i κ_k )v_kds,

here ♢ denotes the complex conjugate of the complex quantity ♢, the Lagrange multipliers

    α_kj = h^-1 + κ_kj, β_kj = h^-1κ_kj^-2 + κ_kj^-1  κ_kj = (κ_k + κ_j)/2,

and

    θ_k = h^-1κ_k^-2 + κ_k^-1.


This discrete variational problem (<ref>) is derived by the minimization of a quadratic functional,

 and the basic idea of the minimization problem is to find a function in V_p(𝒯_h) so that it can satisfy the external boundary conditions and the interface conditions as far as possible <cit.>. From Theorem 3.1 of <cit.>, we can see that the continuous variational problem associate with (<ref>) is equivalent to the reference problem (<ref>) and (<ref>).


It is clear that a(·,·) is sesquilinear and Hermitian, and similar to  the proof of Theorem 3.1 in <cit.>, we can see that a(v,v) ≥ 0, and a(v,v) = 0 for v ∈ V(𝒯_h) if and only if v = 0. Therefore,  a(·,·) is Hermitian positive definite (HPD).


Due to the linear system obtained in (<ref>) are large and highly ill-conditioned when the
wave number is large, it's necessary to study a fast solver for this system.
Adaptive BDDC algorithm is a novel domain decomposition (DD) method with enriched coarse spaces <cit.>,
this algorithm has been successfully applied to solve discrete systems obtained by various discretization methods, such as conforming Galerkin <cit.>,
discontinuous Galerkin <cit.>, and mortar methods <cit.> and so on.  However, adaptive BDDC
algorithm for PWLS discretizations has not previously been discussed in the literature.
Here we will extend the adaptive BDDC algorithm to PWLS discretizations with high and various wave numbers.

















§ ADAPTIVE BDDC PRECONDITIONER




 §.§ Domain decomposition and Schur complement problem


Differ from the discretizations which dofs are defined on the vertices or edges of the mesh,
the dofs in the PWLS discretization are defined on the elements,
therefore, we need to introduce a special interface and domain decomposition firstly.

Let {D_r}_r=1^N_d be a non-overlapping subdomain partiaon of Ω and each D_r consists of several complete elements and part of the elements in 𝒯_h (see
Figure <ref>).




Let 𝒯_d denote the coarse partition associated with the subdomains D_1, D_2, ⋯, D_N_d.



For any given r = 1,⋯, N_d, define

    ℱ_I^(r) = {γ̃_kj: γ̃_kj = γ_kj|_D̅_r, ∀γ_kj∈ℱ_I}, ℱ_B^(r) = {γ̃_k: γ̃_k = γ_k|_D̅_r, ∀γ_k∈ℱ_B}.

From this definition, the sesquilinear form a(·,·) defined in (<ref>) can be rewrite as

    a(u,v) = ∑_r=1^N_d a_r (u,v),

where

    a_r(u,v)    =   ∑_γ̃_kj∈ℱ_I^(r) (α_kj∫_γ̃_kj (u_k-u_j) ·(v_k-v_j)ds + β_kj∫_γ̃_kj
    (∂_ n_ku_k+∂_ n_ju_j) ·(∂_ n_kv_k+∂_ n_jv_j)ds)
      
          + ∑_γ̃_k ∈ℱ_B^(r)θ_k ∫_γ̃_k ((∂_ n+iw )u_k) ·(∂_ n+i ω)v_kds.

and it is easy to verify that a_r(·,·) is Hermitian positive semi-definite.


If ∂ D_r ∩∂ D_j(r ≠ j) is a common edge, we call it an interface Γ_rj (or Γ_jr).
Let Γ=∪Γ_rj, N_f denotes the number of interfaces.
Since there is one-to-one correspondence between
any given Γ_rj and the interface, we denote {Γ_rj} by {F_k: k=1,⋯,N_f} for convenience.
Denote {v_k: k=1,⋯,N_v} be the set of the vertices corresponding to 𝒯_d.
For each D_r (r = 1,⋯,N_d), let






    ℳ_r := {k: F_k ⊂∂ D_r \∂Ω,  1≤ k ≤ N_f},
    
          ℳ_r^c := {k: v_k ⊂∂ D_r \∂Ω,  1≤ k ≤ N_v}.

For sake of argument, we set the elements of ℳ_r as {r_m: 1≤ m ≤ f_r}, where f_r denotes the size of ℳ_r, i.e.
the number of interfaces on ∂Ω_r \∂Ω.


Set

    𝒟_r    := {m: Ω_m ∩ D_r ≠∅,  1 ≤ m ≤ N_h}, 1 ≤ r ≤ N_d,
    
    ℐ_r    := {m: Ω_m ⊂ D_r,  1 ≤ m ≤ N_h}, 1 ≤ r ≤ N_d,
    
    𝒱_k    := {m: Ω_m ∩ v_k ≠∅,  1 ≤ m ≤ N_h}, 1 ≤ k ≤ N_v,
    
    ℱ_k    := {m: Ω_m ∩ F_k ≠∅,  1 ≤ m ≤ N_h  m ∉𝒱_l, 1 ≤ l ≤ N_v}, 1 ≤ k ≤ N_f.

and (see Figure <ref> for illustration)

    D_v_k = ⋃_m ∈𝒱_kΩ_m  D_F_k = ⋃_m ∈ℱ_kΩ_m.

Then the special interface subdomain can be defined as

    D_Γ = (⋃_k=1^N_v D_v_k) ⋃(⋃_k=1^N_f D_F_k).











Based on the aforementioned domain decomposition, the Schur complement problem of (<ref>) can be introduced.


Let the local spaces

    V_I = ⊕_r=1^N_dV^(r)_I,  V^(r)_I = span{φ_m,l: 1 ≤ l ≤ p, ∀ m ∈ℐ_r}, 1 ≤ r ≤ N_d,
    
           V_k,c = span{φ_m,l: 1 ≤ l ≤ p, ∀ m∈𝒱_k}, 1 ≤ k ≤ N_v,

and

    V_k = span{φ_m,l: 1 ≤ l ≤ p, ∀ m ∈ℱ_k}, 1 ≤ k ≤ N_f.



For any interface F_k, we always assume that it is a common edge of two subdomains D_r and D_j.
Let n_k = p |ℱ_k|, where |♢| denotes the size of set ♢. Define a vector

    Φ^k = (ϕ^k_1,⋯,ϕ^k_n_k)^T:=(ϕ^k_m,1,⋯,ϕ^k_m,p, ∀ m ∈ℱ_k)^T,

where ϕ^k_m,l∈ V_I^(r)⊕ V_I^(j)⊕ V_k (1 ≤ l ≤ p, ∀ m∈ℱ_k) satisfy

    {[  a(ϕ^k_m,l, v) = 0,   ∀ v ∈ V_I,; ϕ^k_m,l|_D̅_F_k = φ_m,l|_D̅_F_k. ].


For any vertex v_k, let

    𝒩_k = {r: ∂ D_r ∩ v_k ≠∅, 1≤ r≤ N_d},  n_k^c = p|𝒱_k|,

we can define another vector

    Ψ^k = (ψ^k_1,⋯,ψ^k_n_k^c)^T:=(ψ^k_m,1,⋯,ψ^k_m,p, ∀ m∈𝒱_k)^T,

where ψ^k_m,l∈ (⊕_r ∈𝒩_k V_I^(r)) ⊕ V_k,c (1 ≤ l ≤ p, ∀ m∈𝒱_k) satisfy

    {[  a(ψ^k_m,l, v) = 0,   ∀ v ∈ V_I,; ψ^k_m,l|_D̅_v_k = φ_m,l|_D̅_v_k. ].












Utilizing these two vectors, the function space related to the Schur complement problem can be defined as

    Ŵ = (⊕_k=1^N_f W_k) ⊕ (⊕_k=1^N_v W_k,c),

where

    W_k = span{ϕ^k_1,⋯,ϕ^k_n_k}  W_k,c = span{ψ^k_1,⋯,ψ^k_n_k^c}.


Let Ŝ: Ŵ→Ŵ be the Schur complement operator defined by

    (Ŝû, v̂) = a(û, v̂),  ∀û, v̂∈Ŵ.

Due to the coercive of the restriction of a(·,·) on Ŵ, Ŝ is a HPD operator.

Then the Schur complement variational form of (<ref>) can be expressed as:
find ŵ∈Ŵ such that

    (Ŝŵ, v̂) = ℒ(v̂),  ∀v̂∈Ŵ.




In order to give an adaptive BDDC preconditioner for solving the Schur complement problem (<ref>),
we will give some preparations in next subsection.



 §.§ Some auxiliary spaces and dual-primal bases


The BDDC algorithm consist of several independent subdomain problems and one global coarse problem.
The function spaces associated with these problems are spanned by the so-called dual basis functions and primal basis functions.
In the following, we will introduce some auxiliary spaces and dual-primal bases.

For each subdomain D_r, let truncated basis functions

    φ_m,l^(r)()=
    {[          φ_m,l()      ∈Ω̅_m ∩D̅_r;                0 ∈Ω\ (Ω̅_m ∩D̅_r) ]. ,  1≤ l ≤ p, m ∈𝒟_r.

Define the corresponding local spaces

    V_k,c^(r) = span{φ_m,l^(r): 1 ≤ l ≤ p, ∀ m ∈𝒱_k},  k ∈ℳ_r^c,
    
        V_k^(r) = span{φ_m,l^(r): 1 ≤ l ≤ p, ∀ m ∈ℱ_k},  k ∈ℳ_r,

and

    V^(r) = V^(r)_I ⊕ (⊕_k∈ℳ_r V_k^(r)) ⊕ (⊕_k∈ℳ_r^c V_k,c^(r)),

where V^(r)_I is defined in (<ref>).

Using the above-mentioned local spaces, some new set of basis functions and corresponding spaces can be constructed firstly.

For any given interface F_k = ∂ D_r ∩∂ D_j, denote

    Φ^k,ν =  (ϕ^k,ν_m,1, ⋯, ϕ^k,ν_m,p, ∀ m∈ℱ_k)^T  Φ̅^k,ν = (ϕ̅^k,ν_m,1, ⋯, ϕ̅^k,ν_m,p, ∀ m∈ℱ_k)^T, ν=r,j,

where ϕ^k,ν_m,l∈ V_I^(ν)⊕ V_k^(ν) and ϕ̅^k,ν_m,l∈ V^(ν) separately satisfy

    {[ a_ν(ϕ^k,ν_m,l, v) = 0,   ∀ v ∈ V^(ν)_I;  ϕ^k,ν_m,l|_D̅_F_k = φ^(ν)_m,l|_D̅_F_k ].,  1 ≤ l ≤ p, m ∈ℱ_k,

and

    {[ a_ν(ϕ̅^k,ν_m,l, v) = 0,   ∀ v ∈ V^(ν)\ V_k^(ν);         ϕ̅^k,ν_m,l|_D̅_F_k = φ^(ν)_m,l|_D̅_F_k ].,  1 ≤ l ≤ p, m ∈ℱ_k.



For any given vertex v_k, denote

    Ψ^k,ν = (ψ^k,ν_m,1, ⋯, ψ^k,ν_m,p, ∀ m∈𝒱_k)^T, ν∈𝒩_k,

where ψ^k,ν_m,l∈ V_I^(ν)⊕ V_k,c^(ν) satisfy

    {[ a_ν(ψ^k,ν_m,l, v) = 0,   ∀ v ∈ V^(ν)_I;  ψ^k,ν_m,l|_D̅_v_k = φ^(ν)_m,l|_D̅_v_k ].,  1 ≤ l ≤ p, m ∈𝒱_k.

and 𝒩_k is defined in (<ref>).



From the following lemma, we know that the basis functions ϕ^k,ν_m,l and ϕ̅^k,ν_m,l (1≤ l ≤ p, ∀ m ∈ℱ_k) are available.

The restriction of a_r(·,·) (r=1,⋯,N_d) on V^(ν)\ V_k^(ν) or V_I^(ν) is HPD.



First, from the definition (<ref>) of a_r(·,·), it's easy to verify that a_r(u,u) ≥ 0 for any u ∈ U (U = V^(ν)\ V_k^(ν)  V_I^(ν)).

Next, we show a_r(u,u) = 0 holds if and only if u=0. For simplicity, we only give the proof of the case U = V^(ν)\ V_k^(ν),
it can be proved in a similar way if U = V^(ν)_I. 

In fact, a_r(u,u) = 0 is equivalent to

    ∑_γ̃_mj∈ℱ_I^(r) (α_mj∫_γ̃_mj |u_m-u_j|^2 ds
    + β_mj∫_γ̃_mj |∂_ n_mu_m+∂_ n_ju_j|^2 ds) 
        + ∑_γ̃_m ∈ℱ_B^(r)θ_m ∫_γ̃_m |(∂_ n+iw )u_m|^2 ds = 0.

Note that α_mj, β_mj, θ_m >0, the above equality implies that

    ∫_γ̃_mj |u_m-u_j|^2 ds = 0, ∫_γ̃_mj |∂_ n_mu_m+∂_ n_ju_j|^2 ds = 0, ∫_γ̃_m |(∂_ n+iw )u_m|^2 ds = 0.


From the definition of space V^(ν)\ V_k^(ν), we have

    u|_D_F_k = 0, ∀ u ∈ V^(ν)\ V_k^(ν).

It implies that if γ̃_mj⊂∂ D_F_k and  m ∈ℱ_k, we have

    ∫_γ̃_mj |u_m - u_j|^2 ds = 0 ⇔ ∫_γ̃_mj |u_j|^2 ds = 0.

Namely,

    u_j = ∑_l=1^p u_j,lφ_j,l = 0,  γ̃_mj,

where u_j,l∈ℂ (l=1,⋯,p).

From this and based on the linear independence of basis functions φ_j,l(l=1,⋯,p), we obtain

    u_j,l = 0, l=1,⋯,p ⇔ u_j = 0,  Ω_j,

where j satisfies γ̃_mj⊂∂ D_F_k\∂ D_r and m ∈ℱ_k.

Furthermore, applying this process recursively, we can prove that u = 0  D_r.

Therefore we complete the proof of the coercive of a_r(·,·) in V^(ν)\ V_k^(ν).


From (<ref>), (<ref>), (<ref>) and (<ref>), we can easily verify that

    ϕ^k,ν_m,l|_D̅_ν = ϕ^k_m,l|_D̅_ν  ψ^k,ν_m,l|_D̅_ν = ψ^k_m,l|_D̅_ν.





Similar to Φ^k and Ψ^k expressed in (<ref>) and (<ref>), we denote

    Φ^k,ν = (ϕ^k,ν_1, ⋯, ϕ^k,ν_n_k)^T, Φ̅^k,ν = (ϕ̅^k,ν_1, ⋯, ϕ̅^k,ν_n_k)^T, ν=r,j,

and

    Ψ^k,ν = (ψ^k,ν_1, ⋯, ψ^k,ν_n_k^c)^T, ν∈𝒩_k.



Using the basis functions defined above, denote the auxiliary spaces

    W_k^(ν)=span{ϕ^k,ν_1, ⋯, ϕ^k,ν_n_k},  W̅_k^(ν)=span{ϕ̅^k,ν_1, ⋯, ϕ̅^k,ν_n_k},  ν=r,j,

and

    W_k,c^(ν)=span{ψ^k,ν_1, ⋯, ψ^k,ν_n_k^c},  ∀ν∈𝒩_k.


For any given subdomain D_s, let

    |·|^2_a_s := a_s(·, ·).

The following lemma can be proved.

For any given subdomain D_s and vectors {w⃗_m ∈ℂ^n_m, m ∈ℳ_s}, there have

    ∑_m ∈ℳ_s |w̅_m^(s)|_a_s^2 ≤ f_s |w^(s)|_a_s^2, s = 1,⋯,N_d,

where

    w^(s) = ∑_m ∈ℳ_s w_m^(s) + ∑_m ∈ℳ_s^c w_m,c^(s),
     w_m^(s) = (w⃗_m^(s))^T Φ^m,s,  w_m,c^(s)∈ W_m,c^(s), w̅_m^(s) = (w⃗_m^(s))^T Φ̅^m,s,

here ℳ_s, ℳ_s^c are separately defined in (<ref>) and (<ref>), f_s = |ℳ_s|.


For each m ∈ℳ_s, by using (<ref>), we obtain

    w^(s)
    = w_m^(s) + ∑_α∈ℳ_s α≠ m w_α^(s) + ∑_α∈ℳ_s^c w_α,c^(s)
    = w̅_m^(s) + η^(s),

where

    η^(s) = ∑_α∈ℳ_s α≠ m w_α^(s) + ∑_α∈ℳ_s^c w_α,c^(s) + (w_m^(s) - w̅_m^(s)) ∈ V^(s)\
    V_m^(s).


Therefore, using (<ref>), (<ref>), (<ref>), (<ref>) and Lemma <ref>, we have

    a_s(w^(s), w^(s))
    = a_s(w̅_m^(s) + η^(s),w̅_m^(s) + η^(s))
    =a_s(w̅_m^(s),w̅_m^(s)) + a_s(η^(s), η^(s))
    ≥ a_s(w̅_m^(s),w̅_m^(s)).

From this and the definition (<ref>) of |·|^2_a_s, we can complete the proof of (<ref>).







Then, for any given ν=r,j, we introduce the scaling operator D_F_k^(ν): U → U (U = W_k, W_k^(r)  W_k^(j)) or scaling matrix D⃗_F_k^(ν)∈ℂ^n_k× n_k, which satisfy that for all w = w⃗^T Ψ with w⃗∈ℂ^n_k and Ψ = Φ^k, Φ^k,r  Φ^k,j, we have

    D_F_k^(ν) w = w⃗^T (D⃗_F_k^(ν))^T Ψ,

where D⃗_F_k^(ν) is nonsingular, and

    D_F_k^(r) + D_F_k^(j)  = I,

here I is the identity operator.



Two commonly used scaling matrices D⃗_F_k^(ν)(ν=r,j) 
are the multiplicity scaling matrices

    D⃗^(r)_F_k=1/2I⃗, D⃗^(j)_F_k=1/2I⃗,

and the deluxe scaling matrices (<cit.>)

    D⃗^(r)_F_k=(S⃗^(r)_F_k +
    S⃗^(j)_F_k)^-1S⃗^(r)_F_k, D⃗^(j)_F_k=(S⃗^(r)_F_k +
    S⃗^(j)_F_k)^-1S⃗^(j)_F_k,

where I⃗ denotes the n_k × n_k identity matrix, and

    S⃗^(ν)_F_k  = (a_l,m^(ν))_n_k× n_k, a_l,m^(ν) = a_ν(ϕ^k,ν_m, ϕ^k,ν_l), l,m=1,⋯, n_k, ν=r,j.



Since the spectral condition number of the plane wave discretizations of the Helmholtz equations with high wave numbers
grows like h^-q, where q is proportional to the number of plane wave bases in each element p (see <cit.>),
some local techniques are introduced, for example,
by using the incomplete LU factorization preconditioner to get the deluxe scaling matrices.


Using the above-mentioned scaling operators or matrices, a new set of basis functions of W_k^(ν) (ν = r,j) can be defined as





 
    Φ^k,r_D=D_F_k^(j)Φ^k,r=(D⃗_F_k^(j))^T
    Φ^k,r,  Φ^k,j_D=D_F_k^(r)Φ^k,j=(D⃗_F_k^(r))^T Φ^k,j.





In order to defined the so-called dual-primal basis functions,
we need to use the function spaces W_k^(ν) and W̅_k^(ν) (ν = r,j) defined in (<ref>),
and the scaling operators defined in (<ref>), to introduce a class of transformation operators (or matrices).


Set n_k = n_Δ^k + n_Π^k, where the integer n_Δ^k, n_Π^k ≥ 0. Let
n_k-order complex nonsingular matrix

    T⃗_F_k=(T⃗_Δ^F_k, T⃗_Π^F_k),

where T⃗_Δ^F_k and T⃗_Π^F_k are separately n_k× n_Δ^k and n_k× n_Π^k matrices.

For any ν=r,j, using the matrix T⃗_F_k, introduce the linear operators T_F_k.
These operators transform the basis vectors Φ̅^k,ν and Φ^k,ν_D into

    Φ̅^k, ν_T    = T_F_kΦ̅^k,ν:=(
    [ Φ̅^k, ν_Δ;  Φ̅^k,ν_Π ]), Φ^k, ν_T_D=T_F_kΦ^k,ν_D:=(
    [ Φ^k, ν_D,Δ;  Φ^k,ν_D,Π ]),

where

    Φ̅^k,ν_χ = (ϕ̅^k,ν_χ,1,⋯,ϕ̅^k,ν_χ,n_χ^k)^T  = T_χ^F_kΦ̅^k,ν=
     (T⃗_χ^F_k)^T  Φ̅^k,ν,
     Φ^k,ν_D,χ=T_χ^F_kΦ^k,ν_D= (T⃗_χ^F_k)^T  Φ^k,ν_D, χ=Δ,
    Π

are the so-called dual-primal basis functions.

For any given real number Θ≥ 1, the above operator T_F_k must satisfies
 
    |w_D,Δ^k,r|^2_a_r +
    |w̃_D,Δ^k,j|^2_a_j≤Θ
    |w̅_k,Δ^(r) +w̅_k,Π^(r)|^2_a_r,

where

    w_D,Δ^k,r =(w⃗_Δ)^T
    Φ^k,r_D,Δ, w̃_D,Δ^k,j =(w⃗_Δ)^T
    Φ^k,j_D,Δ, w̅_k,χ^(r)=
    (w⃗_χ)^T  Φ̅^k,r_χ, χ=Δ, Π,

here w⃗_Δ∈ℂ^n_Δ^k, w⃗_Π∈ℂ^n_Π^k are any given vectors.

The inequality (<ref>) plays a crucial role in the estimation of the condition number of the adaptive BDDC algorithm, and it is always
be replaced by lazy eigenanalysis.  
Following <cit.>, we introduce the matrices

    S⃗̅⃗^(ν)_F_k  = (b_l,m^(ν))_n_k× n_k, b_l,m^(ν) = a_ν(ϕ̅^k,ν_m, ϕ̅^k,ν_l), l,m=1,⋯, n_k, ν=r,j,

where the sesquilinear form a_ν(·,·) and the basis functions {ϕ̅^k,ν_l}_l=1^n_k are seperately defined in (<ref>) and (<ref>).


Then considering a generalized eigenvalue problem (see <cit.>)

    A⃗_F_k^D v⃗
     = λB⃗_F_kv⃗,

where


    A⃗_F_k^D = (D⃗_F_k^(r))^H S⃗_F_k^(j)D⃗_F_k^(r) + (D⃗_F_k^(j))^H S⃗_F_k^(r)D⃗_F_k^(j),
      B⃗_F_k = S⃗̅⃗_F_k^(r):S⃗̅⃗_F_k^(j),

here ♢^H denotes the conjugate transpose of ♢,  v⃗∈ℂ^n_k, λ∈ℂ,
D⃗_F_k^(ν)(ν=r,j), S⃗^(ν)_F_k(ν=r,j) are separately defined in (<ref>) and (<ref>),
and the parallel sum


    B⃗_F_k =
    S⃗̅⃗_F_k^(j)(S⃗̅⃗_F_k^(r)+S⃗̅⃗_F_k^(j))^†S⃗̅⃗_F_k^(r),

here
(S⃗̅⃗_F_k^(r)+S⃗̅⃗_F_k^(j))^†
is a pseudo inverse of the matrix
S⃗̅⃗_F_k^(r)+S⃗̅⃗_F_k^(j).

Since S⃗̅⃗_F_k^(ν)(ν=r,j) are both Hermitian positive semi-definite,  S⃗̅⃗_F_k^(r):S⃗̅⃗_F_k^(j)
is also Hermitian positive semi-definite and satisfies the following spectral inequalities <cit.>

    B⃗_F_k≤S⃗̅⃗_F_k^(ν), ν=r,j.



Let

    |λ_1| ≤ |λ_2| ≤⋯≤ |λ_n_Δ^k| ≤Θ≤ |λ_n_Δ^k+1| ≤⋯≤ |λ_n_k|,

where λ_k (k=1,⋯,n_k) is the eigenvalue of (<ref>), n_Δ^k is a non-negative integer, and Θ≥ 1 is given in (<ref>).

Denote T⃗_Δ^F_k and T⃗_Π^F_k in the n_k × n_k transform matrix T⃗_F_k defined in (<ref>) as

    T⃗_Δ^F_k:= (v⃗_1,⋯,v⃗_n_Δ^k),  T⃗_Π^F_k:= (v⃗_n_Δ^k+1,⋯,v⃗_n_k),

here v⃗_l(l=1,⋯,n_k) are the generalized eigenvectors  of (<ref>) corresponding to λ_l and
satisfy

    v⃗_l^H A⃗_F_k^D v⃗_m = v⃗_l^H B⃗_F_kv⃗_m = 0,  l ≠ m.













From this, we have

    (T⃗_Δ^F_k)^H C⃗T⃗_Δ^F_k= diag
    (v⃗_1^H C⃗v⃗_1,⋯, v⃗_n_Δ^k^H C⃗v⃗_n_Δ^k),  C⃗ = A⃗_F_k^D, B⃗_F_k,

and

    (T⃗_Π^F_k)^H B⃗_F_kT⃗_Δ^F_k=0,  (T⃗_Δ^F_k)^H B⃗_F_kT⃗_Π^F_k = 0.

Using (<ref>), we can prove that

    v⃗_l^H A⃗_F_k^D v⃗_l = |λ_l| v⃗_l^H B⃗_F_kv⃗_l, l=1,⋯,n_Δ^k.

































 Using the above matrix T⃗_F_k, we can obtain the operator
T_F_k defined in (<ref>). Next, we verify that it satisfies (<ref>).

By (<ref>) and (<ref>), we can rewrite the functions in (<ref>) as

    w_D,Δ^k,r =(w⃗_Δ)^T
    (T⃗_Δ^F_k)^T  (D⃗_F_k^(j))^T Φ^k,r, w̃_D,Δ^k,j =(w⃗_Δ)^T
    (T⃗_Δ^F_k)^T  (D⃗_F_k^(r))^T Φ^k,j,

and

    w̅_χ^k,r=
    (w⃗_χ)^T (T⃗_χ^F_k)^T  Φ̅^k,ν, χ=Δ, Π.


From (<ref>), (<ref>), the property of the sesquilinear form a_r(·,·) and (<ref>), it is easy to verify that

    |w_D,Δ^k,r|^2_a_r = a_r(w_D,Δ^k,r, w_D,Δ^k,r)
    
    
    = w⃗_Δ^H(T⃗_Δ^F_k)^H (D⃗_F_k^(j))^HS⃗_F_k^(r)D⃗_F_k^(j)T⃗_Δ^F_kw⃗_Δ.

Similarly, we have
 
    |w̃_D,Δ^k,j|^2_a_j = w⃗_Δ^H (T⃗_Δ^F_k)^H (D⃗_F_k^(r))^H S⃗_F_k^(j)D⃗_F_k^(r)T⃗_Δ^F_kw⃗_Δ.


By using (<ref>), (<ref>), (<ref>) and (<ref>),
we can obtain

    |w_D,Δ^k,r|^2_a_r +
    |w̃_D,Δ^k,j|^2_a_j = w⃗_Δ^H
     (T⃗_Δ^F_k)^H A⃗_F_k^D  T⃗_Δ^F_kw⃗_Δ
    = w⃗_Δ^H diag{v⃗_1^H A⃗_F_k^D v⃗_1,⋯, v⃗_n_Δ^k^H A⃗_F_k^D v⃗_n_Δ^k}w⃗_Δ.


From this, and using (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and note that B⃗_F_k is Hermitian positive semi-definite, we known that

    |w_D,Δ^k,r|^2_a_r +
    |w̃_D,Δ^k,j|^2_a_j
       = w⃗_Δ^H diag{|λ_1| v⃗_1^H B⃗_F_kv⃗_1,⋯, |λ_n_Δ^k| v⃗_n_Δ^k^H B⃗_F_kv⃗_n_Δ^k}w⃗_Δ
    
    
    
    
    
       ≤Θw⃗_Δ^H diag{v⃗_1^H B⃗_F_kv⃗_1,⋯, v⃗_n_Δ^k^H B⃗_F_kv⃗_n_Δ^k}w⃗_Δ
       = Θw⃗_Δ^H (T⃗_Δ^F_k)^H B⃗_F_kT⃗_Δ^F_kw⃗_Δ
    
    
       ≤Θ (T⃗_Δ^F_kw⃗_Δ + T⃗_Π^F_kw⃗_Π)^H B⃗_F_k (T⃗_Δ^F_kw⃗_Δ + T⃗_Π^F_kw⃗_Π)
    
    
       ≤Θ (T⃗_Δ^F_kw⃗_Δ + T⃗_Π^F_kw⃗_Π)^H S⃗̅⃗_F_k^(r) (T⃗_Δ^F_kw⃗_Δ + T⃗_Π^F_kw⃗_Π)
    
    
    
       = Θ
    a_r(w̅_k,Δ^(r) + w̅_k,Π^(r), w̅_k,Δ^(r) + w̅_k,Π^(r))
    
    
       = Θ |w̅_k,Δ^(r) + w̅_k,Π^(r)|_a_r^2

Then (<ref>) holds. 

Further, utilizing the linear operator T_F_k (or matrix T⃗_F_k) defined above, we can transform the basis function vector Φ^k of W_k into the so-called new dual-primal basis function vector

    Φ_T^k = T_F_kΦ^k :=
    ([ Φ^k_Δ; Φ^k_Π ]),

where

    Φ^k_χ = (ϕ^k_χ,1,⋯,ϕ^k_χ,n_χ^k)^T = T_χ^F_kΦ^k:=
     (T⃗_χ^F_k)^T Φ^k, χ=Δ,
    Π.


From this, we can decompose W_k into

    W_k = W_k,Δ⊕  W_k,Π,

where the function spaces W_k,Δ and W_k,Π (k=1,⋯, N_f) are separately formed by the component functions of Φ^k_Δ and Φ^k_Π.

Then, using (<ref>), a decomposition of the space Ŵ can be obtained as follows

    Ŵ = W_Δ⊕ W_Π,

where the dual and primal variable space 

    W_Δ=⊕_k=1^N_f W_k,Δ , W_Π = (⊕_k=1^N_f W_k,Π) ⊕ (⊕_k=1^N_v W_k,c).


















































Similarly, using the linear operator T_F_k (or matrix T⃗_F_k), we can transform the basis function vector Φ^k,ν of
W_k^(ν) into the dual-primal basis function vector

    Φ_T^k,ν = T_F_kΦ^k,ν :=
    ([ Φ^k,ν_Δ; Φ^k,ν_Π ]),

where

    Φ^k,ν_χ = (ϕ^k,ν_χ,1,⋯,ϕ^k,ν_χ,n_χ^k)^T = T_χ^F_kΦ^k,ν:=
     (T⃗_χ^F_k)^T Φ^k,ν, χ=Δ,
    Π.



By using the definitions (<ref>),
(<ref>) of {ϕ^k_χ,l} and {ϕ^k,ν_χ,l} (χ = Δ,Π, ν=i,j), and (<ref>), we get

    ϕ^k_χ,l|_D̅_ν = ϕ^k,ν_χ,l|_D̅_ν, l=1,⋯,n_χ^k, χ = Δ,Π.


Decompose W_k^(ν) into

    W_k^(ν) = W_k,Δ^(ν)⊕
    W_k,Π^(ν), ν=i,j,

where the basis function vector of W_k,Δ^(ν) and W_k,Π^(ν) are Φ^k,ν_Δ and Φ^k,ν_Π, respectively.


Using the above decomposition, the partially coupled function space which is relied on the adaptive BDDC preconditioner can be obtained and expressed as

    W̃ = W̃_Δ⊕ W_Π,

where
W_Π is defined in (<ref>), and

    W̃_Δ = ⊕_r=1^N_d W_Δ^(r),
     W_Δ^(r) =⊕_k∈ℳ_r
    W_k,Δ^(r) = ⊕_l=1^f_r
    W_r_l,Δ^(r), r=1,⋯, N_d,

here ℳ_r is defined in (<ref>), f_r denotes the size of ℳ_r, and the subspace W_r_l, Δ^(r) is defined in (<ref>).





















































Further, for the need of the theoretical analysis, we can present the decomposition of the auxiliary space W̅_k^(ν) based on the dual-primal basis function vector as

    W̅_k^(ν) = W̅_k,Δ^(ν)⊕W̅_k,Π^(ν),

where the basis function vectors Φ̅_Δ^k,ν and Φ̅_Π^k,ν of the subspaces W̅_k,Δ^(ν) and
W̅_k,Π^(ν) are defined in (<ref>), respectively.

In the following, the dual-primal basis function vectors {Φ_Δ^k} and {Φ_Π^k} for the Schur complement space Ŵ will be adopted.
And we will design and analyze the adaptive BDDC preconditioner for the corresponding Schur complement system  (<ref>).








 §.§ BDDC preconditioner


We focus on the two-level adaptive BDDC preconditioner firstly.
In order to describe this preconditioner, we need to introduce some commonly used linear operators firstly.

Let R_U^V be the restriction operator from the Hilbert space U onto its subspaces V,
and I_V^U be the interpolation operator from V to U (see reference <cit.>).
In particular, when V=U, I_U^U (or I_V^V) is an identity operator.

Denote the Hilbert spaces Z = span{ϕ_1^Z,⋯,ϕ_n^Z} and W = span{ϕ_1^W,⋯,ϕ_n^W}, define the linear basis transformation operator T_Z^W: Z→ W such that

    T_Z^W ϕ^Z_l =ϕ^W_l, l=1,⋯, n.

In particular, for any k=1, ⋯, N_f and ν=r, j, we have

    T_W_k^W_k^(ν)ϕ^k_χ, l= ϕ^k,
    ν_χ,l, 
    T^W_k_W_k^(ν)ϕ^k, ν_χ, l =
    ϕ^k_χ, l, 
    l=1,⋯, n_χ^k, χ = Δ,Π,
           T_W_k^(r)^W_k^(j)ϕ^k,r_χ,l= ϕ^k,j_χ,l, T_W_k^(ν)^W̅_k^(ν)ϕ^k,ν_χ,l = ϕ̅^k,ν_χ,l, 
    l=1,⋯, n_χ^k, χ = Δ,Π,

and for any k =1,⋯,N_v and r ∈𝒩_k, we have

    T_W_k,c^W_k,c^(r)ψ^k_l= ψ^k,r_l, l=1,⋯,n_k^c.








For a given linear operator L from the Hilbert space U to the Hilbert space V, the complex conjugate transpose
operator L^H: V → U is defined by

    (L^H v, u) = (v, Lu),  ∀ u ∈ U, v∈ V.



Then, using the basis transformation operators T_W_k^W_k^(r) and T_W_k,c^W_k,c^(r), another sesquilinear form on W̃ can be introduced.

For any ũ,ṽ∈W̃, using the decomposition (<ref>), (<ref>) and (<ref>) of W̃, we have

    ζ̃ = ∑_r=1^N_d∑_k∈ℳ_rζ̃_k,Δ^(r) + ∑_k=1^N_fζ̃_k,Π + ∑_k=1^N_vζ̃_k,c,  ζ̃ = ũ, ṽ,

where ζ̃_k,Δ^(r)∈ W_k,Δ^(r),
ζ̃_k,Π∈ W_k,Π and ζ̃_k,c∈ W_k,c.
From this, we can define a sesquilinear form
ã(·,·) and its corresponding HPD operator S̃: W̃→W̃ such that

    (S̃ũ, ṽ) :=ã(ũ,
    ṽ) = ∑_r=1^N_d a_r(ũ^(r),
    ṽ^(r)),  ∀ũ, ṽ∈W̃,

where

    ζ̃^(r) = ∑_k ∈ℳ_r
    (ζ̃_k,Δ^(r) + T_W_k^W_k^(r)ζ̃_k,Π) + ∑_k∈ℳ_r^cT_W_k,c^W_k,c^(r)ζ̃_k,c, ζ̃ = ũ, ṽ.




For any r=1,⋯, N_d, we define a linear operator Î_W_Δ^(r)^D: W_Δ^(r)→Ŵ such that

    Î_W_Δ^(r)^D
    =∑_k ∈ℳ_r T^W_k_W_k^(r)
    D_F_k^(r) R_W_Δ^(r)^W_k,Δ^(r),

where R_W_Δ^(r)^W_k,Δ^(r) is a restriction operator from W_Δ^(r) to its subspace W_k,Δ^(r).
It's easy to verify that for any given k ∈ℳ_r, we can obtain

    Î_W_Δ^(r)^D w =  T^W_k_W_k^(r)
    D_F_k^(r) w, ∀ w  ∈ W_k,Δ^(r).









Using Î_W_Δ^(r)^D (r=1,⋯,N_d) and the restriction operators R_W̃^W_Δ^(r) (r=1,⋯,N_d) and R_W̃^W_Π, we can introduce an average operator E_D: W̃→Ŵ, which satisfy

    E_D = Q_D+ R_W̃^W_Π,

where

    Q_D=∑_r=1^N_dÎ_W_Δ^(r)^D
    R_W̃^W_Δ^(r).


With the above-mentioned preparations, by using the sesquilinear form ã(·,·) and the average operator E_D,
the adaptive BDDC operator M_BDDC^-1: Ŵ→Ŵ for solving the Schur
complement system (<ref>) can then be given in the following algorithm.


From this algorithm, using the definition (<ref>) of S̃ and note that S̃ is Hermitian positive definite, it is easy to verify that M_BDDC^-1 can be written as

    M_BDDC^-1 = E_D S̃^-1 (E_D)^H.



In order to facilitate parallel programming, we can give an equivalent description of Algorithm <ref>.
For this purpose, some operators are introduced firstly.


Using the prolongation operators I_W_Δ^(r)^W̃ (r=1,⋯,N_d), a linear operator from W̃ to W̃_Δ can be defined as

    S̃_Δ^-1= ∑_r=1^N_d
    (S̃^(r)_ΔΔ)^-1(I_W_Δ^(r)^W̃)^H = ∑_r=1^N_d
    I_W_Δ^(r)^W̃
    (S̃^(r)_ΔΔ)^-1(I_W_Δ^(r)^W̃)^H,

where

    S̃^(r)_ΔΔ =
    (I_W_Δ^(r)^W̃)^H S̃
    I_W_Δ^(r)^W̃.


Further, utilizing S̃_Δ^-1, we can introduce linear operator O_Π̃: W_Π→W̃ as

    O_Π̃= I_W_Π^W̃ - S̃_Δ^-1S̃
    I_W_Π^W̃ = (I_W̃^W̃ - S̃_Δ^-1S̃) I_W_Π^W̃,

where I_W_Π^W̃ is the prolongation operator from W_Π to W̃, I_W̃^W̃ is an identity operator on W̃.


Thus, by using the expression (<ref>) of the adaptive BDDC operator M_BDDC^-1,
refer to the derivation process of Theorem 4.1 in <cit.>, we can see that Algorithm <ref> can be described as equivalent as follows.


Since Algorithm <ref> can be considered as a two-level algorithm, we also call Algorithm <ref> or Algorithm <ref> two-level adaptive BDDC algorithm.


























































































































































































































































































































Further, by using Algorithm <ref> or Algorithm <ref>,
we can get the following algorithm for solving the original variational problem (<ref>).








From the results of numerical experiments for the PWLS discretizations of the Helmholtz equations, we find that
the number of primal unknowns increases as the number of wave numbers and subdomains increase, and the corresponding coarse problem will become too large and hard to solve directly. This leads to multilevel extension of this algorithm naturally. Multilevel BDDC algorithm were first presented in <cit.>, and further developed in <cit.>. In addition, the multilevel preconditioners for solving the systems arising from the plane wave discretizations for Helmholtz equations with large wave numbers were constructed by Hu and Li in <cit.>.

Following <cit.>, in our multilevel algorithm, we denote the sth level mesh by 𝒯_h^s (s=0,⋯,L-1),
where L is the total number of levels, 𝒯_h^0 = 𝒯_h is the finest level,
and a subdomain at a finer level is considered as an element of a coarser mesh.
Let V_p^s(𝒯_h^s) be the finite element spaces of the original problem associated with 𝒯_h^s,
and set V_p^s+1(𝒯_h^s+1) := W_Π^s,
where W_Π^s is the coarse space at level s.
In particular, V_p^0(𝒯_h^0) = V_p(𝒯_h) and V_p^1(𝒯_h^1) = W_Π are defined in (<ref>) and (<ref>), respectively. Noticing that Algorithm <ref> gives an iteration process from level s to level s+1 (s=0), and by using the solution of the Schur complete problem, the solution of the original problem at level 0 can be obtained.
Iterating this procedure until s < L-1, and we can arrive at our multilevel adaptive BDDC algorithm.
A relation between the condition number of the multilevel BDDC algorithm with corner coarse function
in 2D and each level problem was presented in <cit.>.


In the next section, we will provide the condition number estimate of the two-level adaptive BDDC preconditioned operator.




§ ANALYSIS OF CONDITION NUMBER BOUND



We first establish the relation between the operators S̃ and Ŝ, which are defined in (<ref>) and (<ref>), respectively. For this purpose, we introduce a subspace of W̃ such that
 
    W̅̃̅ =W̅̃̅_Δ⊕ W_Π,

where W_Π is defined in (<ref>),  and

    W̅̃̅_Δ = ⊕_k=1^N_fW̅̃̅_k,Δ,  W̅̃̅_k,Δ = span{ϕ^k,r_Δ,1
    +ϕ^k,j_Δ,1,⋯,ϕ^k,r_Δ,n_Δ^k
    +ϕ^k,j_Δ,n_Δ^k},

here {ϕ^k, ν_Δ, l} is defined in (<ref>).

Using (<ref>), we can define another linear basis transformation operator T_Ŵ^W̅̃̅: Ŵ→W̅̃̅, which satisfies that for any k=1, ⋯, N_f

    T_Ŵ^W̅̃̅ϕ^k_Π, l=
    ϕ^k_Π, l, l=1,⋯,
    n_Π^k; T_Ŵ^W̅̃̅ϕ^k_Δ, l=
    ϕ^k,r_Δ,l
    +ϕ^k,j_Δ,l, l=1,⋯, n_Δ^k,

and for any k=1,⋯,N_v

    T_Ŵ^W̅̃̅ψ^k_l=
    ψ^k_l, l=1,⋯,n_k^c,

where {ϕ^k_Δ, l} and {ϕ^k_Π, l} are defined in (<ref>), and {ψ^k_l} are defined in (<ref>).


From (<ref>) and the definition
(<ref>) of the scaling operators
D_F_k^(r) (k ∈ℳ_r, 1 ≤ r ≤ N_d), we can easily prove
that









    T_Ŵ^W̅̃̅ T_W_k^(r)^W_k  D_F_k^(r) w
    = D_F_k^(r)(w + T^W_k^(j)_W_k^(r) w), ∀ w
    ∈ W_k,Δ^(r),

here we have used the assumption that F_k = ∂ D_r ∩∂ D_j.

It follows from (<ref>), (<ref>), (<ref>) and (<ref>) that

    Ŝ = (T_Ŵ^W̅̃̅)^H S̃
     T_Ŵ^W̅̃̅.


Combining (<ref>) and (<ref>), we can obtain the preconditioned operator as

    Ĝ =M_BDDC^-1Ŝ = E_D S̃^-1 (E_D)^H
    (T_Ŵ^W̅̃̅)^H S̃
    T_Ŵ^W̅̃̅.



In the following, we will derive the upper bound of the condition number of Ĝ.

For any w̃∈W̃, using (<ref>), we have

    w̃ = ∑_r=1^N_d∑_k ∈ℳ_r w_k,Δ^(r) + w_Π,  w_Π := ∑_k=1^N_f w_k,Π + ∑_k=1^N_v w_k,c,

where

    w_k,Δ^(r) = (w⃗_k,Δ^(r))^T Φ_Δ^k,r∈ W_k,Δ^(r),  w_k,Π = (w⃗_k,Π)^T Φ_Π^k∈ W_k,Π,  w_k,c = (w⃗_k,c)^T Ψ_Π^k∈ W_k,c,

here w⃗_k,Δ^(r)∈ℂ^n_Δ^k, w⃗_k,Π∈ℂ^n_Π^k and w⃗_k,c∈ℂ^n_k^c.


Following Theorem 1 in <cit.>, we need to estimate the bound

    ã(P_D w̃, P_D w̃) ≤ C ã(w̃, w̃), ∀w̃∈W̃,

where P_D: W̃→W̃ is a jump operator defined as

    P_D = I_W̃^W̃ - T_Ŵ^W̅̃̅ E_D,

here I_W̃^W̃ is the identity operator on W̃, the linear basis transformation operator T_Ŵ^W̅̃̅ and the average operator E_D are separately defined in (<ref>) and (<ref>).





For any w̃∈W̃, we have

    P_D w̃
    = ∑_r=1^N_d∑_k ∈ℳ_r
    (w_D,Δ^k,r - w̃_D,Δ^k,r),

where

    w_D,Δ^k,r:= (w⃗_k,Δ^(r))^T
    Φ_D,Δ^k,r,  w̃_D,Δ^k,r :=
    (w⃗_k,Δ^(j))^T Φ_D,Δ^k,r,

here the basis function vector Φ_D,Δ^k,r is defined in (<ref>).



By the definition (<ref>) of P_D, (<ref>) of E_D, (<ref>) of T_Ŵ^W̅̃̅, and using the decomposition (<ref>),
P_D w̃ can be rewritten as follows:

    P_D w̃   =   
    w̃ - T_Ŵ^W̅̃̅E_D w̃
    
       =   w̃ - T_Ŵ^W̅̃̅(Q_D +  R_W̃^W_Π) w̃
    
       =    (∑_r=1^N_d∑_k ∈ℳ_r w_k,Δ^(r) + w_Π) - (T_Ŵ^W̅̃̅Q_D w̃ + w_Π)
    
    
       =   ∑_r=1^N_d∑_k ∈ℳ_r w_k,Δ^(r)
       - T_Ŵ^W̅̃̅Q_D w̃.



By using (<ref>), (<ref>), (<ref>)
 and (<ref>), we find the second term in
(<ref>) satisfies (F_k = ∂ D_r ∩∂ D_j)

    T_Ŵ^W̅̃̅ Q_D w̃   =    T_Ŵ^W̅̃̅∑_r=1^N_dÎ_W_Δ^(r)^D ∑_k
    ∈ℳ_r w_k,Δ^(r)
       =   ∑_r=1^N_d∑_k ∈ℳ_r
    T_Ŵ^W̅̃̅  T^W_k_W_k^(r) D_F_k^(r)
    w_k,Δ^(r)
       =   ∑_r=1^N_d∑_k ∈ℳ_r
    D_F_k^(r)(w_k,Δ^(r) + T_W_k^(r)^W_k^(j)
    w_k,Δ^(r)).



Combining (<ref>) and (<ref>),  and using
(<ref>), (<ref>), (<ref>),
(<ref>), (<ref>) and
(<ref>), we have

    P_D w̃   =   ∑_r=1^N_d∑_k ∈ℳ_r w_k,Δ^(r) - ∑_r=1^N_d∑_k ∈ℳ_r D_F_k^(r)(w_k,Δ^(r)
    + T_W_k^(r)^W_k^(j) w_k,Δ^(r))
       =   ∑_r=1^N_d∑_k ∈ℳ_r
    (D_F_k^(j)w_k,Δ^(r) - D_F_k^(r)
    T_W_k^(r)^W_k^(j) w_k,Δ^(r))
       =   ∑_k=1^N_f((D_F_k^(j)
    w_k,Δ^(r) - D_F_k^(r) T_W_k^(r)^W_k^(j)
    w_k,Δ^(r)) + (D_F_k^(r) w_k,Δ^(j) -
    D_F_k^(j) T_W_k^(j)^W_k^(r)
    w_k,Δ^(j)))
       =   ∑_r=1^N_d∑_k ∈ℳ_rD_F_k^(j)(w_k,Δ^(r) -
    T_W_k^(j)^W_k^(r) w_k,Δ^(j))
       =   ∑_r=1^N_d∑_k ∈ℳ_r
    (w_D,Δ^k,r - w̃_D,Δ^k,r)

where w_D,Δ^k,r and w̃_D,Δ^k,r are defined in (<ref>). That is to say (<ref>) holds.



For a given tolerance Θ≥ 1, we can obtain the following estimate for w̃∈W̃

    ã(P_D w̃, P_D w̃) ≤ C Θã(w̃, w̃),

where C = 2 C_F^2 and C_F = max_r{f_r}, here f_r denotes the number of interface on ∂ D_r. 



In view of the definition (<ref>) of the sesquilinear form ã(·,·), (<ref>) of the semi-norm |·|_a_r, and the decompositions (<ref>) and (<ref>), it is equivalent to show that

    ∑_r=1^N_d |(P_D w̃)^(r)|_a_r^2 ≤ C
    Θ∑_r=1^N_d |w̃^(r)|_a_r^2,

where

    w̃^(r) = ∑_k ∈ℳ_r (w_k,Δ^(r) + w_k,Π^(r)) + ∑_k∈ℳ_r^c w_k,c^(r),  (P_D w̃)^(r) =  ∑_k ∈ℳ_r (w_D,Δ^k,r - w̃_D,Δ^k,r),

here

    w_k,Π^(r) = (w⃗_k,Π)^T Φ_Π^k,r∈ W_k,Π^(r),  w_k,c^(r) = (w⃗_k,c)^T Ψ^k,r∈ W_k,c^(r).




By (<ref>),

C_F = max_r{f_r},
the conditions  (<ref>) associate with T_F_k (k∈ℳ_r, 1≤ r≤ N_d), and  (<ref>), we have 

    ∑_r=1^N_d |(P_Dw̃)^(r)|_a_r^2
       =   ∑_r=1^N_d |∑_k ∈ℳ_r(w_D,Δ^k,r -
    w̃_D,Δ^k,r)|_a_r^2
       ≤    2 C_F
    ∑_r=1^N_d∑_k ∈ℳ_r(|w_D,Δ^k,r|_a_r^2 +
    |w̃_D,Δ^k,r|_a_r^2 )
       =    2 C_F ∑_r=1^N_d∑_k ∈ℳ_r(|w_D,Δ^k,r|_a_r^2 +
    |w̃_D,Δ^k,j|_a_j^2)
       ≤    2 C_F Θ∑_r=1^N_d∑_k ∈ℳ_r |w̅_k,Δ^(r) +
    w̅_k,Π^(r)|_a_r^2
       ≤    2 C_F^2 Θ∑_r=1^N_d
    |w̃^(r)|_a_r^2

where

    w̃_D,Δ^k,j := (w⃗_k,Δ^(r))^T Φ_D,Δ^k,j,  w̅_k,Δ^(r) = (w⃗_k,Δ^(r))^T Φ̅_Δ^k,r∈W̅_k,Δ^(r),  w̅_k,Π^(r)  = (w⃗_k,Π)^T Φ̅_Π^k,r∈W̅_k,Π^(r).



Finally, the estimate (<ref>) follows from (<ref>).



By Lemma <ref> and following Theorem 1 in <cit.>,  we obtain:

For a given tolerance Θ≥ 1, the condition number bound of the two-level adaptive BDDC preconditioned systems Ĝ

    κ(Ĝ) ≤ C Θ,

where  C is a constant which is just depending on the maximum number of interfaces per subdomain.






§ NUMERICAL RESULTS


In this section, we present examples with constant and variable wave number κ to investigate the convergence properties of the adaptive BDDC preconditioners proposed in this paper.

In the following numerical tests, the given region Ω is divided into N_h quadrilateral mesh 𝒯_h,
where h denotes the size of elements.
Let p denote the number of plane wave bases in each element. Therefore, it is easy to known that the number of dofs
in the plane wave space V_p(𝒯_h) is N_h × p.

Since the special interface is needed in this paper, and to guarantee the load balance,
we decompose 𝒯_h into some subdomains which satisfy that the number of complete elements in each subdomain is the same and denotes as n.
For example if we set N_h = 35^2 and the number of subdomains N_d = 4^2, then n = 8^2.


Since the stiffness matrix of the PWLS method is HPD,
we solve the Schur complement system (<ref>) by PCG method, and the iteration is stopped either the relative residual is reduced by the factor of 10^-5 or the iteration counts are greater than 100. These algorithms are implemented using Matlab
and run in a machine with Intel(R) Xeon(R) CPU E5-2650 v2 2.60 GHz and 96 GB memory.

Numerical experiment results show that the modulus of the eigenvalue λ of the generalized eigenvalue problems (<ref>) satisfies |λ| ≥ 1, hence, in our adaptive BDDC algorithms, we set the tolerance Θ = 1+ log(min_1 ≤ r ≤ N_d{n_x^(r), n_y^(r)}) for a given mesh partition,
where n_x^(r), n_y^(r) are separately the number of complete and part elements in x and y direction of subdomain D^(r) (r=1,⋯,N_d).
The transform matrix T⃗_F_k in each interface F_k is defined in (<ref>).
Therefore, the algorithm is uniquely determined by the scaling matrices D⃗_F_k^(ν)(ν=r,j) for each interface.
In the following experiments, we separately denote the algorithm with D⃗_F_k^(ν)(ν=r,j, k=1,⋯,N_f) defined in (<ref>) and (<ref>) as method1 and method2.
We will apply both methods to three typical examples to investigate the scalability of these methods
measured by the mesh size, number of subdomains and angular frequency (or wave number).
Since the wave number κ = ω/c, where ω and c are separately the angular frequency and the wave speed,
we can react the variousness of the wave number to the angular frequency and the wave speed.


In all of the tables in this section,  iter is the number of iterations for the PCG algorithms, λ_min and λ_max separately denote the minimum  and maximum eigenvalues of the preconditioned system,  pnum is the number of primal unknowns, the average number of primal unknowns per interface are given in the parentheses, and the proportion of the total number of primal unknowns to the total number of interface dofs are denoted as ppnum.

Before numerical studies are performed, three typical examples are given here.

(Constant medium)
<cit.> Consider model problem (<ref>), where Ω = (0,2) × (0,1), the wave speed c ≡ 1, the exact solution of the problem can be expressed as

    u_ex = cos(12 π y)(A_1 e^-i ω_x x + A_2 e^i ω_x x),

here ω_x = √(ω^2 - (12 π)^2), and coefficients A_1 and A_2 satisfy the equation

    (
    [                ω_x               -ω_x; (ω - ω_x)e^-2i ω_x  (ω + ω_x)e^2i ω_x ]) ([ A_1; A_2 ])
    = ([ -i;  0 ]).



(Piecewise constant medium)
Consider model problem (<ref>), which is a variant of the Marmousi model in <cit.> or <cit.>,
where Ω = (0,7200)×(0,3600), and the wave speed c is defined by (see figure <ref> for illustration)

    c(x,y) = {[            1800    y ∈ [0,1200];            3600 y ∈ [1200,2400];            5400 y ∈ [2400,3600] ].  x∈ [0,7200],

and g = x^2 + y^2.



(Random medium)
Consider Example <ref>, where c(x,y)  is chosen randomly from [1500,5500] for each grid element, as shown in Figure <ref>,
which can be seen as a complex version of the Marmousi model in <cit.>.






 §.§ Two-level results




  §.§.§ Scalability study with respect to angular frequency


In this subsection, we investigate the influence of the angular frequency (or the wave number) to the efficiency of our adaptive BDDC algorithms.
This numerical studies are carried on a N_d = 4^2 subdomain partition with different frequency ω.
The relative L^2 error presented in Table <ref> is defined as

    u_ex - u_h_L^2(Ω)/u_ex_L^2(Ω),

where u_ex and u_h are separately the exact solution and the numerical solution of this problem.
We kept ω h as a constant and slightly increase p or decrease h to control the relative error less than 5× 10^-3.
In this article, since the errors of the approximate solutions are not our main interest, we only list the relative L^2 errors in Table <ref>
for Example <ref>.

As seen in Table <ref>, with the same stopping criterion, the relative L^2 errors of the approximate solutions of Example <ref> generated by method1 and method2 do not have large differences.











































































































The dependence of the iteration counts on the angular frequency is presented in Table <ref>, Table <ref> and Table <ref> for various examples, and a weak dependency relationship is clearly shown.
In addition, from these tables, we can see that the total number of primal unknowns (pnum) increased with the increase of the angular frequency for both method1 and method2, especially when ω = 320π in Example <ref>, the dofs of the coarse problem has reached more than thirty thousand in method1.
The number of primal unknowns also highlights the superiority of method2
over method1.













  §.§.§ Scalability study with respect to number of subdomains


In the second subsection, we set the number of the complete elements in each subdomain n is 8^2, p=10 and
the frequency to ω = 20π. The effect of number of the subdomains is explored
here by adding N_d from 3^2 to 6^2.

Results of these experiments are shown in Table <ref>, Table <ref>
and Table <ref>. From these tables, we can see that the iteration counts are almost independent on the number of subdomains, and as the number of subdomains increases, the total number of primal unknowns increase, but the average number of primal unknowns per interface almost remains the same.
Apart from this, one can see that for these three examples, method2 still has an absolute advantage in the number of the primal dofs over method1,
in particular, less than half of the primal dofs are required for method2 to ensure the convergence stability.









  §.§.§ Scalability study with respect to mesh size


Finally, in this subsection, we consider the scalability with respect to mesh size, where p=9, the number of the subdomains and the angular frequency are both held constant (N_d = 4^2, ω = 20π), but the number of complete elements in each subdomain varied from n=6^2 to n=24^2.

The iteration efficiency as a function of the mesh size are reported in Table <ref>, Table <ref> and Table <ref> for method1 and method2. We can see from these tables that the iteration counts are mildly dependent on the mesh size. The minimum eigenvalues of the preconditioned systems are larger than 1 and the maximum eigenvalues are mildly dependent on the mesh size.
As the mesh size decrease, there is a great difference between method1 and method2 in the number of the primal unknowns even for constant medium, which is different from the adaptive BDDC algorithms for the two-order elliptic problems.

























































The results listed in the above-mentioned tables have verified
the correctness of the theoretical results,
and indicate that the proposed two-level adaptive BDDC algorithms are very effective for solving Helmholtz equation with large and various wave numbers.
However, one fact worth noticing is that the number of primal unknowns increase as the wave numbers or the number of subdomains increase,
the corresponding coarse problem will become too large and hard to solve directly, and it is more obvious for the multiplicity scaling matrices.































 §.§ Multilevel results





In this subsection, we try to apply the multilevel adaptive BDDC algorithm with multiplicity scaling matrices to Example <ref> to show its performance for large wave numbers.

For this test, the PCG algorithm is stopped either the iteration counts are greater than 100 or the relative residual is reduced by the factor of 10^-5 at level 0 and 10^-2 at other levels. We consider a higher wave number model and set the angular frequency ω = 320π, the number of plane wave bases in each element p=15.
In the algorithm, four subdomains at the finer level are treated as a coarser subdomain.








Table <ref> shows the performance of our multilevel adaptive BDDC algorithm for different N_d and n, where Tdofs is the total number of dofs, Fpnum and Cpnum separately denote the number of dofs at level 1 and the coarsest level.
From this table, we can see that the number of dofs at level 1 (Fpnum) increases as the number of subdomains increases,
and the scale of the coarse problem which need to be compute directly are hard to be accepted when we use a two-level adaptive BDDC algorithm.
With the increase of the number of levels (L), the number of dofs at the coarsest level (Cpnum) is reduced to about 30%, and the iteration counts remains almost the same. That is to say, the multilevel algorithm is effective for reducing the number of dofs at the coarse problem and it is also efficient for solving large wave number problems.














§ CONCLUSIONS


In this paper, by introducing some auxiliary spaces, dual-primal basis functions and operators with essential
properties, BDDC algorithms with adaptive primal unknowns are developed and analyzed for the PWLS discritizations of the Helmholtz equations with high and various wave number.

Since the dofs of the PWLS discritization are defined on elements rather than vertices or edges,
we introduce a special “interface" and the corresponding sesquilinear form for each subdomains,
which is different from the interface in the existing BDDC algorithms.

As the eigenvalues of the local generalized eigenvalue problems are complex, we choose the primal constraints which are formed by the eigenvectors with
their complex modulus of eigenvalues greater than a user-defined tolerance value Θ,
which is used to construct the transformation operators for selecting dual-primal basis functions and control the condition number.

We prove that the condition number of the preconditioned system is bounded above by C Θ.

A multilevel algorithm is attempted to resolve the bottleneck in large scale coarse problem.
Numerical results are presented to verify the robustness and efficiency of the proposed approaches.
Further, we will extend this method to three-dimensional
case. There is fundamentally different from two-dimensional case and three-dimensional case,
such as we need to consider the construction of the generalized eigenvalue problem for each edge which is shared by more than two subdomains.



§ ACKNOWLEDGMENTS


The authors would like to thank Professor Qiya Hu for his suggestions in the PWLS discritization.
This work is supported by the National Natural Science Foundation of China (Grant Nos. 11571293, 11201398, 11601462), Hunan Provincial Natural
Science Foundation of China (Grant No. 2016JJ2129), General Project of Hunan Provincial Education Department of China
(Grant No. 17C1527), Open Foundation of Guangdong Provincial Engineering Technology Research Center for Data Science(Grant No. 2016KF07), and Hunan Provincial Civil-Military Integration Industrial Development Project “Adaptive Multilevel
Solver and Its Application in ICF Numerical Simulation".



§ REFERENCES





elsarticle-num



99

HMP2016
Hiptmair R, Moiola A, Perugia I. A survey of Trefftz methods for the Helmholtz equation, in the book “Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations" (Lecture Notes in Computational Science and
Engineering, Volume 114), Editors: G. R. Barrenechea, F. Brezzi, A. Cangiani and E. H.
Georgoulis, Springer, 2016; 237-279.


FHH2003
Farhat C, Harari I, Hetmaniuk U.
A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the mid-frequency regime.
Computer Methods in Applied Mechanics and Engineering 2003; 192:1389-1419.

FWT2004
Farhat C, Wiedemann-Goiran P, Tezaur R.
A discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured meshes.
Wave Motion 2004; 39(4):307-317.










CD1998
Cessenat O, Després B.
Application of an ultra weak variational formulation of elliptic PDEs to the
two-dimensional Helmholtz problem.
SIAM Journal on Numerical Analysis 1998; 35(1):255-299.


HKM2004
Huttunen T, Kaipio JP, Monk P.
The perfectly matched layer for the ultra weak variational formulation of the 3D Helmholtz equation.
International Journal for Numerical Methods in Engineering 2004; 61(7):1072-1092.


L1996
Ladevèze P.
A new computational approach for structure vibrations in the medium frequency range.
Comptes Rendus de l'Académie des Sciences de Paris, Sér. II 1996; 322(12):849-856.

LA2001
Ladevèze P, Arnaud L, Rouch P, Blanzé C.
The variational theory of complex rays for the calculation of medium-frequency vibrations.
Eegineering Computations 2001; 18193-214.


BS1997
Babuška IM, Sauter SA.
Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?
SIAM Journal on Numerical Analysis 1997; 34(6):2392-2423.

DB1999
Deraemaeker A, Babuška I, Bouillard P.
Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions.
International journal for numerical methods in engineering 1999; 46(4):471-499.




TP1995
Thompson LL, Pinsky PM.
A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation.
International Journal for Numerical Methods in Engineering 1995; 38(3):371-397.

S1998
Stojek M.
Least-squares Trefftz-type elements for the Helmholtz equation.
International Journal for Numerical Methods in Engineering 1998; 41(5):831-849.


MW1999
Monk P, Wang DQ.
A least-squares method for the Helmholtz equation.
Computer Methods in Applied Mechanics and Engineering 1999;  175:121-136.

HY2014
Hu Q, Yuan L.
A weighted variational formulation based on plane wave basis for discretization of Helmholtz equations.
International Journal of Numerical Analysis and Modeling 2014; 11(3):587-607.


























HL20171
Hu Q, Li X.
Novel multilevel preconditioners for the systems arising from plane wave discretization of Helmholtz equations with large wave numbers.
SIAM Journal on Scientific Computing 2017; 39(4):A1675-A1709.

HL20172
Hu Q, Li X.
Efficient Multilevel Preconditioners for Three-Dimensional Plane Wave Helmholtz Systems with Large Wave Numbers.
Multiscale Modeling & Simulation 2017; 15(3):1242-1266.


CW1992
Cai XC, Widlund OB.
Domain decomposition algorithms for indefinite elliptic problems.
SIAM Journal on Scientific and Statistical Computing 1992; 13(1):243-258.



LX2014
Liu Y, Xu X.
A Robin-type Domain Decomposition Method with Red-Black Partition.
SIAM Journal on Numerical Analysis 2014; 52(5):2381-2399.

CLX2016
Chen W, Liu Y, Xu X.
A robust domain decomposition method for the Helmholtz equation with high wave number.
ESAIM: Mathematical Modelling and Numerical Analysis 2016; 50(3):921-944.

HZ2016
Hu Q, Zhang H.
Substructuring preconditioners for the systems arising from plane wave discretization of Helmholtz equations.
SIAM Journal on Scientific Computing 2016; 38(4):A2232-A2261.




FMT1999
Farhat C, Macedo A, Tezaur R.
FETI-H: A scalable domain decomposition method for high frequency exterior Helmholtz problems.
Eleventh International Conference on Domain Decomposition Method 1999; 231-241.

TMF2001
Tezaur R, Macedo A, Farhat C.
Iterative solution of large-scale acoustic scattering problems with multiple right hand-sides by a domain decomposition method with Lagrange multipliers.
International Journal for Numerical Methods in Engineering 2001; 51(10):1175-1193.


M2002
Mandel J.
An iterative substructuring method for coupled fluid-solid acoustic problems.
Journal of Computational Physics 2002; 177(1):95-116.

FAT2005
Farhat C, Avery P, Tezaur R, Li J.
FETI-DPH: a dual-primal domain decomposition method for acoustic scattering.
Journal of Computational Acoustics 2005; 13(3):499-524.


Dohrmann2003
Dohrmann CR.
A preconditioner for substructuring based on constrained energy minimization.
SIAM Journal on Scientific Computing 2003; 25:246-258.


LW2006
Li J, Widlund OB.
FETI-DP, BDDC, and block Cholesky methods.
International Journal for Numerical Methods in Engineering 2006; 66:250-271.

BS2007
Brenner SC, Sung LY.
BDDC and FETI-DP without matrices or vectors.
Computer Methods in Applied Mechanics and Engineering 2007; 196:1429-1435.


LT2008
Li J, Tu X.
Convergence analysis of a balalncing domain decomposition method for solving interior Helmholtz equations.
Lawrence Berkeley National Laboratory 2009.

TL2009
Tu X, Li J.
BDDC for nonsymmetric positive definite and symmetric indefinite problems.
Domain Decomposition Methods in Science and Engineering XVIII 2009; 70:75.

MS2007
Mandel J, Sousedík B.
Adaptive selection of face coarse degrees of freedom in the BDDC and FETI-DP iterative substructuring methods.
Computer Methods in Applied Mechanics and Engineering 2007; 196(8):1389-1399.



DP2012
Dohrmann CR, Pechstein C.
Constraint and weight selection algorithms for BDDC.
No. SAND2012-5038C. Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States), 2012.





PD2013
Pechstein C, Dohrmann CR.
Modern domain decomposition solvers-BBDC, deluxe scaling, and an algebraic approach.
Talk by Pechstein at RICAM, Linz, Austria. 2013.

DSO2015
Oh DS, Widlund OB, Zampini S, Dohrmann CR.
BDDC algorithms with deluxe scaling and adaptive selection of primal constraints for Raviart-Thomas vector fields.
Mathematics of Computation 2018; 87(310):659-692.

KC2015
Kim HH, Chung ET.
A BDDC algorithm with enriched coarse spaces for two-dimensional elliptic problems with oscillatory and high contrast coefficients.
Multiscale Modeling & Simulation 2015; 13(2):571-593.

KCW2017
Kim HH, Chung ET, Wang J.
BDDC and FETI-DP preconditioners with adaptive coarse spaces for three-dimensional elliptic problems with oscillatory and
high contrast coefficients.
Journal of Computational Physics 2017;  349:191-214.

Z2016
Zampini S.
PCBDDC: a class of robust dual-primal methods in PETSc.
SIAM Journal on Scientific Computing 2016; 38(5):S282-S306.


KRR2016
Klawonn A, Radtke P, Rheinbach O.
Adaptive coarse spaces for BDDC with a transformation of basis.
Domain Decomposition Methods in Science and Engineering XXII.
Springer International Publishing, 2016; 301-309.

JGC2016
Calvo JG, Widlund OB.
An adaptive choice of primal constraints for BDDC domain decomposition algorithms.
Electronic Transactions on Numerical Analysis 2016; 45:524-544.  

PC2016
Pechstein C, Dohrmann CR.
A unified framework for adaptive BDDC.
Electronic Transactions on Numerical Analysis 2017; 46:273-336.


KCX2017
Kim  HH, Chung ET, Xu C.
A BDDC algorithm with adaptive primal constraints for staggered discontinuous Galerkin approximation of elliptic problems with highly oscillating coefficients.
Journal of Computational and Applied Mathematics 2017;  311:599-617.


DV2017
Da Veiga LB, Pavarino LF, Scacchi S, Widlund OB, Zampini S.
Adaptive selection of primal constraints for isogeometric BDDC deluxe preconditioners.
SIAM Journal on Scientific Computing 2017; 39(1):A281-A302.

PSW2017
Peng J, Shu S, Wang J.
An adaptive BDDC algorithm in variational form for mortar discretizations.
Journal of Computational and Applied Mathematics 2018; 335:185-206.







T20071
Tu X.
Three-level BDDC in three dimensions.
SIAM Journal on Scientific Computing 2007; 29:1759-1780.

T20072
Tu X.
Three-level BDDC in two dimensions.
International Journal for Numerical Methods in Engineering 2007; 69:33-59.

MSD2008
Mandel J, Sousedík B, Dohrmann CR.
Multispace and multilevel BDDC.
Computing 2008; 83:55-85.

ZT2017
Zampini S, Tu X.
Multilevel balancing domain decomposition by constraints deluxe algorithms with adaptive coarse spaces for flow in porous media.
SIAM Journal on Scientific Computing 2017; 39(4):A1389-A1415.

S2013
Stolk CC.
A rapidly converging domain decomposition method for the Helmholtz equation.
Journal of Computational Physics 2013; 241:240-252.

RF1999
Rixen DJ, Farhat C.
A simple and efficient extension of a class of substructure based preconditioners to heterogeneous
structural mechanics problems.
International Journal for Numerical Methods in Engineering 1999; 44(4):489-516.

 1969AD
Anderson WN, Duffin RJ.
Series and parallel addition of matrices.
Journal of Mathematical Analysis and Applications 1969; 26(3):576-594.






































































