Goal-oriented error estimation and h-adaptivity for Maxwell’s equations

Goal-oriented error estimation and h-adaptivity for Maxwell’s equations

Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616 www.elsevier.com/locate/cma Goal-oriented error estimation and h-adaptivity for MaxwellÕs equ...

735KB Sizes 1 Downloads 94 Views

Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616 www.elsevier.com/locate/cma

Goal-oriented error estimation and h-adaptivity for MaxwellÕs equations P. Ingelstr€ om *, A. Bondeson Department of Electromagnetics, Chalmers University of Technology, SE-412 96 G€oteborg, Sweden Received 13 November 2002; received in revised form 28 February 2003

Abstract Goal-oriented error estimates, based on dual solutions, are derived for the S-parameters of a waveguide cavity resonator and a computationally cheap method is proposed to compute these. Numerical results show that the error estimators are relatively accurate, usually well within a factor 2 from the exact errors. Adaptive mesh-refinement based on these estimates recovers optimal convergence rates for complete and incomplete, first and second order, curl-conforming finite elements of Nedelec type, despite singularities at reentrant corners. It is also shown numerically that the convergence rate is independent of the element-order when singularities are present and the mesh is refined uniformly. Ó 2003 Elsevier Science B.V. All rights reserved. Keywords: MaxwellÕs equations; Adaptivity; Goal-oriented error estimation; Hierarchic basis functions; Waveguide simulation

1. Introduction In this article we derive, apply and evaluate the performance of goal-oriented error estimation [1], to a finite element scheme for solving the two-dimensional vector wave equation derived from MaxwellÕs equations. To avoid spurious solutions that occur with node-based vector elements, we use the tangentially continuous (curl-conforming) basis functions of Nedelec-type [2,3]. To our knowledge, this is the first work on goal-oriented error estimation with finite elements of Nedelec-type. A common goal for computations is to give as accurate results as needed with as little work as possible. Today, several first order methods are very popular, e.g. the Finite-Difference Time-Domain (FDTD) method [4] and the Finite Element (FE) method [5]. For many smaller problems, where an accuracy of about 1% is enough, these methods work well. However, most first order methods suffer from considerable dispersion errors that make them less attractive for large problems, i.e. problems where the size of the objects are many wavelengths. They also converge rather slowly to the correct solution, as h2 / N 2=d if no singularities are present, where h is the typical cell size, N is the number of degrees of freedom and d is the dimension of the problem, in this case d ¼ 2. *

Corresponding author. E-mail address: [email protected] (P. Ingelstr€ om).

0045-7825/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0045-7825(03)00295-0

2598

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

Higher-order methods, e.g. the FE method which we use in this article, are often better alternatives when high accuracy is needed. Under good conditions they converges as N 2p=d , where p is the order of the method. The problem is that most real electromagnetic applications involve non-smooth, non-convex geometries (i.e. with reentrant corners) where the fields will be singular. In these cases, if meshes with uniform cell sizes are used, the convergence rate will be limited to ha , for some a < 2 depending on the geometry of the singularities, but not on the order p. Hence there is little use for higher-order methods unless we combine them with adaptive mesh refinement––h-adaptivity. We note that while higher-order methods are more severely affected by singularities, also first order methods will benefit from h-adaptivity. Another attractive strategy is to use hp-adaptivity, where both the cell size h and the order of the method p can vary. Work on hp-methods in electromagnetics can be found in [6–8]. If the hp-adaptivity works correctly, exponential convergence rates should be obtained. However, it seems not yet clear how the element sizes and orders should be determined to obtain optimal convergence for the curl-conforming elements used in electromagnetics. A general strategy for the construction of adaptive meshes is to use local a posteriori error estimators and refine the mesh where these are large [9]. Despite a substantial amount of theoretical work on a posteriori error estimators for MaxwellÕs equations, [10–13] computational studies in this area have only begun. It is common practice to estimate the error in global energy (or L2 ) norm. Using goal-oriented error estimators we instead estimate, and compute, the error of some functional of the solution. A survey of goaloriented error estimators can be found in [14]. There are mainly two advantages of goal-oriented error estimators. Firstly, we are usually more interested in knowing some functional depending on the solution rather than the solution itself. For example, the work presented here aims at finding the S-parameters of a microwave component. The detailed behaviour of the field in every point in space is not of great interest. Secondly, when refining the mesh, we obtain a mesh that is ‘‘optimal’’ for computing the desired goalfunctional instead of a mesh that minimizes the energy-error everywhere. The paper is organized as follows. In Section 2 we introduce our model problem: To compute the Sparameters for a two-dimensional waveguide cavity resonator. The S-parameters are used to characterize microwave components of various kinds, and to compute these is a common and important task in microwave applications. In Section 2 we also derive error estimates for the S-parameters based on solutions of dual problems and give several alternatives for how to obtain approximate dual solutions with a small amount of additional work. In Section 3, we give results from convergence studies of the S-parameters and compare the error estimates obtained using the different dual approximations and different goal-functionals. In these studies, we use incomplete [2] and complete order [3] curl-conforming triangular finite elements of orders 1, 2 and 3 for uniform refinement and of orders 1 and 2 for adaptive refinement.

2. Problem description The model problem we have chosen is to compute S11 (reflection coefficient of the lowest waveguide mode) and S12 (transmission coefficient for the lowest waveguide mode) for an air-filled, waveguide cavity resonator with metal walls. We also compute error estimates, e11 and e12 , for S11 and S12 respectively. The resonator contains two irises with eight reentrant corners where the field will be singular. This makes it a suitable test case for an algorithm with adaptive h-refinement. 2.1. Geometry The geometry of the waveguide in the xy-plane is shown in Fig. 1 and the dimensions are given in Table 1. We denote the unspecified height of the waveguide in the z-direction by hz . In the rest of this article we will let X symbolize the computational domain, which is enclosed by the metal walls ðC0 Þ, Port 1 ðC1 Þ and Port 2

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2599

