Adaptive and frugal BDDC coarse spaces for virtual element discretizations of a Stokes problem with heterogeneous viscosity
Tommaso Bevilacqua Axel Klawonn Martin Lanser


Adaptive and frugal BDDC coarse spaces for virtual element discretizations of a Stokes problem with heterogeneous viscosity

BDDC coarse spaces for virtual element discretizations of Stokes equations





Tommaso Bevilacqua^1, Axel Klawonn^1,2, Martin Lanser^1,2 ^1Department of Mathematics and Computer Science, Division of Mathematics, University of Cologne, Weyertal 86-90, 50931 Cologne, Germany, tommaso.bevilacqua@uni-koeln.de, axel.klawonn@uni-koeln.de, martin.lanser@uni-koeln.de, url: <https://www.numerik.uni-koeln.de> ^2Center for Data and Simulation Science, University of Cologne, Germany, url: <https://www.cds.uni-koeln.de>






*
    Tommaso Bevilacqua Axel Klawonn Martin Lanser
    September 9, 2024
=================================================




The virtual element method (VEM) is a family of numerical methods to discretize partial differential equations on general polygonal or polyhedral computational grids. However, the resulting linear systems are often ill-conditioned and robust preconditioning techniques are necessary for an iterative solution. Here, a balancing domain decomposition by constraints (BDDC) preconditioner is considered. Techniques to enrich the coarse space of BDDC applied to a Stokes problem with heterogeneous viscosity are proposed. In this framework a comparison between two adaptive techniques and a computationally cheaper heuristic approach is carried out. Numerical results computed on a physically realistic model show that the latter approach in combination with the deluxe scaling is a promising alternative.



§ INTRODUCTION



The virtual element method (VEM) <cit.> is a finite element discretization approach for partial differential equations (PDEs) which can deal with very general polygonal/polyhedral computational grids. The linear systems that arise from  these discretizations of PDEs are then generally worse conditioned than in the case of  standard finite element methods (FEMs) for which recent studies proposed robust BDDC methods <cit.>.
In this work we analyze a Stokes problem with high heterogeneity in the viscosity function;  thus an adequatly enriched coarse space is needed.
In the adaptive framework  considered here, this is done by solving a generalized eigenvalue problem on each subdomain edge and by adding the solutions to the coarse space in an appropriate way.  In our numerical simulations we make use of two coarse spaces analyzed for standard low order FEM in the paper <cit.>, identifying with "First" the approach present in Section 4.5 and with "Second" the one in Section 5.
The first approach, was originally introduced in <cit.> and already successfully used for the  VEM in <cit.>. The second one, extensively used in dual-primal finite element tearing and interconnecting (FETI-DP) and BDDC in <cit.>, has been also recently extended to the  VEM for diffusion and linear elasticity problems in <cit.>.

An alternative approach to enrich the coarse space denoted as frugal, has been introduced in <cit.> and already succesfully used for the VEM for stationary diffusion and linear elasticity in <cit.>. This is a heuristic and cheaper technique, that does not involve the solution of eigenvalue problems. This often allows to construct robust coarse spaces in a computationally efficient way when it is sufficient to approximate the largest, or smallest, eigenvalues depending on the chosen coarse space.

In the present work we extend  the first and second adaptive coarse space approaches as well as the frugal coarse space to the virtual element discretization of a Stokes problem with a heterogeneous viscosity function.




§ CONTINUOUS PROBLEM AND VIRTUAL ELEMENT DISCRETIZATION


Let Ω⊆ℝ^2 be a bounded Lipschitz domain, with Γ = ∂Ω, and consider the stationary Stokes problem with homogeneous Dirichlet boundary conditions: find (𝐮,p) s.t.  

    -ν(𝐱) Δ𝐮 - ∇ p = 𝐟        in Ω
    div 𝐮 = 0          in Ω
    𝐮=0         on Γ.

Here 𝐮 and p are respectively the velocity and the pressure fields, 𝐟∈ [H^-1(Ω)]^2 represents the external force and ν∈ L^∞(Ω), ν(𝐱) >0, 𝐱∈Ω, is the heterogeneous viscosity function. 

