Domain decomposition and upscaling technique for metascreens

Abstract

Purpose

This paper aims to consider a multiscale electromagnetic wave problem for a housing with a ventilation grill. Using the standard finite element method to discretise the apertures leads to an unduly large number of unknowns. An efficient approach to simulate the multiple scales is introduced. The aim is to significantly reduce the computational costs.

Design/methodology/approach

A domain decomposition technique with upscaling is applied to cope with the different scales. The idea is to split the domain of computation into an exterior domain and multiple non-overlapping sub-domains. Each sub-domain represents a single aperture and uses the same finite element mesh. The identical mesh of the sub-domains is efficiently exploited by the hybrid discontinuous Galerkin method and a Schur complement which facilitates the transition from fine meshes in the sub-domains to a coarse mesh in the exterior domain. A coarse skeleton grid is used on the interface between the exterior domain and the individual sub-domains to avoid large dense blocks in the finite element discretisation matrix.

Findings

Applying a Schur complement to the identical discretisation of the sub-domains leads to a method that scales very well with respect to the number of apertures.

Originality/value

The error compared to the standard finite element method is negligible and the computational costs are significantly reduced.

Citation

Leumüller, M., Hollaus, K. and Schöberl, J. (2022), "Domain decomposition and upscaling technique for metascreens", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 41 No. 3, pp. 938-953. https://doi.org/10.1108/COMPEL-03-2021-0073

1. Introduction

In the context of electromagnetic compatibility, ventilation grills are also called metascreens (Holloway and Kuester, 2018). Such geometries contain regions with a locally periodic structure, for example, openings or scattering objects. This property will be called quasi-periodic. Around these objects a fine mesh is needed in finite element (FE) simulations. The goal of this work is to improve the computational times in electromagnetic wave simulations by reducing the degrees of freedom.

Many different methods have been considered in simulating shieldings with apertures. In the time domain setting, the finite difference time domain method (FDTD) is prominent and has been considered, for example, in Georgakopoulos et al. (2001) and Jiao et al. (2006). For simulations in the frequency domain, the method of moments (MoM) (Araneo and Lovat, 2009), the finite element method (FEM) (Kubík and Skála, 2015; Carpes et al., 2000; Boubekeur et al., 2014) or transmission line models (TLM) (Carpes et al., 2000; Nie et al., 2011) are common. The cited works above mainly focus on shieldings with a single or few apertures. The present work deals with simulations for shieldings with many apertures. In Ali et al. (2005) and Li et al. (2000) the FDTD and the MoM are considered for such examples. Another approach for large numbers of quasi-periodic apertures is homogenisation which has been studied in Bardi et al. (2006) and Cruciani et al. (2015). In this work, the idea of the Nitsche-type mortaring finite element method for the vector potential wave equation considering non-matching meshes introduced in Hollaus et al. (2010) and Heinrich and Nicaise (2001) is altered to fit the quasi-periodic setting with many apertures.

The method will be used to establish an upscaling by non-matching meshes. In (Heinrich and Nicaise, 2001) the structure of the resulting system matrix is exploited to reduce the computational effort of solving the system of linear equations by applying a Schur Complement (SC). This work applies the method to the electromagnetic wave problem with multiple sub-domains. The main contribution of this work is the efficient implementation of a SC for quasi-periodic structures. This is achieved by discretising each aperture identically.

The weak formulation for the vector potential wave equation in the time harmonic setting is derived in Section 2. The quasi-periodic structure is exploited via a domain decomposition (DD) approach where many of the generated sub-domains are identical. This property is further used by introducing non-matching meshes, leading to an upscaling (Graham and Scheichl, 2007) and applying the Nitsche-type mortar finite element method (NMFEM) (Hollaus et al., 2010) to cope with the non-matching mesh in Section 3. In Section 4 the resulting system matrix for the NMFEM and its special structure are highlighted. To exploit this structure a SC is applied, to effectively eliminate all the sub-domains from the computation reducing the computation time, in Section 5. The proposed method is applied to various electromagnetic wave problems for 2 D geometries as well as 3 D geometries, in Section 6. The computational times for the NMFEM and the FEM are compared in Section 6.

2. Domain decomposition with upscaling for the time harmonic wave equation

The vector potential wave equation in the time harmonic setting using a magnetic vector potential A is given as:

(1) curlμ−1curl A−κ2A=0 in Ω,

(2) A×n=0 on ΓD,

(3) μ−1curl A×n=0 on ΓN,

(4) μ−1curl A×n+jωZ(n×A)×n=0 on Γt,