Fig. 1. Geometry of the waveguide cavity resonator. The shaded region correspond to metal walls.

Table 1 Dimensions of the waveguide cavity resonator shown in Fig. 1 Symbol

Length (mm)

Description

d l w t s

36 24 39 3 21

Distance from the irises to the ports Length of the resonating cavity Width of the waveguide Thickness of the iris walls Width of the iris slots

ðC2 Þ. On the boundary we define a normal unit vector n^ pointing outwards from X and a tangential unit vector ^t which has X on its left-hand side so that ^t ¼ ^z  n^. 2.2. The time harmonic vector wave equation For computations in the frequency domain, we assume that the electric and magnetic fields in the timedomain are time-harmonic with the angular frequency x and can be represented by corresponding complex fields in the frequency domain, E ¼ Eðx; y; z; xÞ and H ¼ Hðx; y; z; xÞ, E TD ðx; y; z; tÞ ¼ RefEðx; y; z; xÞejxt g; H TD ðx; y; z; tÞ ¼ RefHðx; y; z; xÞejxt g; where ÔjÕ is the imaginary unit. The frequency domain amplitudes satisfy the Maxwell equations r  H ¼ þjxE; r  E ¼ jxlH; where the material parameters  and l are the electric permittivity and the magnetic permeability respectively. Combining the two equations above, we eliminate E and obtain the vector wave equation, 1 r  r  H  x2 lH ¼ 0: 

ð1Þ

The geometry of the waveguide is independent of z (for 0 < z < hz ). We will assume that also the excitation is independent of z so that the magnetic field varies only in two dimensions, H ¼ Hðx; y; xÞ. Further, due to

2600

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

boundary conditions, the z-component of the magnetic field, Hz must be zero, hence H ¼ Hx ðx; y; xÞ^ xþ Hy ðx; y; xÞ^ y . Similarly E / r  H can only have a z-component, E ¼ Ez ðx; y; xÞ^z. 2.3. Boundary conditions Inside an air-filled waveguide, which is uniform in the longitudinal ðxÞ direction, the magnetic field can be expressed as a sum of modal waves with different amplitudes, see e.g. [15]. For the waveguide considered in this paper, the magnetic field for a modal wave with unit amplitude propagating in the positive x-direction is given from cm x Hðx; yÞ ¼ H þ ; m ðyÞe

m ¼ 1; 2; . . . ;

where H m is a vector phasor which we will refer to as the ‘‘modal field’’ for the mth mode, cm is the corresponding propagation constant and the coordinates x and y are defined as in Fig. 1. Similarly a modal wave propagating in the negative x-direction is given from þcm x Hðx; yÞ ¼ H  ; m ðyÞe

m ¼ 1; 2; . . .

The modal fields are given from  mpy   mpy  mp ^ lH ¼ sin cos y þ x^; m w wcm w where w is the width of the waveguide, and the propagation constants are rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  mp 2  x2 l: cm ¼ w When deriving the boundary conditions for the ports, we will assume that the waveguide is operated at a frequency f ¼ 2px such that cm is imaginary for the first mode ðm ¼ 1Þ and real for all other modes. Hence only the first mode will propagate and all other modes will decay exponentially. To stress this we set c1 ¼ jb from now on. We also assume that the ports are sufficiently far from any discontinuities in the waveguide so that the amplitudes of the higher modes at the ports are negligible. The magnetic field close to the ports can then be expressed as a sum of two first order modal waves jbx þjbx H Aþ H þ þ A H  ; 1e 1e

where the coefficients Aþ and A are the amplitudes of the waves propagating in the positive and negative xdirections respectively. To obtain boundary conditions suitable for implementation, we use combinations of r  H and H  n^ to eliminate the normal components and the unknown amplitudes. Port 1. Here n^ ¼ ^ x, ^t ¼ ^ y , and we have the following relations:  1 py  þ jbx x2 n^  r  H ¼  sin þ A eþjbx Þ ^t; ðA e  w jb  py  n^  lH  n^ ¼  sin ðAþ ejbx  A eþjbx Þ^t: w Multiplying the second equation with x2 =jb and then adding the two equations we obtain the boundary condition for Port 1,  py  1 x2 x2 ^t; sin n^  r  H þ n^  lH  n^ ¼ 2Aþ 1  w jb jb þ jbx jPort 1 . where Aþ 1 ¼ A e

ð2Þ

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

Port 2. Here n^ ¼ x^, ^t ¼ y^, and in the same way we obtain the boundary condition  py  1 x2 x2 ^t; sin n^  r  H þ n^  lH  n^ ¼ 2A 2  w jb jb

2601

ð3Þ

 þjbx where A jPort 2 . 2 ¼ A e

Metal walls. At the metal walls of the waveguide, which we assume to be perfect electric conductors, we have the homogeneous Neumann boundary condition, 1 n^  r  H ¼ 0; 

ð4Þ

which states that the tangential component of E, which is proportional to r  H, is zero on the boundary. 2.4. The original problem and its weak formulation The partial differential equation is the vector wave equation in Eq. (1) together with the boundary conditions defined by Eqs. (2)–(4). The waveguide is fed with a wave of unit amplitude through Port 1,  which implies that Aþ 1 ¼ 1 and A2 ¼ 0. 1 r  r  H  x2 lH ¼ 0 

ð5Þ

in X;

 py  1 x2 2x2 ^t on C1 ; sin n^  r  H þ n^  lH  n^ ¼   w jb jb

ð6Þ

1 x2 n^  r  H þ n^  lH  n^ ¼ 0  jb

ð7Þ

1 n^  r  H ¼ 0 

on C2 ;

ð8Þ

on C0 :