Introducing 𝐕:=[H^1_0(Ω)]^2 and Q:=L^2_0(Ω)= { q∈ L^2(Ω)  s.t.  ∫_Ω q =0 }, the standard variational formulation reads: find (𝐮,p) ∈𝐕× Q s.t.

    a(𝐮,𝐯)+b(𝐯,p)=(𝐟,𝐯)    for all 𝐯∈𝐕, 
    
    		b(𝐮,q)=0    for all  q ∈ Q,

where a(𝐯,𝐰) := ∫_Ων(𝐱) ∇ 𝐯 : ∇ 𝐰 for all 𝐯,𝐰∈𝐕, b(𝐯,q) := ∫_Ωdiv 𝐯  q for all 𝐯∈𝐕, q ∈ Q and (𝐟,𝐯):=∫_Ω𝐟·𝐯  for all 𝐯∈𝐕.

The discretization of problem (<ref>) is based on a virtual element space which is designed to solve a Stokes problem. In the following we present the basic elements of this discretization; we refer  to <cit.> for further details.

Let {_h}_h be a sequence of triangulations of Ω into general polygonal elements K with h_K:=diameter(K) and h:=sup_K∈_hh_K.
We suppose that, for all h, each element K∈_h satisfies, for some γ>0 and c>0, the following assumptions

0em 
	
  * K is star-shaped with respect to a ball of radius greater or equal than γ h_k,
	
  * the distance between any two vertices of K is greater or equal than c h_K,
 

For k ∈ℕ, we then define the spaces: ℙ_k(K) the set of polynomials on K of degree smaller or equal than k, 𝔹_k(K):={v ∈ C^0(∂ K )  s.t.  v_|e∈ℙ_k(e)   ∀ edge  e ∈∂ K }, 𝐺_k(K):=∇(ℙ_k+1(K)) ⊆ [ℙ_k(K)]^2 and its L^2-orthogonal complement 𝐺_k(K)^⊥⊆ [ℙ_k(K)]^2.
The local virtual element spaces are defined, for k≥2, on each K ∈𝑇_h as

    𝐕_h^K:={𝐯∈ [H^1(K)]^2   | 𝐯_| ∂ K∈ [𝔹_k(∂ K)]^2,
    		 -νΔ𝐯 - ∇ s ∈𝐺_k-2(K)^⊥, 
    div 𝐯∈ℙ_k-1(K),-20mu s ∈ L^2(K)},
        Q_h^K:=ℙ_k-1(K),

and the global virtual element spaces are	𝐕_h:={𝐯∈ [H^1_0(Ω)]^2   |  𝐯_|K∈𝐕^K_h,  ∀ K∈_h} and Q_h:={q∈ L^2_0(Ω)    |   q_|K∈ Q^K_h,  ∀ K∈_h}.

In the VEM framework the basis function are never explicitly computed since it would be necessary solve the PDE given in 𝐕_h^K for each element. Alternatively, they are defined by using of some polynomial projection operators (see <cit.>) and suitable degrees of freedom (dofs).

Given 𝐯∈𝐕_h^K we take the following linear operators 𝐃_𝐕, split into 

0em 
	
  * 𝐃_𝐕^1: the values of 𝐯 at the vertices of the polygon K,
	
  * 𝐃_𝐕^2: the values of 𝐯 at k-1 internal points of the (k+1)-Gauss-Lobatto quadrature rule in e∈∂ K,
	
  * 𝐃_𝐕^3: the moments of 𝐯: ∫_K𝐯·𝐠^⊥_k-2  for all 𝐠^⊥_k-2∈𝐆_k-2(K)^⊥;
	
  * 𝐃_𝐕^4: the moments of div𝐯: ∫_Kdiv 𝐯 q_k-1  for all  q_k-1∈ℙ_k-1(K)/ℝ.
 
Furthermore, for the local pressure, given q∈ Q^K_h, 
 the linear operators 𝐃_𝐐:

	
  * 𝐃_𝐐: the moments of q:
	∫_Kq p_k-1 ,	  for any  p_k-1∈ℙ_k-1(K).