(5) μ−1curl A×n+jωZ(n×A)×n=h×n on Γe

with κ2=jωσ+ω2ε and Ω=∪i=0nΩi⊂ℝ3. The physical parameters μ, ε and σ are the magnetic permeability, the electric permittivity and the conductivity, respectively, ω = 2πf is the angular frequency, j is the imaginary unit, h is a prescribed boundary field and Z=με is the wave impedance. Sub-domains are defined in Figure 1. To derive the, for FEM needed, weak formulation the following function spaces are needed:

(6) L2(Ω):={v:Ω→ℂ;∫Ω|v|2 dx<∞},

(7) H(curl,Ω):={v∈[L2(Ω)]3;curlv∈[L2(Ω)]3},

(8) H0(curl,Ω):={v∈H(curl,Ω);v×n|ΓD=0},

(9) L2(Γ):={v:Γ→ℂ;∫Γ|v|2 ds<∞}.

Equation (1) is multiplied by v∈H0(curl,Ω) and then integration by parts is carried out. This leads to the weak formulation: Find A∈H0(curl,Ω) so that:

(10) ∫Ωμ−1curl A·curlv−κ2A·v dx−jωZ∫Γe∪ΓtA×n·v×n ds=∫Γeh×n·v×n ds

for all v∈H0(curl,Ω).

To be able to exploit the quasi-periodic structure the whole domain is decomposed into the exterior domain Ω₀ and the sub-domains Ωi,i∈{1,…,n}, see Figure 1. The domains are non-overlapping and share the interface Γ. This interface is not only between Ω₀ and the Ω_i, but also between the sub-domains. The Ω_i are chosen such that each represents exactly one aperture and Ω₀ encloses all Ω_i. In the decomposition the sub-domains are geometrically identical; therefore, each can be discretised in the same way. This means that, independent of the number of apertures, only one Ω_i has to be meshed and the FE system matrix for just a single sub-domain needs to be assembled.

The DD can be used to address a second issue. The small geometric structure around an aperture leads to a very fine discretisation. This entails a fine discretisation in the vicinity of the aperture due to the necessarily regular, conforming mesh for the FEM. A sufficiently accurate solution could be represented by a coarser mesh; therefore, the fine mesh unnecessarily increases the computational costs. The goal is to discretise the exterior domain as coarse as possible such that the wave propagation is appropriately resolved and at the same time discretise the sub-domains fine enough to grasp the local, small-scale behaviour of the solution around the apertures. This leads to an upscaling between the solutions on the Ω_i and the solution on Ω₀. The resulting non-conforming mesh over Γ can be seen in Figure 2 and 3. The FEM is not capable to cope with non-conforming meshes. To address this issue the NMFEM is applied.

3. Nitsche-type mortar finite element method

The general idea of the NMFEM is to break the strong continuity of the FE solution on domain boundaries and reinforce it in a weak sense. This is achieved by introducing a FE space X(Γ)=[L2(Γ)]3 which is defined on Γ, see Figure 3. Using the NMFEM, the FEs on different sub-domains are not coupled directly, but only via the interface FEs on Γ. In an usual discontinuous Galerkin method the continuity of the numerical flux is enforced by testing with test functions from both side of the interface, which leads to a coupling between elements. This is avoided by the NMFEM, because the continuity of the normal flux across interfaces is tested by new interface variables. Therefore FEs on different sub-domains do not couple directly (Hollaus et al., 2010; Cockburn et al., 2009).

3.1 NMFEM for the time harmonic wave equation

The vector potential wave equation in the time harmonic setting using local magnetic vector potentials Ai for each sub-domain is given as:

(11) curlμ−1curlAi−κ2Ai=0 in Ωi,

(12) Ai×ni=−Aj×nj=:A^×ni on Γ,

(13) μ−1curlAi×ni+μ−1curlAj×nj=0 on Γ,

(14) Ai×ni=0 on ΓD,

(15) μ−1curlA0×n0=0 on ΓN,

(16) μ−1curlA0×n0+jωZ(n0×A0)×n0=0 on Γt,

(17) μ−1curlA0×n0+jωZ(n0×A0)×n0=h×n0 on Γe.

The vector ni is the outward pointing normal vector on the boundary of Ω_i. The tangential continuity of the Ai over Γ is enforced by equation (12). The vector field A^ represents the tangential component of the Ai and is needed to derive the weak formulation for the NMFEM. The tangential continuity of the flux is ensured by equation (13).