In order to derive a weak form of Eq. (5) we define the symmetric inner products ZZ Z u  v dx dy; hu; viC ¼ uv dl ðu; vÞX ¼ X

ð9Þ

C

and the space 2

2

Hðcurl; XÞ ¼ fu : ðu; u ÞX þ ðr  u; r  u ÞX < 1g;

ð10Þ

where u denotes the complex conjugate of u. The derivation follows the standard procedure (see e.g. [5]): (1) Multiply Eq. (5) with a test function V 2 Hðcurl; XÞ. (2) Use the vector identity u  r  v ¼ v  r  u  r  ðu  vÞ, with u ¼ V and v ¼ 1 r  H to move one of the curl-operators from H to V. (3) Integrate the equation over the domain X. (4) Apply GaussÕ divergence theorem to change the divergence term into a boundary integral. (5) Finally, use the boundary conditions to change the boundary integral into a convenient form. The obtained weak formulation reads: Find H 2 Hðcurl; XÞ such that aðH; VÞ ¼ b1 ðVÞ 8V 2 Hðcurl; XÞ;

ð11Þ

2602

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

where the operators aðH; VÞ and b1 ðVÞ are defined by   1 x2 r  H; r  V aðH; VÞ ¼  x2 ðlH; VÞX  hlH  ^t; V  ^tiC1 þC2  jb X

ð12Þ

and bn ðVÞ ¼

E 2x2 D  py  sin ; V  ^t ; w jb Cn

n ¼ 1; 2:

ð13Þ

The functional b1 ðÞ corresponds to the source, or load, of the problem, and the subindex indicates from which port the waveguide is fed. The reason for this notation is that the dual problems that occur when the goal-functional is S11 or S12 can be fed either through Port 1 or Port 2. From now on, whenever the subindex n occurs it can take either the value 1 or 2. Note that we could equally well use Hermitian inner products, where the second argument is complex conjugated, instead of the symmetric definition (9). The advantage of using the symmetric formulation is only formal; it removes the complex conjugates in the dual problems. 2.5. The finite element formulation We use curl-conforming finite elements of Nedelec-type similar to those developed by Webb in [16]. The only difference is that we have modified the basis functions associated with edges so that their tangential components along the edges are Legendre polynomials. This makes the mass matrix slightly more sparse and makes it easier to implement the boundary conditions. The advantage with WebbÕs basis functions is that they are hierarchical, which is essential for the algorithm we will use later to compute the dual solutions. In the present work we use a direct solver and have not considered the conditioning of the resulting matrices. For iterative solvers, other sets of basis functions might be preferable. Other curl-conforming hierarchical basis functions on triangles can be found in e.g. [8]. Let N ph denote either the space of incomplete [2] or complete order [3] curl-conforming finite elements of order p on a triangulation of X with typical element size h. The dimension Nhp of the space N ph is shown in Table 2 and depends on the order p and the number of edges Ned and elements Nel in the triangular mesh. The finite element formulation corresponding to Eq. (11) now reads: Find H ph 2 N ph such that aðH ph ; VÞ ¼ b1 ðVÞ

8V 2 N ph ;

ð14Þ

H ph

is the discrete solution of order p that approximates the true solution H. Expanding where functions V i ,

H ph

in basis

p

H ph

¼

Nh X

ðhph Þi V i ;

ð15Þ

i¼1

where hph is a vector of expansion coefficients, Eq. (14) gives a linear system of equations, Aph hph ¼ bph;1 ;

ð16Þ

Table 2 Dimension of the space N ph p

Nhp (incomplete order)

Nhp (complete order)

1 2 3

Ned 1:5Nel 2Ned þ 2Nel 5:0Nel 3Ned þ 6Nel 10:5Nel

2:0Ned 3Nel 3Ned þ 3Nel 7:5Nel 4Ned þ 8Nel 14:0Nel

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2603

where the elements in the matrix Aph and the vector bph;1 are defined from ðAph Þij ¼ aðV i ; V j Þ; ðbph;n Þj ¼ bn ðV j Þ;

i; j ¼ 1; 2; . . . ; Nhp ; j ¼ 1; 2; . . . ; Nhp :

The second subscript of bph;n labels the feeding port. 2.6. Definition of S11 and S12 S11 is defined as the ratio between the amplitudes of an inserted modal wave at Port 1 and its reflected wave at Port 1. Similarly S12 is the ratio between the inserted modal wave at Port 1 and its transmitted wave at Port 2. We will now derive a mathematical expression for these. Remembering that the fields at the ports can be described by two first order modal waves, we write the magnetic field at the two ports as þ þ     HjPort 1 ¼ Aþ 1 H 1 þ A1 H 1 ¼ H 1 þ A1 H 1 ; þ   þ þ HjPort 2 ¼ Aþ 2 H 1 þ A2 H 1 ¼ A 2 H 1 ; þ   where we have used the boundary conditions Aþ 1 ¼ 1 and A2 ¼ 0. The unknown amplitudes A1 and A2 can be computed by taking the scalar product of sinðpy=wÞ and the tangential component of the equations above,

hsinðpy=wÞ; H  ^tiC1 ¼ ð1  A 1 Þhsinðpy=wÞ; sinðpy=wÞiC1 ; hsinðpy=wÞ; H  ^tiC2 ¼ Aþ 2 hsinðpy=wÞ; sinðpy=wÞiC2 : Now using the definition of bn ðÞ in Eq. (13), the S-parameters can be written as A 1 ¼ A 1 ¼ 1  Cb1 ðHÞ; Aþ 1 Aþ S12 ¼ 2þ ¼ Aþ 2 ¼ Cb2 ðHÞ; A1

S11 ¼

ð17Þ

where C ¼ jb=wx2 . We define the error e of the finite element solution as e ¼ H  H ph and write the errors in S11 and S12 as e11 ¼ ð1  Cb1 ðHÞÞ  ð1  Cb1 ðH ph ÞÞ ¼ Cb1 ðeÞ; e12 ¼ Cb2 ðHÞ  Cb2 ðH ph Þ ¼ Cb2 ðeÞ:

ð18Þ

2.7. Error estimates involving dual solutions We use a dual, or adjoint, problem to derive the desired error estimates. Because MaxwellÕs equations are self-adjoint, the dual problem is the same as the original except for the boundary conditions. The dual problem is fed through Port 1 if we want an error estimate for S11 and fed through Port 2 if we want an error estimate for S12 . The weak form of the dual problem reads: Find G n 2 Hðcurl; XÞ such that aðG n ; WÞ ¼ bn ðWÞ

8W 2 Hðcurl; XÞ;

ð19Þ

where G n is the solution to the dual problem fed through port n. Note that G 1 is identical to H. This is a manifestation of reciprocity of backscattering for the self-adjoint MaxwellÕs equations, see e.g. [17].

2604

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

In order to derive the error estimates, we let W take the value of the error of the finite element solution, W ¼ e ¼ H  H ph , and then use the symmetry of the operator að; Þ to get aðH; G n Þ  aðH ph ; G n Þ ¼ bn ðeÞ:

ð20Þ

Since G n 2 Hðcurl; XÞ we can use Eq. (11) with V ¼ G n to eliminate H, b1 ðG n Þ  aðH ph ; G n Þ ¼ bn ðeÞ:

ð21Þ aðH ph ; VÞ

¼ 0 for all V 2 From Eq. (14) we know that b1 ðVÞ  we get explicit formulas for the error in S11 and S12 ,

N ph .

Combining this with Eqs. (18) and (21)

e11 ¼ Cðb1 ðG 1  V ph;1 Þ  aðH ph ; G 1  V ph;1 ÞÞ;

ð22Þ

e12 ¼ þCðb1 ðG 2  V ph;2 Þ  aðH ph ; G 2  V ph;2 ÞÞ:

The choices of V ph;1 2 N ph and V ph;2 2 N ph are unimportant at this point since they do not affect e11 and e12 . However, in the following sections they are important as we shall see. 2.8. Approximations of the dual solution If we had the exact dual solutions G n , Eq. (22) would give the exact errors. However, the exact dual solutions are unknown, and instead we have to use some approximations for them. While it is important that the approximate dual solutions are accurate enough to give good error estimates, it is also very important that they are cheap to compute; ideally the computational cost should be small compared to the cost for solving the original problem. In this section we will present three different methods for computing approximate dual solutions, and we ðBÞ will denote the solutions obtained from these methods with G ðAÞ and G ðCÞ n , Gn n . Eq. (22) shows that we cannot use the finite element method with elements of order p and cell size h (where p now denotes the order used to solve the original problem) to compute the approximate dual solutions, since if G ðA;B;CÞ 2 N ph , the n resulting error estimates will vanish. Instead we choose to work with approximate solutions in the (complete or incomplete) Nedelec space of the next higher order, N hpþ1 . pþ1 Method A. Here we solve the dual problem(s) with finite elements of order p þ 1: Find G ðAÞ such that n 2 Nh

aðG ðAÞ n ; VÞ ¼ bn ðVÞ

8V 2 N pþ1 h :

ð23Þ

The computational cost for this method is, of course, much higher than the cost for the original problem since it involves solving similar problems but with higher order elements. Hence this method would not help in achieving the best accuracy possible with a certain computational cost. In this article we use it mainly as a ‘‘reference’’ method with which the other cheaper and less accurate methods can be compared. Method B. Here we use the same idea as in [12,18]. First we compute approximate dual solutions G ph;n 2 N ph pþ1 and then extend these to functions G ðBÞ n 2 N h . Seeking the dual solutions in the same space as the original solution has the advantage that we can use the same matrix to solve for both the original solution and the dual solutions. If we want error estimates for both S11 and S12 we solve a system with two right-hand sides p

Aph ½gh;1 p gh;1

p gh;2  ¼ ½ bph;1 p gh;2

bph;2 ;

ð24Þ G ph;1

G ph;2 reH ph ‘‘for

and are vectors containing the expansions coefficients (cf. Eq. (15)) for and where spectively. Remember that G ph;1 is identical to H ph , so if we solve this system of equations, we get free’’. Especially when using a direct solver, an extra right-hand side adds very little to the computational cost since the LU -factorization is only done once.

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2605

In the following we use a discrete interpolation operator rph which projects a field in N hpþ1 onto N ph . We define this by saying that the degrees of freedom (i.e. the expansion coefficients for the basis functions) associated with N ph should be preserved. To extend the solutions into N hpþ1 , we note that it is not really the accuracy of G nðBÞ itself that is important p p p ðBÞ ðBÞ (see Eq. (22)) but rather the accuracy of the difference DðBÞ n ¼ G n  V h;n . If we choose V h;n ¼ rh G n , we see that p ðBÞ pþ1 ðBÞ n N ph : DðBÞ n ¼ G n  rh G n 2 N h

Noting that the actual value of rph G ðBÞ is not of primary importance, this lead us to the idea of putting n p ðBÞ rph G ðBÞ n ¼ G h;n . We then compute Dn by solving a defect problem, p aðDðBÞ n ; VÞ ¼ bn ðVÞ  aðG h;n ; VÞ

8V 2 N hpþ1 n N ph :

In order to do this we write the matrix Ahpþ1 as " # e p;pþ1 Aph A pþ1 h Ah ¼ ; e pþ1;p A e pþ1;pþ1 A h

where Aph is space N hpþ1 ,

h

e -matrices are associated with the additional basis functions in the defined as before and the A

e p;pþ1 Þ ¼ ð A e pþ1;p Þ ¼ aðV i ; V jþN p Þ; ðA h h ij ji h e pþ1;pþ1 Þ ¼ aðV iþN p ; V jþN p Þ; ðA h ij h h