We note that b(𝐯_h,q_h) for all 𝐯_h ∈𝐕_h and q_h ∈ Q_h can be computed directly from the 𝐃^1_𝐕, 𝐃^2_𝐕 and 𝐃^3_𝐕. While a(𝐯_h,𝐰_h) and (𝐟,𝐯_h) for all 𝐯_h, 𝐰_h ∈𝐕_h can
not be exactly computed. It is therefore necessary to introduce approximations a_h(𝐯_h,𝐰_h) and (𝐟_h,𝐯_h) making use of suitable polynomial projections. Further details of the construction of these bilinear forms and their related theoretical estimates can be found in <cit.>.
The discrete virtual element problem states: find (𝐮_h,p_h) s.t.

    a_h(𝐮_h,𝐯_h)+b(𝐯_h,p_h)=(𝐟_h,𝐯_h)     for all 𝐯_h ∈𝐕_h
    
    		b(𝐮_h,q_h) = 0     for all q_h ∈ Q_h
    	.




§ DOMAIN DECOMPOSITION, BDDC PRECONDITIONER, AND COARSE SPACES

We decompose 𝒯_h into N non-overlapping subdomains Ω_i with characteristic size H_i as 𝒯̅_h = ⋃_i=1^N Ω̅_i where each Ω_i is the union of different polygons of the tessellation  𝒯_h and we define Γ = ⋃_i≠ j∂Ω_i ∩∂Ω_j as interface among the subdomains. We assume that the decomposition is shape-regular in the sense of <cit.> Section 3.

We refer to the edges of the subdomains Ω_i as macro edges and we denote them with ℰ_i, moreover ℰ_ij denote the macro edge shared by the subdomains Ω^i and Ω^j.

From now we omit the underscore h since we will always refer to the finite-dimensional space and we write 𝐕× Q instead of 𝐕_h× Q_h.
We split the velocity dofs into interface (Γ) and internal (I) dofs. In particular, the 𝐃^1_V and 𝐃^2_V that belongs to a single subdomain Ω_i are classified as internal, while the ones that belong to more than a single subdomain as interface dofs. All the dofs 𝐃^3_V and 𝐃^4_V are classified as internal ones.
Following the notations introduced in <cit.> and <cit.>, we decompose the discrete velocity and pressure space 𝐕 and Q into 𝐕 = 𝐕_I ⊕𝐕_Γ and Q = Q_I⊕ Q_0, with Q_0:=∏_i=1^N {q∈Ω_i | q  is constant in Ω_i}.

𝐕_Γ is the continuous space of the traces on Γ of functions in 𝐕, 𝐕_I = ⊕_i=1^N𝐕_I^(i) and Q_I = ⊕_i=1^N Q_I^(i) are direct sums of subdomain interior velocity and pressure spaces.

We also define the space of interface velocity variables of the subdomain Ω_i by 𝐕_Γ^(i) and the associated product space by 𝐕_Γ = ∏_i=1^N𝐕_Γ^(i).
The discrete global saddle-point problem (<ref>) can be written as: find (𝐮_I,p_I,𝐮_Γ,p_0) ∈ (𝐕_I,Q_I,𝐕_Γ,Q_0) s.t.

    [
    	[    A_II  B_II^T A_Γ I^T       0;    B_II       0    B_IΓ       0;   A_Γ I  B_IΓ^T    A_ΓΓ  B_0Γ^T;       0       0  B_0Γ^T       0;         ]]
    	
    	[
    	[ 𝐮_I; p_I; 𝐮_Γ; p_0;     ]]
    	=
    	
    	[
    	[ 𝐟_I;   0; 𝐟_Γ;   0;     ]],

where the blocks · related to the continuous interface velocity are assembled from the corresponding subdomain submatrices. By static condensation one eliminates the interior variables and obtains the global interface saddle point problem

    Su = 
    	[
    	[    S_Γ B_0Γ^T;   B_0Γ      0;        ]]
    	[
    	[ 𝐮_Γ; p_0;     ]]
    	=
    	[
    	[ 𝐠_Γ;   0;     ]]
    	= 𝐠,