To obtain the weak formulation for the NMFEM, equation (11) is multiplied with a test function vi∈H0(curl,Ωi) on each domain and then integration by parts is carried out. This leads to the initial weak formulation:

(18) ∑i∫Ωiμ−1curlAi·curlvi−κ2Ai·vi dx−∫Γ(μ−1curlAi×ni)·vi ds−jωZ∫Γe∪ΓtA0×n0·v0×n0 ds=∫Γeh×n0·v0×n0 ds.

This closely resembles equation (10) with the difference that the volume term is separated into the sub-domains and due to the fact that the solution is not continuous over Γ, the boundary integrals appear. This weak formulation is incomplete in the sense that the continuity over Γ is not yet considered. To be able to enforce the field trace continuity in equation (12) and equation (13), a vector valued space X(Γ) of traces of Nédélec kind [Schöberl and Zaglmayr (2005)] are needed. Multiplying equation (13) with v^∈X(Γ) and integrating over Γ leads to the term:

(19) 0=∫Γ(μ−1curlAi×ni+μ−1curlAj×nj)·v^ ds,

which is zero if the flux of the solution is continuous. Adding this term to equation (18) yield:

(20) ∑i∫Ωiμ−1curlAi·curlvi−κ2Ai·vi dx−∫Γ(μ−1curlAi×ni)·(vi−v^) ds−jωZ∫Γe∪ΓtA0×n0·v0×n0 ds=∫Γeh×n0·v0×n0 ds.

Using:

(21) (μ−1curlAi×ni)·(vi−v^)=−μ−1curlAi·(vi−v^)×n

for the boundary term and adding:

(22) ∑i∫Γμ−1curlvi·((Ai−A^)×ni) ds=0

to symmetrise and:

(23) ∑i∫Γαp2μh((Ai−A^)×ni)·((vi−v^)×ni) ds=0

to stabilise the formulation, which are zero because of equation (12), leads to:

Find Ai∈H0(curl,Ωi),i∈{1,…,n},A^∈X(Γ) so that:

(24) ∑i∫Ωiμ−1curlAi·curlvi−κ2Ai·vi dx+∫Γμ−1curlAi·((vi−v^)×ni) ds+∫Γμ−1curlvi·((Ai−A^)×ni) ds+∫Γαp2μh((Ai−A^)×ni)·((vi−v^)×ni) ds−jωZ∫Γe∪ΓtA0×n0·v0×n0 ds=∫Γeh×n0·v0×n0 ds.

for all vi∈H0(curl,Ωi),i∈{1,…,n},v^∈X(Γ). The scalar α > 0 represents a sufficiently large stabilisation factor of the weak formulation which is independent of the polynomial degree p of the finite element space and the mesh size h.

In simulations with a geometry where one dimension is considered as infinitely long the original vector valued problem is reduced to the scalar Helmholtz problem defined on a 2 D geometry. The steps for deriving the weak formulation for the Helmholtz equation are similar to the steps for the vector potential wave equation. The biggest difference is that instead of the vector valued space X(Γ) the scalar valued space L²(Γ) is needed to reinforce the field trace continuity over Γ.

3.2 Interface finite element spaces

Finite elements for the space H0(curl,Ωi) as well as for the interface spaces L²(Γ) and X(Γ) are needed. The first space is a FEM space described in (Schöberl and Zaglmayr, 2005; Zaglmayr, 2006). For the space L²(Γ) defined on lines in 2 D, Legendre polynomials of a fixed degree are used. To discretise X(Γ) on rectangles in 3 D, tensor products of Legendre polynomials are considered. The polynomial degree for L²(Γ) and X(Γ) can be chosen independently of the polynomial degree of H0(curl,Ωi).

4. System matrix for NMFEM

The system of linear equations for the NMFEM:

(25) Mx=f

has the following block structure:

(26) M:=(M00M0Γ00⋯0MΓ0MΓΓMΓ1MΓ2⋯MΓn0M1ΓM110⋯00M2Γ0M22⋱⋮⋮⋮⋮⋱⋱00MnΓ0…0Mnn),

(27) x:=(x0,xΓ,x1,x2,…,xn)⊤,

(28) f:=(f0,fΓ,f1,f2,…,fn)⊤.