i ¼ 1; 2; . . . ; Nhp ;

j ¼ 1; 2; . . . ; Nhpþ1  Nhp ;

i; j ¼ 1; 2; . . . ; Nhpþ1  Nhp :

pþ1 are defined as The vectors gnðBÞ and bh;n " p # " p # bh;n gh;n pþ1 ðBÞ gn ¼ and bh;n ¼ pþ1 ðBÞ b~h;n g~n

where g~nðBÞ are the additional expansion coefficients for the dual solution G nðBÞ and the additional elements in the load vector are defined as pþ1 ðb~h;n Þj ¼ bn ðVjþNhp Þ;

j ¼ 1; 2; . . . ; Nhpþ1  Nhp :

The system of equations we solve to compute g~nðBÞ is e pþ1;pþ1 g~ðBÞ ¼ A e pþ1;p gp þ b~pþ1 : A h;n h;n n n n

ð25Þ

Since most of the global behaviour of the solution is determined by the values of the lower order degrees of freedom, especially now when these are fixed, this system is easy to solve. We use the Jacobi method, where about 5–15 iterations are enough for a tolerance of 102 . This tolerance affects only the relative accuracy of the error estimators, and 102 should be sufficient in most situations. Note that the expansion coefficients p ðBÞ ðBÞ ~nðBÞ alone, for DðBÞ n ¼ G n  rh G n are determined by g " p #  p  " # g 0 g h;n h;n  ; dnðBÞ ¼ ¼ g~nðBÞ 0 g~nðBÞ where 0 indicates vectors with appropriate number of zeros. Method C. Method 3 is a simplification of method 2 where we now make only one Jacobi iteration. This is indeed a very cheap operation. The extra work needed to obtain the error estimators for Method C is

2606

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

(1) The additional work associated with an extra right-hand side in the original problem––if an error estimate for S12 is required (cf. Eqs. (16) and (24)). e pþ1;p and the diagonal of A e pþ1;pþ1 . (2) The work needed to assemble A n n (3) One or two matrix–vector multiplications for the Jacobi iteration, cf. Eq. (25) e pþ1;pþ1 Þ1 ð A e pþ1;p gp þ b~pþ1 Þ: g~nðCÞ ¼ diagð A h;n h;n n n (4) Element-wise evaluation of the inner products in Eq. (22). With the possible exception of the first step, if an iterative solver is used, all of the operations above require work that increases linearly with the number of degrees of freedom, Nhp . The work required to solve the original linear system of equations usually increases faster than that, and therefore the additional work should be less and less important when the problem size increases. 2.9. Computing errors and residuals Approximate error estimates corresponding to Eq. (22) are now obtained from ðX Þ

ðX Þ

ðX Þ

ðX Þ

ðX Þ

ðX Þ

e11 ¼ Cðb1 ðD1 Þ  aðH ph ; D1 ÞÞ; e12 ¼ þCðb1 ðD2 Þ  aðH ph ; D2 ÞÞ; ðX Þ

ðX Þ

ð26Þ

ðX Þ

where D1 ¼ G 1  rph G 1 2 N pþ1 n N ph and X ¼ A, B or C. h To determine which elements to refine, we use the element-wise residuals. These are conventionally computed by splitting the total error in Eq. (26) into element-wise contributions and then integrating by parts over each element, see e.g. [14]. However, from the integration by parts over an element K we get         1 1 1 1 p p p Þ ðX Þ ðX Þ ðX Þ ^ ¼ r  ; D  n  ; D ¼ r  ; D r  H ph ; r  DðX r  H r  H r  H h h h n n n n     K K @K K since DnðX Þ contain only basis functions of order p þ 1, which are Legendre polynomials along the edges and therefore orthogonal, with respect to h; iedge , to lower order functions. Consequently integration by parts is unnecessary and the element-wise residuals can be computed using ðX Þ

n

Þ gK;n ¼ ð1Þ CðbK ðDnðX Þ Þ  aK ðH ph ; DðX n ÞÞ;

where aK ð; Þ and bK ðÞ are the operators corresponding to að; Þ and bðÞ but restricted one single element K instead of the entire domain X. This is advantageous since these quantities are partly known from the asðX Þ sembling of the matrix and right-hand side vector corresponding to að; Þ and bðÞ. The total errors, e11 and ðX Þ ðX Þ ðX Þ e12 are the sum of all the element-wise residuals gK;1 and gK;2 respectively.

3. Results Ref In this section we will compare the computed S-parameters with a pair of reference solutions, S11 and We obtained these reference solutions by extrapolating the results we got using complete second order elements and h-adaptivity to h ¼ 0 (or equivalently Nhp ¼ 1). We estimate that the error of the reference solutions is at least one order of magnitude smaller than the most accurate computed results we compare with. We have also used a mode-matching code and an FDTD code to compute the S-parameters, but these methods required a much larger computational effort to reach the same accuracy. When the accuracy is not so important, both the mode-matching code and the FDTD code work well, but when very high accuracy is Ref S12 .

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2607

1 0.9

Ref

|S 11 | & |SRef | 12

0.8 0.7 0.6 0.5 0.4 0.3 0.2

|S11| |S12|

0.1 0

4.5

5

5.5

6

6.5

7

Frequency [GHz] Ref Ref Fig. 2. Absolute values of S11 and S12 for the waveguide cavity resonator.

Ref Ref needed they converge too slowly. Fig. 2 shows the magnitudes of S11 and S12 in the frequency range 4.3 and 7.0 GHz. The problem is weakly resonant with a resonance frequency of about 5.78 GHz. All simulations start with a uniform triangulation of the domain X shown in Fig. 3. The longest edge in the mesh is w=8, which corresponds to a little less than one tenth of the free-space wavelength at the resonant frequency 5.78 GHz. When comparing the convergence rates of different methods we will plot the errors versus the number of degrees of freedom, Nhp . However, we stress that while the size of the matrix Aph is determined by Nhp , this is not a correct measure for the work required to solve the linear system Aph hph ¼ bph;1 . One must take into account that the number of entries in each row of Aph depend on the order of the method. For iterative solvers also the conditioning of the matrix affects the computational cost. Higher order methods give denser matrices and their conditioning depends on how the element bases are constructed. The asymptotic number of entries per row for different element types are given in Table 3.