where S_Γ is the Schur complement of the submatrices constituted by the top 3 × 3 block of the left-hand side matrix in (<ref>) and 𝐠 the corrispective right-hand side.

We introduce 𝐕_Γ = 𝐕_Π⊕𝐕_Δ = 𝐕_Π⊕( ∏_i=1^N 𝐕_Δ^(i)), a partially assembled interface velocity space where 𝐕_Π is the continuous coarse-level primal velocity space, which dofs are shared by neighboring subdomains and the complementary space 𝐕_Δ, that is the direct sum of the subdomain dual interface velocity spaces 𝐕_Δ^(i).

In particular, the primal dofs of our problem are represented by the nodal evaluation of both components of the velocity in the vertices of the subdomain and one extra dof for each subdomain edge to satisfy the no-net-flux condition ∫_∂Ω_i𝐯_Δ^(i)·𝐧 = 0   ∀𝐯_Δ∈𝐕_Δ, where 𝐧 is the outward normal of ∂Ω_i. This last condition is needed to ensure that the operator of the preconditioned system (<ref>) with the BDDC is symmetric and positive definite on some particular subspaces <cit.>, so the preconditioned conjugate gradient (CG) method can be used for solution.

In our study we use a generalized transformation of basis approach <cit.> such that each primal basis function corresponds to an explicit dof. Firstly, on each ℰ_ij, we assume that the velocity vector should fulfill N_ij constraints c_ij^l for l = 1,...,N_ij, i.e., c_ij^l^T𝐯_Γ|ℰ^(i) = c_ij^l^T𝐯_Γ|ℰ^(j). Then, we compute the orthonormal trasformation with a modified Gram-Schmidt algorithm. Finally, since the transformations are independent of each other, we construct the resulting block diagonal global transformation. The constraints c_ij^l are established via the no-net-flux condition and the techniques to enrich the coarse space. More details  can be found in <cit.>.

For each Ω_i we introduce the scaling matrices D. These can be chosen in different ways and can be either diagonal or not, but they always must provide a partition of unity, i.e., R_D,Γ^T R_Γ = R_Γ^T R_D,Γ = I, where R_Γ,R_D,Γ :𝐕_Γ→𝐕_Γ are respectively a restriction operator, and its scaled version. We use the standard multiplicity-scaling and a variant of the deluxe-scaling to preserve the normal fluxes <cit.>.

We then define the average operator E_D = RR_D^T, which maps 𝐕_Γ× Q_0, with generally discontinuous interface velocities, to elements with continuous interface velocities in the same space. Here R and R_D^T are simply the two previous operators extended by identity to the space of piecewise constant pressures.

The preconditioner for solving the global saddle-point problem (<ref>) is then M^-1=R_D^T S^-1R_D, where S is the Schur complement system that arises using the partially assembled velocity interface functions. Theoretical estimates show that the condition number is bounded by the norm of the average operator <cit.>.


Adaptive and frugal coarse spaces

The coarse spaces is alternatively enriched by two adaptive techniques or a heuristic one. The idea is to detect the largest, or smaller, eigenvalues on each macro edge and then include the corresponing eigenvectors in the coarse space as primal constraints c_ij with the transormation of basis approach we saw before. 

We define here the frugal coarse space for our model problem. Like in the linear elasticity case, when applying the BDDC method to the Stokes problem in two dimensions, we need three constraints for each edge to control the three (linearized) rigid-body motions (two translations and one rotation). Given two subdomains Ω_l, l=i,j with diameter H_l, we have

    𝐫_1:=
    [
    [ 1    0;        ]],
      𝐫_2:=
    [
    [ 0    1;        ]],
      𝐫_3:=1/H_l[
    [ x_2 - x_2    -x_1 + x_1;                         ]],

where 𝐱∈Ω_l is the center of the rotation. Differently from the approach in <cit.>, we do not rescale the rigid body modes and we define the "approximate" eigenvector

    v(𝐱)^(m,l)_ℰ_ij:=
        
            r(𝐱)_m^(l),    𝐱∈ℰ_ij,
    
            0 ,    𝐱∈∂Ω_l ∖ℰ_ij,