The right hand side vector can be split into f₀ corresponding to the excitation on Ω0, fΓ to Γ and f_i to the sub-domain Ω_i. In a similar way the unknowns x are split. Sub-blocks with mixed indices in equation (26) indicate the coupling matrices between different FE spaces. For example, the matrix M_0Γ represents the coupling between the FE spaces defined on Ω₀ and Γ. Due to the identical discretisation of Ω_i, the matrices M_ii are identical for i∈{1,…,n}. As can be seen in equation (24) the weak formulation is symmetric; therefore, the system matrix is also symmetric and MiΓ=MΓi. These properties, as well as the structure of the system matrix, will be exploited so that the whole system matrix does not need to be assembled.

5. Schur complement

To efficiently use the introduced technique, a SC is applied to solve the resulting system of linear equations. In the SC, the inverse of the system matrix of a single sub-domain can be used for all sub-domains, reducing the computational costs. Due to the small size of the system matrix of a single sub-domain, the inverse can be cheaply calculated. The sparsity of the system matrix is retained, because the FEs on Ω_i only couple with a few degrees of freedom (DoF) on Γ.

The SC of a matrix M is defined as:

(29) M=(ABCD)=(IBD−10I)(A−BD−1C00D)(I0D−1CI)

and therefore the inverse of M can be written as:

(30) M−1=(I0−D−1CI)((A−BD−1C)−100D−1)(I−BD−10I),

if D is invertible. The matrix A – BD^{– 1}C is always invertible if the matrices M and D are invertible. For M in equation (26) the blocks are:

(31) A=(M00M0ΓMΓ0MΓΓ),

(32) B=(00⋯0MΓ1MΓ2⋯MΓn),

(33) C=(0M1Γ0M2Γ⋮⋮0MnΓ),

(34) D=(M110⋯00M22⋱⋮⋮⋱⋱00…0Mnn).

Using this approach the system of linear equations can be solved in the three steps given below for the case of two sub-domains.

5.1 Pre-processing

The first step has the form:

(35) (I−BD−10I)f=(I0000I−MΓ1(M11)−1−MΓ2(M22)−100I0000I)(f0fΓf1f2)=(f0f˜Γf1f2)⊤=:f˜.

As can be seen, only f_Γ is altered. To assemble the necessary matrices MΓi(Mii)−1 systems of linear equations with M_ii have to be solved. The total number of needed solutions is equal to the number of degrees of freedom (N_DOF) of the boundary space on Γ ∩ ∂Ω_i. This might be time expensive depending on the size of the sub-domain discretisation matrix and on the N_DOF on the boundary, but due to the identical discretisation of each Ω_i and of the boundary spaces around the sub-domains, the matrices MΓi(Mii)−1 are, up to permutations, identical for all sub-domains and therefore only one needs to be assembled.

5.2 Solving

The second step is given by:

(36) ((A−BD−1C)−100D−1)f˜=(S−1000(M11)−1000(M22)−1)(f0f˜Γf1f2)=(x0xΓx˜1x˜2)=:x˜,

Where:

(37) S:=(M00M0ΓMΓ0MΓΓ−∑i=12MΓi(Mii)−1MiΓ).

Instead of one large system of linear equations, smaller systems of linear equations on the individual sub-domains can be considered in parallel. On the sub-domains solutions with M_ii have to be calculated. On the exterior domain a system of linear equations with S has to be solved. To assemble the matrices MΓi(Mii)−1MiΓ in S the matrix MΓi(Mii)−1, which has already been assembled in the first step, can be used.

5.3 Post-processing

The third step consists of:

(38) (I0−D−1CI)x^=(I0000I000−(M11)−1M1ΓI00−(M22)−1M2Γ0I)(x0xΓx˜1x˜2)=(x0xΓx1x2).

In this step only the solution vectors of the sub-domains are altered. Due to the symmetry of the system matrix, the matrix assembled in the first step can be used for the third step.

In case of no excitations on the sub-domains, fi=0,i∈{1,…,n}, the first step can be omitted. If the solution on the sub-domains is not needed, the third step can be neglected.

6. Numerical example

To show the feasibility of the introduced method, 2 D and 3 D ventilation grills have been considered.

6.1 2 D Ventilation grill with one aperture

As a 2 D example, a ventilation grill with one aperture, see Figure 4, has been chosen. The shielding Γ_D has been modelled by homogeneous Dirichlet boundary conditions. Homogeneous Neumann boundary conditions have been applied on Γ_N. The bottom boundary Γ_t has been considered as transparent and modelled by first order absorbing boundary conditions (ABC) (Chatterjee et al., 1993). On the boundary Γ_e at the top a plane wave and a first order ABC have been applied via a Robin boundary condition with h=(jωZ,0)⊤. The permeability and permittivity of vacuum have been used and a frequency f = 10 GHz is considered. The reference solution with FEM and the solution with NMFEM have been compared. The logarithmically scaled norm of the solutions can be seen in Figure 4. The relative L² error on Ω₀ of the NMFEM solution is:

(39) ||Aref−ANMFEM||L2(Ω0)||Aref||L2(Ω0)·100%=0.06%.

6.2 3 D Ventilation grill with 25 apertures

As a more realistic 3 D example, a ventilation grill with 5 × 5 apertures, see Figure 5 and 6, has been considered.

The boundary value problem equation (24) has been solved at the frequency f = 10 GHz. For the permeability μ and the permittivity ε vacuum has been chosen. The excitation has been set to:

(40) h=(exp⁡(−(x−20)21002),0,0)⊤.

Although, for instance, perfectly matched layers (Leumüller et al., 2019) would be more accurate, for the sake of simplicity first-order absorbing conditions have been prescribed on the outer boundary Γt∪Γe. The shielding Γ_D has been modelled by homogeneous Dirichlet boundary conditions. Finally, this work focuses on upscaling and DD. The solution with the NMFEM on the exterior domain can be seen in Figure 7. All the numerical examples have been calculated with the aid of the open source FE software Netgen/NGSolve (Schöberl, 2022).

The computational times for the FEM simulation and the NMFEM simulation have been compared. A constant maximum size for the mesh elements has been chosen in all simulations. The meshes around the apertures are still very fine in the simulations. For the H(curl) FE spaces the polynomial degree has been set to p = 3 in all simulations. The polynomial degree of the boundary FE space on Γ has been set to p_Γ = 6. All calculations have been carried out on 16 cores in parallel. As a solver for the systems of linear equations, PARDISO (De Coninck et al., 2016; Verbosio et al., 2017; Kourounis et al., 2018) has been used. The number of sub-domains N_SD has been varied between N_SD ∈ {1, 4, 9, 16, 25}. A list of the N_DOF for the different simulations can be seen in Table 1. In Figure 8 the simulation times for NMFEM and FEM are shown. For small N_SD the computation time for the sub-domain matrices in the SC, e.g. MΓi(Mii)−1, is dominant, but for a more realistic number of apertures the factorisation and solving with PARDISO dominate. In the FEM simulation the computation time is always dominated by the factorisation and solving with PARDISO. The meshing and assembling times for the system matrices are small in comparison to the factorisation times. The maximal used memory for NMFEM and FEM with respect to N_SD can be seen in Figure 9. Using NMFEM the maximal memory demands are enormously reduced.

7. Conclusion

A NMFEM with non-matching meshes for 2 D and 3 D ventilation grills with quasi-periodic apertures has been introduced. The method leads to a large reduction of the computational time and the memory demands through upscaling and SC reduction. This has been achieved by an efficient implementation of the SC. The identical discretisation of each sub-domain has been essential to decrease the N_DOF in the simulations.

Figures

Figure 1.

Domain decomposition of the domain (left) into the exterior domain Ω₀ and the sub-domains Ω_i (right)

Figure 2.

Coarse mesh on the exterior domain (left), fine mesh on the sub-domain (right)

Figure 3.

A detail of the interface Γ between Ω₀ and Ω_i on which the interface finite element space is defined

Figure 4.

Reflection and transmission of an incident plane wave from above of the aperture, reference solution with fine mesh (left), domain decomposition solution with upscaling i.e. fine and coarse mesh (right)

Figure 5.

Geometry of the exterior domain of a ventilation grill

Figure 6.

Geometry of the periodic sub-domains with quadratic openings; boundaries between the sub-domains are indicated by dashed lines, front view (above) and top view (below)

Figure 7.

Incident, reflected and transmitted field of a plane wave impinging a shielding with an array of 5 × 5 apertures

Figure 8.

Computation time for solving with the NMFEM and the FEM with respect to the number of sub-domains

Figure 9.

Maximal memory requirement in gigabyte (GB) for solving with the NMFEM and the FEM with respect to the number of sub-domains

Table 1.

The NDOF for NMFEM and FEM with respect to the number of sub-domains

NSD	1	4	9	16	25
NMFEM
NDOF,Ωe	21,908	40,384	80,348	126,092	164,432
NDOF,Γ	432	1,440	3,024	5,184	7,920
NDOF,Ωi	19,060	19,060	19,060	19,060	19,060
NDOF,Γ∩∂Ωi	432	432	432	432	432
FEM NDOF	66,000	187,824	408,048	707,364	1,073,208