Fig. 3. Initial triangulation of X used for all computations.

Table 3 Asymptotic number of elements (boundary effects ignored) per row in the coefficient matrix Aph for different element types p 1 2 3

Entries per row in Aph Incomplete

Complete

5.0 11.6 19.6

10.0 16.6 22.6

2608

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

3.1. Uniform refinement First, we use uniform refinement to compute the convergence rates for both incomplete and complete, Ref first, second and third order, finite elements. In Fig. 4 we show how the error, jS12  S12 j, depend on the number of degrees of freedom at the frequency 5.50 GHz. All elements, except for the complete first order 2=3 elements, show a convergence rate close to h4=3 / ðNhp Þ . This rate is due to the right-angled reentrant corners, more generally we would get the convergence rate h1=ð1/=2pÞ for a metal edge that subtends the angle / in two dimensions. The strange convergence behaviour of the complete first order elements is most likely caused by cancellation between errors coming from different phenomena, e.g. numerical dispersion and singularities at the reentrant corners. We expect that if we could refine the mesh even more, also these elements would reach a convergence rate of ðNhp Þ2=3 . For other frequencies, and for S11 as well, the general pattern is the same; the complete first order elements give accurate results but show a strange convergence behaviour, while the other elements give convergence rates as expected. Note that if an accuracy of 104 is required, which can be the case in the design of e.g. satellite components, about 107 degrees of freedom would be required for the complete third order elements, and more than 108 degrees of freedom for the ordinary incomplete first order elements. This is clearly unrealistic, and would be even more so for three-dimensional problems. 3.2. Adaptive refinement

|S12 – SRef | [–] 12

Ref In Fig. 5 we show how the error jS12  S12 j depends on the number of degrees of freedom ðNhp Þ for the different element types and for the different methods to compute the error estimates. The computations are made at the frequency 5.50 GHz, but the behaviour is similar for all frequencies. We note that all methods (A–C) give almost identical results, even though Methods B and C are much cheaper to compute. Comparing the different element types, we see that the incomplete element types give an initial rapid decrease of the error while the initial uniform mesh is refined to get meshes with correct element size 1 2 distribution, then constant convergence rates of ðNhp Þ / h2 and ðNhp Þ / h4 for the first and second order

10

–1

10

–2

p = 1, incompl. p = 1, compl. p = 2, incompl. p = 2, compl. p = 3, incompl. p = 3, compl.

–3

10

3

10

4

10

5

p

10

Degrees of freedom, Nh [–]

Fig. 4. Total error for S12 when complete and incomplete, first, second and third order, finite elements are used. The dotted line corresponds to N 2=3 .

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616 1:st order incomplete elements

–1

10

2609

Method A Method B Method C

–2

|S12 – S Ref | [–] 12

10

1:st order complete elements

–3

10

2:nd order incom– plete elements

–4

10

–5

10

2:nd order complete elements –6

10

3

4

10

5

10

10

Degrees of freedom, N ph [–]

Fig. 5. Convergence behaviour for the different element types when Methods A, B and C are used to refine the mesh. The dotted lines indicates the slopes ðNhp Þ1 and ðNhp Þ2 .

elements respectively. We note that the convergence rates are higher than the nominal convergence rates in global energy norm, which are h1 and h2 for the two element types respectively. Currently we cannot fully explain this. However, the superconvergence is connected to the anisotropic nature of the elements. If we integrate the normal instead of the tangential component of H over the ports, the convergence rates decrease to h1 and h2 respectively. For the first order complete elements the error actually increases during the initial refinement. The explanation for this is probably that when the elements are refined close to the corner singularities, the cancellation of different errors mentioned in Section 3.1 is lost. However, the correct element size-distri1 bution is reached quite quickly and then the convergence rate is ðNhp Þ as it should. The second order complete elements show a more strange behaviour; initially the error decreases very rapidly until the number of elements, or degrees of freedom, are about twice of what they were in the initial 2 mesh, but after this also these elements give the correct convergence rate of ðNhp Þ . The same errors as in Fig. 5 are shown in the complex plane in Fig. 6, where one can see that the error for the complete second order elements changes ‘‘signs’’ during the refinement process. Thus the strongly non-asymptotic initial convergence seen in Fig. 5 is clearly the result of cancellation between different types of errors.

90

90 complete 1:st order

complete 2:nd order

180 –2 –1 0

270

incomplete 0180 1:st order

0 –5

in –4 2: com nd p –3 or lete de r –2

–1

270

Ref Ref Fig. 6. The error S12  S12 for the frequency 5.50 GHz shown in the complex plane. We have log10 ðjS12  S12 jÞ on the radial axis and Ref Þ on the angular axis. The origin corresponds to 103 in the left diagram and 106 in the right. \ðS12  S12

2610

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616 2.5 Method A Method B Method C

Efficiency index [–]

2 1:st order complete elements

1.5

1 1:st order incomplete elements

0.5

0

3

10

4

5

10

10

Degrees of freedom, N ph [–]

Fig. 7. Efficiency indices for Methods A, B and C when first order elements are used.