for m=1,2,3 and v(𝐱)^(m) T_ℰ_ij:=[v(𝐱)^(m,i) T_ℰ_ij, -v(𝐱)^(m,j) T_ℰ_ij]. 

The three frugal edge constraints are then obtained by c_ij := B_D_ij^T S_ij P_D_ij  v(𝐱)^(m) T_ℰ_ij, where S_ij := diag(S_Γ^(i),S_Γ^(j)), B_D_ij is a scaled jump operator, B_D,ℰ_ij its restriction to the edge ℰ_ij, and P_D_ij = B_D,ℰ_ij^T B_ℰ_ij.



§ NUMERICAL RESULTS





We solve a lid-driven cavity benchmark problem on the unit square domain Ω = [0,1]×[0,1], applying Dirichlet boundary conditions on the whole ∂Ω and using a VEM implementation of degree k = 2. 
The heterogeneity is introduced to physically represent a practical example where drops (or sinkers) of a high viscosity material are spread in the fluid, in particular this is modeled defining ν(𝐱) as a continuous function that exhibits sharp gradients (Figure <ref>) <cit.>. 
These inclusions of equal size are placed randomly in the unit square domain so that they can overlap and intersect the boundary. For 𝐱∈Ω, the viscosity ν(𝐱) ∈ [ν_min,ν_max], 0<ν_min<ν_max<∞, is defined as ν(𝐱) := (ν_max-ν_min)(1-χ_n(𝐱))+ν_min. Here, χ_n(𝐱) ∈ C^∞ is an indicator function χ_n(𝐱)∈ [0,1] that accumulates n sinkers defined as χ_n(𝐱) := ∏_i=1^n 1-exp (-δmax(0,|𝐜_i-𝐱|-Ω/2)^2 ), where 𝐜_i ∈Ω, i = 1,...,n are the centers of the sinkers, Ω≥ 0 is their diameter and  δ>0 a parameter that controls the exponential decay. By choosing δ = 2000, Ω = 0.05, λ_min = 10^-3 and λ_max = 10^3, we ensure that the viscosity exibits sharp gradients. 
The right hand side is defined as 𝐟(𝐱):=(0,β(χ_n(𝐱-1))), with β = 10 to simulate gravity that takes down the high viscosity material.
In our experiments we use meshes with a Centroid Voronoi Tassellation (CVT) and Random meshes (RND), while the subdomain partitioning is performed by METIS. We compare the two adaptive coarse spaces, with TOL = 100, and the frugal one by applying the two different type of scaling mentioned before. Our numerical simulation have been performed with MATLAB R2023A  therefore no CPU time analysis is provided.

tr0.4


  

    
    < g r a p h i c s >

  



	METIS decomposition of a RND mesh. ν = 1e3 on yellow elements, ν = 1e-3 on blue ones.


In the following tables, we report the number of iterations to solve the global interface saddle-point problem (<ref>) with the PCG method, accelerated by a BDDC preconditioner, where we set the tolerance for the relative residual error to 10^-6. 

In the tables we use the following notation: nSub = number of subdomains, nSink = number of sinkers, n_Π = number of primal constraints, it = iteration count (CG), k_2 = condition number, Frugal = frugal coarse space, First =  first adaptive technique, Second = second adaptive technique. 




We consider two different tests. We first set a configuration with nSink = 11 and we increase the number of the subdomains; see Table <ref>. For both the type of the mesh considered the adaptive coarse spaces with the multiplicity scaling respect our expectations, while the frugal approach exibits a high condition number since the heuristic coarse space is not able to catch all the largest eigenvalues. Introducing the deluxe scaling we see that the number of primal constraints is drastically reduced in the adaptive coarse spaces. The frugal approach is then able to control the largest eigenvalues and performs well. 
In Table <ref> we instead keep fixed the number of subdomains at 4×4 and we increase the number of the inclusions. Again, the adaptive coarse spaces are robust and also when introducing the deluxe scaling the frugal one shows a good improvement presenting itself as a valid alternative.

spmpsci