Further we note that the number of degrees of freedom needed to get an accuracy of 104 is about 104 for the second order complete elements. Comparing this with what was required when uniform refinement was applied we see that the problem size is reduced by at least a factor 103 . We conclude that higher order elements and h-adaptivity is absolutely necessary when very high accuracy is needed. In Figs. 7 and 8 we show the efficiency indices for first and second order elements respectively. The efficiency index is defined as (the magnitude of) the ratio between the estimated error and the actual error. If the efficiency index is less than one, the error is under-estimated. Studying the efficiency indices for the first order methods, we see that the error is over-estimated for the complete elements and underestimated for the incomplete elements. For the incomplete elements the indices lies in the range 0.75–1.00, which must be considered very good. Also the estimate for the complete elements, which is at most twice the exact error, is good except for small meshes where cancellation of errors occurs. However, we note that Methods A, B and C give somewhat different efficiency indices. Method A gives better estimates, even though the difference is not very big, and the estimates from Methods B and C do not seem to converge to the correct values even though they are quite close. 2.5 Method A Method B Method C

Efficiency index [–]

2

2:nd order complete elements

1.5

1 2:nd order incomplete elements

0.5

0

4

10

Degrees of freedom, Nph [–]

5

10

Fig. 8. Efficiency indices for Methods A, B and C when second order elements are used.

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2611

As shown by Fig. 8, the efficiency index for the incomplete second order elements is less than one, reaching the minimum value of a little less than 0.5 in the beginning of the refinement process and then increasing towards one. For most applications this is good enough, and though the estimates from method B and C do not quite converge to one, they still give good information about the accuracy. The error estimates for the second order complete elements show a more complex behaviour. First they underestimate the error almost by an order of magnitude at the peak. Then the efficiency index rises above one, and when the solution reaches asymptotic behaviour (for Nhp 2  104 ) it approaches unity again. 3.3. Comparison between adaptivity based on e11 and e12 In the previous section we used the local contributions from both e11 and e12 to refine the mesh. Here we will examine how the accuracy of S11 and S12 vary when the mesh refinement is based on e11 only, e12 only or both together. We study two frequencies, 5.50 GHz, where the transmission is quite large and 4.30 GHz, where the transmission is quite low. Method C was used for all mesh refinements in this section.

complete 2:nd order

–4

incomplete 2:nd order

10

–2

complete 1:st order

|S 11 – S Ref 11 | [–]

10

incomplete 1:st order

3.3.1. High transmission Fig. 9 shows the errors of the computed values for S11 at 5.50 GHz. As expected, all element types give better results when the refinement is based on e11 only. However, the difference is small, and this can be understood by noting that since the transmission is quite big, the dual solutions G 1 and G 2 have similar magnitudes everywhere. In Fig. 10 we show jr  Hj ¼ jr  G 1 j (a similar plot for jr  G 2 j is obtained by simply rotating the figure 180 degrees around the center). We note that although the fields are not symmetric with respect to the plane centered between Port 1 and Port 2, they have similar magnitudes on both

e11–refinement e12–refinement (e11+e12)–refinement

–6

10

3

4

10

5

10

Degrees of

freedom, N ph

10

[–]

Fig. 9. Convergence rates for S11 at 5.50 GHz when the mesh refinement is based on e11 only, e12 only and both.

Fig. 10. Contour plot of jr  Hj ¼ jr  G 1 j at 5.50 GHz. The contour levels are separated by 25 A/m2 .

2612

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

Fig. 11. Mesh refined by local contributions to e11 at 5.50 GHz using first order incomplete elements. The mesh consists of 6598 triangles.

sides, and this gives meshes with similar element sizes on both sides. A mesh resulting from e11 -based refinement is shown in Fig. 11, and in this case meshes refined by e12 are very similar. Consequently the accuracy does not depend much on which goal-functional we use.

–2

–6

complete 2:nd order

incomplete 2:nd order

–4

10

complete 1:st order

incomplete 1:st order

|S11 – S Ref | [–] 11

10

e –refinement 11 e –refinement 12 (e +e )–refinement

10

11

12

3

4

10

5

10

Degrees of

10

p freedom, Nh

[–]

Fig. 12. Convergence rates for S11 at 4.30 GHz when the mesh refinement is based on e11 only, e12 only and both.

–2

incomplete 2:nd order

complete 2:nd order

–6

10

complete 1:st order

–4

10

incomplete 1:st order

|S12 – S Ref | [–] 12

10

e –refinement 11 e –refinement 12 (e +e )–refinement 11

12

–8

10

3

10

4

5

10

10 p

Degrees of freedom, N h [–]

Fig. 13. Convergence rates for S12 at 4.30 GHz when the mesh refinement is based on e11 only, e12 only and both.

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2613

Fig. 14. Contour plot of jr  Hj ¼ jr  G 1 j at 4.30 GHz. The contour levels are separated by 25 A/m2 .

3.3.2. Low transmission Figs. 12 and 13 show the errors for S11 and S12 respectively at 4.30 GHz. Unlike the results at 5.50 GHz, there is now a significant difference in accuracy depending on how the meshes were refined. Refinement based on local contributions to e11 give the best results for S11 , but they give the worst results for S12 . Similarly the opposite is true for mesh-refinement based on e12 . The refinement based on both e11 and e12 , seem to be a good compromise in all cases. Note that no matter what we base our refinements on, we finally achieve the correct order of convergence, ðNhp Þ1 or ðNhp Þ2 , but the errors can differ by more than a factor of 5. The difference between the results for different refinement criteria can be understood by looking at the solution in Fig. 14. Both H and G 1 are concentrated close to Port 1, and therefore adaptive meshes based on e11 (which depend on H and G 1 ) are refined mainly close to Port 1. On the other hand, since G 2 is

Fig. 15. Mesh refined by local contributions to e11 (top), e12 (middle), and both e11 and e12 (bottom) at 4.30 GHz. First order incomplete elements were used for all meshes, and they contain 6721, 6981 and 6858 elements respectively.

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

Relative number of elements [–]

2614

1 0.8 0.6 0.4 p = 1, incompl. p = 1, compl. p = 2, incompl. p = 2, compl.

0.2 0

0

0.2

0.4

0.6

0.8

1

Relative area [–]

Fig. 16. Element size distribution: the fraction of the total number of elements, sorted by ascending size is plotted versus the fraction of the total area occupied by these elements.

concentrated close to Port 2, adaptive meshes based on e12 (which depend on H and G 2 ) are refined equally on both sides of the waveguide. Fig. 15 shows meshes refined by e11 , e12 and both. 3.4. Element size distribution

Relative number of elements [–]

Different choices of element type leads to significantly different adaptively refined meshes. The major difference is that meshes for second order elements are more heavily refined close to the singularities. In Fig. 16 we have plotted the fraction of the total number of elements, sorted by ascending size, versus the fraction of the total area occupied by these elements. The dotted line correspond to the same element sizes throughout the whole mesh. The more ‘‘bent’’ a curve is, the more non-uniform are the element sizes for the corresponding method. For example we see that the 20% smallest elements occupy about 0.27%, 0.23%, 1.6% and 5.7% of the total area for incomplete second order, complete second order, incomplete first order and complete first order elements respectively. Evidently, second order elements give meshes with more non-uniform element sizes than first order elements. Meshes for complete order elements use a larger fraction of the elements for the area very close to the singularities than meshes for incomplete order elements do (see Fig. 17, which is a close up of Fig. 16).

0.25 0.2 0.15 0.1 p = 1, incompl. p = 1, compl. p = 2, incompl. p = 2, compl.

0.05 0

0

1

2

3

Relative area [–]

Fig. 17. Zoom of Fig. 16.

4

5 –3

x 10

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

2615

Thus, the more accurate the method is, the smaller are the smallest elements relative to the largest ones. This indicates that the order of the elements is of less importance very close to the singularities, and that is what hp-methods make use of.

4. Conclusion The simulations with uniform mesh refinement showed convergence rates independent of the order of the elements for this geometry containing reentrant corners. This means higher order elements are not effective unless they are combined with adaptive mesh refinement. The goal-oriented error estimates for S11 and S12 derived here make it possible to recover the optimal convergence rates for first and second order, complete and incomplete, curl-conforming elements. Further, the error estimates generally give accurate predictions of the errors and can therefore be used to assess whether computed results are accurate enough. Comparing computationally more and less intensive methods for obtaining dual solutions, we found that well-working error estimates can be computed with only little extra work. Thus, if we want results with an accuracy of about 104 , we can decrease the memory and time consumption by several orders of magnitude for many problems by means of adaptive, instead of uniform, mesh refinement. Finally, comparing results from simulations where the mesh-refinements were based either on the error for S11 or S12 , we saw that meshes that are refined using the correct goal-functional can give substantial benefits. However, both parameters can be efficiently and accurately calculated by using error indicators based on both scattering parameters.

References [1] R. Becker, R. Rannacher, Weighted a posteriori error control in FE methods, Lecture at ENUMATH 95, Paris, August 1995, in: H.G. Bock (Ed.), Proc. ENUMATH 97, World Scientific, Singapore, 1998, pp. 621–637. [2] J.C. Nedelec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315–341. [3] J.C. Nedelec, A new family of mixed finite elements in R3 , Numer. Math. 50 (1986) 57–81. [4] A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Norwood, 1995. [5] J. Jin, The Finite Element Method in Electromagnetics, Wiley, New York, 1993. [6] L. Vardapetyan, L. Demkowicz, hp-adaptive finite elements in electromagnetics, Comput. Methods Appl. Mech. Engrg. 169 (1999) 331–344. [7] W. Rachowics, L. Demkowicz, An hp-adaptive finite element method for electromagnetics––part II: A 3D implementation, Int. J. Numer. Methods Engrg. 53 (1) (2002) 147–180. [8] M. Ainsworth, J. Coyle, Hierarchic hp-edge element families for MaxwellÕs equations on hybrid quadrilateral/triangular meshes, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6709–6733. [9] K. Eriksson, D. Estep, P. Hansbo, Johnson, Introduction to adaptive methods for differential equations, Acta Numerica 4 (1995) 105–158. [10] P. Monk, A posteriori error indicators for MaxwellÕs Equations, J. Comput. Appl. Math. 100 (1998) 173–190. [11] R. Beck, R. Hiptmair, R.H.W. Hoppe, B. Wohlmuth, Residual based a posteriori error estimators for eddy current computation, M2 AN 34 (1) (1999) 159–182. [12] R. Beck, R. Hiptmair, B. Wohlmuth, Hierarchical error estimator for eddy current computation, in: P. Neittaanm€aki, T. Tiihonen (Eds.), ENUMATH 99––Proceedings of the 3rd European Conference on Numerical Mathematics and Advanced Applications, Jyv€ askyl€ a, Finland, July 26–30, World Scientific, Singapore, 2000, pp. 110–120. [13] M. Salazar-Palma, T.K. Sarkar, L.-E. Garcia-Castillo, T. Roy, A. Djordjevic, Iterative and Self-Adaptive Finite-Elements in Electromagnetic Modelling, Artech House, Norwood, 1998. [14] R. Becker, R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numerica 10 (2001) 1–102. [15] D.K. Cheng, Fundametals of Engineering Electromagnetics, Addison-Wesley, Reading, MA, 1993.

2616

P. Ingelstr€om, A. Bondeson / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2597–2616

[16] J.P. Webb, Hierarchical vector basis functions of arbitrary order for triangular and tetrahedral finite elements, IEEE. Trans. Antennas Propagat. 47 (8) (1999) 1244–1253. [17] A. Bondeson, Y. Yang, P. Weinerfelt, Shape optimization for radar cross sections by a gradient method, Int. J. Numer. Methods Engrg., submitted for publication. [18] R. Bank, Hierarchical bases and the finite element method, Acta Numerica 5 (1996) 1–43.