Numerical subgrid scale models for the Yee scheme

for the Yee scheme

B. Engquist

, J. H¨ aggblad

, and O. Runborg

1 Department of Mathematics and Institute for Computational Engineering and Sciences, The University of Texas at Austin, 1 University Station C1200, Austin TX

78712, U.S.A.

2 Department of Numerical Analysis, CSC, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden.

3 Swedish e-Science Research Centre, Sweden.

4 jonhagg@kth.se

Abstract. The Yee scheme is a very common and practical algorithm for the simulation of wave propagation on uniform grids. We develop numer- ical subgrid scale models in order to incorporate effects of obstacles and holes that are smaller than the grid spacing. The models are based on pre- computing at the microscale, and are thus including the effect of the detailed small scale shape. Numerical examples in 1D, 2D and 3D are given.

1 Introduction

The finite-difference time-domain (FDTD) method was first introduced by Yee in 1966 [17], and is hence sometimes referred to as the Yee scheme. It is a simple and efficient second order accurate method that was originally formulated for solving Maxwell’s equations [12, 13], but can easily be applied to various forms of wave propagation problems [10], such as acoustic scattering from a rough ocean sur- face [11] and room acoustics [3]. It is based on centred difference approximations on a Cartesian staggered grid. This gives the method relatively low dispersion errors, important for wave propagation, and a very low memory footprint since no explicit grid is stored and the field components are not represented at every discretization point.

Consider the system formulation of the two-dimensional wave equation,

∂p

∂t = a∇ · u, ∂u

∂t = b∇p, (1)

where u = (u, v). This equation turns up in various forms in many fields where

wave-like phenomena occur, such as acoustics and electromagnetics. For a = −κ

and b = −1/ρ, we obtain the system of equations governing the propagation of

sound waves in acoustics,

∂p

∂t + κ∇ · u = 0, (2)

ρ ∂u

∂t + ∇p = 0, (3)

where p(x, y, t) represents the pressure and the vector u(x, y, t) the velocity field.

The coefficients specify density ρ, and bulk modulus κ of the medium. The equation (2) is the mass balance, and (3) is the momentum balance.

If on the other hand we set a = 1/, b = 1/µ together with p = E

, u = H

, v = −H

, we get the TM mode,

 ∂E

∂t = ∂

H

∂x − ∂H

∂y , µ ∂H

∂t = − ∂E

∂y , µ ∂H

∂t = ∂E

∂x ,

and with p = H

, u = −E

, v = −E

, we get the TE mode,

µ ∂H

∂t = ∂E

∂y − ∂E

∂x ,

 ∂E

∂t = ∂H

∂y ,

 ∂E

∂t = − ∂H

∂x ,

respectively, of Maxwell’s equations in vacuum. The coefficients  and µ denote the permittivity and permeability, respectively.

We shall use the formulation (1), and refer to it as the acoustic equations. Later we shall also consider the one-dimensional form as well as a generalization to three dimensions.

The Yee scheme for the system (1) becomes

p

2,l+¹₂

= p

2,l+¹₂

+ ∆ta



D

u

+ D

v



, (4)

u

j,l+¹₂

= u

+ ∆tb

2

D

p

j+¹₂,l+¹₂

, (5)

v

2,l

= v

2,l

+ ∆tb

2,l

D

p

1 2

j+¹₂,l+¹₂

, (6)

where

D

u

2

= u

2

− u

h

,

D

v

2,l

= v

2,l+1

− v

h

,

D

p

2,l+¹₂

= p

2,l+¹₂

− p

2,l+¹₂

h

, D

p

2,l+¹₂

= p

2,l+¹₂

− p

h

.

The field coordinates are given by (x

, y

) = (jh

, lh

) and the time levels are at t

= n∆t. The indices are allowed to be half-integer, and correspond to a field point halfway in between the integer points. We will assume that h

= h

= h. Thus we can think of the spatial grid placement as cells with p

2,l+¹₂

placed in the centre.

Having regions in the computational grid with geometric details smaller than the discretization length h gives rise to the problem of subgrid effects. When the subgrid phenomena is of a simple type, analytical approximations are possible.

Common in the computational electromagnetics (CEM) community is to use con- tour path models [16, 14]. Another possibility is local mesh refinement [2, 1, 9, 5, 4, 8]. Local mesh refinement is however typically computationally costly and may give difficulties with stability and mesh interface reflections.

We are here developing methods with the goal of having the accuracy and the generality of local mesh refinements but the efficiency of analytic formulas fitting the Yee scheme. It is based on the framework provided by the heterogeneous multiscale method (HMM) [6]. A microscale simulation is used to generate effective terms in a macroscale scheme. The microscale simulation can be done on-the-fly during the macroscale computation or as in our case for linear problems in a pre-computation.

Our technique can be seen as a type-A HMM method.

In one dimension we consider the case of an obstacle, defined by a local reduction of wave speed, with extent smaller than the discretization size h. In two and three dimensions we consider a small hole of arbitrary shape in a perfectly reflecting plate of suitable dimension.

If we model this in a naive way by letting the size of the obstacle or gap be

as large as the grid spacing, typically either too little or too much energy is trans-

mitted, depending on whether we are dealing with an obstacle or hole. To better

model the effective behaviour of such subgrid geometric features, we take a linear

combination of absorbing and reflecting boundary conditions, together with source

terms on both sides, to transmit a portion of the waves through. We start with the

one-dimensional problem, before considering higher dimensions. In one dimension

we also investigate the behaviour of a first order correction term, which is necessary

if we want a more accurate frequency dependent behaviour.

By pre-computing an effective stencil, we manage to capture the leading order effects accurately with only minor added complexity to the macro scale scheme, and essentially no performance loss. For the case of very low frequencies in one dimension, we also demonstrate the possibility to correctly model higher order behaviour.

2 Subgrid obstacles in one dimension

We start out by considering subgrid size obstacles in one dimensions. Let an ob- stacle of width w

be placed with centre at x = 0. We then define the macroscale to be a coarse discretization of −1 ≤ x ≤ 1. By coarse we mean that the dis- cretization length h = 2/N is h  w

. The corresponding microscale problem is the discretization of −δ ≤ x ≤ δ, such that δ > w

/2 and ˜ h = 2δ/ ˜ N  w

, i.e., the obstacle is sufficiently well resolved to enable an accurate computation of interacting waves.

The effects of the subgrid obstacle are modelled by dividing the macroscale into two independent regions, x < 0 and x > 0. These two we then connect at x = 0 through the use of connected boundary conditions. The boundary conditions are controlled by parameters which we can use to set the amount of reflection and transmission by letting the numerical boundary be a linear combination between reflecting and absorbing boundary conditions. These coefficients are computed sep- arately by solving the corresponding micro problem.

2.1 Effective stencils on the macro scale

First we derive a set of effective stencils that allow us to incorporate the effects computed at the microscale, to the macroscale.

The system (1) when reduced to one space dimension, is given by

p

= au

, (7)

u

= bp

(8)

where we will assume a, b < 0, which is the case in acoustics. The domain is defined as x ∈ [−1, 1]. The characteristic quantities of the continuous system (7)–(8), are

p + Zu, ← (9)

p − Zu, → (10)

where Z = −pa/b. These two expressions travel solely to the left and to the right,

respectively, and can be obtained through the diagonalisation of (7)–(8).

Leading order effect

Assume then that the obstacle placed at x = 0 divides the domain into two regions.

We use the characteristic quantities (9)–(10) to connect the two. Thus if we let the left going wave on the left side be a linear combination of the right going wave, i.e., a wave that hits the obstacle from the left and then reflects, and a transmitted left going wave hitting the obstacle from the right, then we obtain

p

+ Zu



= α p

− Zu



+ β p

+ Zu

 

.

The superscripts refer to which domain the corresponding variables belong to. Here α and β correspond to the magnitude of the reflection and transmission, respec- tively. Hence if we rearrange,

(α + 1) Zu

= (α − 1) p

+ β p

+ Zu

 

, and solve for u

, we get

u

= α − 1 α + 1

1 Z p

+ β α + 1

 1

Z p

+ u

 

. (11)

On the right side we obtain from a similar argument, p

− Zu



= α p

+ Zu



+ β p

− Zu

 

. Thus

− (α + 1) Zu

= (α − 1) p

+ β p

− Zu

 

, and solving for u

gives us

u

= − α − 1 α + 1

1 Z p

− β α + 1

 1

Z p

− u

 

. (12) What remains is to incorporate these expressions, (11), (12), into the discretized formulation. The standard Yee discretization of (7)–(8) is given by

p

= p

+ aλ u

− u

, u

= u

+ bλ 

p

j+¹₂

− p

 ,

where λ = ∆t/h. The domain is discretized using h = 2/N , which gives the co- ordinates x

= jh, j = −N/2, . . . , N/2, for u

, and x

2

= (j + 1/2)h, j =

−N/2, . . . , N/2 − 1 for p

2

. Thus the grid is such that u

is placed at x = 0, while

p

and p

are placed directly to the left and right, respectively. The velocity

component u

then becomes multi-valued, which we denote by superscripts, u

and u

. Hence a natural discretization of (11), (12) is

(u

)

= α − 1 α + 1

1 Z p

1 2

−¹₂

+ β α + 1

 1 Z p

1 2 1 2

+ (u

)



, (13)

(u

)

= − α − 1 α + 1

1 Z p

2

− β

α + 1

 1 Z p

−¹₂

− (u

)



. (14)

Alternatively, we can reorder these expressions to obtain something that resembles a typical derivative discretization,

(u

)

= β

α + 1 (u

)

+ 1 α + 1

1 Z

 βp

1 2 1 2

− (1 − α)p

 , (u

)

= β

α + 1 (u

)

+ 1 α + 1

1 Z



(1 − α)p

1 2 1 2

− βp

 .

These we can then use in the update stencil for the points adjacent to the obstacle, p

−¹₂

= p

−¹₂

+ aλ (u

)

− u

, (15) p

1 2 1 2

= p

1 2 1 2

+ aλ u

− (u

)

. (16) If we discretize implicitly,

(u

)

= 1 Z p

1 2

−¹₂

+ β α + 1

 1 Z p

1 2 1 2

+ (u

)

 , (u

)

= − α − 1

α + 1 1 Z p

2

− β

α + 1

 1 Z p

−¹₂

− (u

)

 , we obtain

(u

)

= θ

− θ

1 − θ

1 Z p

1 2

−¹₂

+ θ

− θ

θ

1 − θ

1 Z p

1 2 1 2

,

(u

)

= − θ

− θ

1 − θ

1 Z p

1 2 1 2

− θ

− θ

θ

1 − θ

1 Z p

1 2

−¹₂

, where

θ

= α − 1

α + 1 , θ

= β α + 1 . This discretization is however only stable when β ≤ α.

Thus we have produced local update stencils for x = 0 that to leading order can

model the reflection and transmission of an obstacle much smaller than the grid

spacing on the macroscale.

Higher order correction

That the reflection and transmission is frequency dependent is well known, see e.g. [7]. Especially, in the limit when the frequency tends to zero or, equivalently, the width of the obstacle w

goes to zero, then to lowest order the reflection and transmission goes to zero and one, respectively. Thus one needs the next term in the expansion in frequency to even capture qualitatively the reflection behaviour.

Hence to model a potential frequency dependency in the reflection and trans- mission we need to include higher order terms, i.e., we need to include derivatives.

First we consider these on their own, without the lowest order expressions, before stating the full expressions. On the left side we have the additional terms

p

+ Zu



= −γ ∂

∂x p

− Zu



+ µ ∂

∂x p

+ Zu

 

, and on the right side we have

p

− Zu



= γ ∂

∂x p

+ Zu



− µ ∂

∂x p

− Zu

 

.

In these expressions γ corresponds to first order reflection and µ to first order transmission.

A straight-forward discretization of these additional terms leads to the expres- sions

p

1 2

−¹₂

+ Z(u

)

= α  p

1 2

−¹₂

− Z(u

)



− γ h

 p

1 2

−¹₂

− p

−³₂

− Z u

− u

 + β 

p

1 2 1 2

+ Z(u

)

 + µ

h

 p

1 2 3 2

− p

+ Z (u

− u

)  ,

p

1 2 1 2

− Z(u

)

= α  p

1 2 1 2

+ Z(u

)

 + γ

h

 p

1 2 3 2

− p

+ Z (u

− u

)  + β 

p

1 2

−¹₂

− Z(u

)



− µ h

 p

1 2

−¹₂

− p

−³₂

− Z u

− u



,

which we can rewrite as (u

)

= α − 1

α + 1 1 Z p

1 2

−¹₂

− γ/h α + 1

 1 Z

 p

−¹₂

− p

−³₂

 − u

+ u



+ β

α + 1

 1 Z p

1 2 1 2

+ (u

)



+ µ/h α + 1

 1 Z

 p

1 2 3 2

− p



+ u

− u

 ,

(17)

(u

)

= − α − 1 α + 1

1 Z p

1 2 1 2

− γ/h α + 1

 1 Z

 p

1 2 3 2

− p

 − u

+ u



− β

α + 1

 1 Z p

1 2

−¹₂

− (u

)



+ µ/h α + 1

 1 Z

 p

−¹₂

− p

 − u

+ u

 .

(18)

As before we use this in the standard update stencil for the pressure (15)–(16).

Including derivatives results in stencils that are wider than the Yee stencil, thus straying slightly outside the focus of the current context of efficient and simple modifications to the standard FDTD-method. Furthermore, the discretization in (17)–(18) is only stable for CFL = 1. Hence its inclusion is mostly to show the po- tential of frequency dependent correction, rather than that of a practical algorithm.

2.2 Computing the effective coefficients

Having obtained effective update stencils for the macro scheme, we turn to the problem of computing the coefficients that the stencils depend on. This is accom- plished by resolving the area around the obstacle in a small micro problem, which we only need to solve once.

Leading order coefficients

Let the domain of the micro problem be −δ ≤ x ≤ δ, where δ is large enough to cover the extent of the obstacle. At both boundaries, x = ±δ, we specify absorbing boundary conditions, and on the left we also add a source term F (x

, t

),

u

= 1 Z p

1 2

−^N₂+¹₂

+ F (x

, t

), u

= − 1 Z p

1 2 N

2−¹₂

.

10⁻³ 10⁻² 10⁻¹ 10⁰ 0

0.2 0.4 0.6 0.8 1

k ⋅ width

α² β²

Fig. 1. The transmitted an reflected amplitude in the micro problem, as a function of frequency multiplied with the width of the obstacle. The obstacle is given by a Gaussian function.

The source term is set to generate a slowly varying harmonic wave, F (x, t) = 2e

.

Thus if we let the computation reach a stationary behaviour we can subtract away the incoming wave F and take the absolute value to reveal the reflected and trans- mitted amplitude. These then specify α and β.

To obtain a specified desired accuracy in the measured α and β, we can specify a limit in the sample variance for both time and space. In the absence of a power source or absorbing materials we also expect to observe energy conservation, which manifests itself by the relation α

+ β

= 1. This can also be used as a criterion for the accuracy of the measured amplitudes.

Higher order terms in the low frequency limit

In the limit when k → 0, then the reflected amplitude goes to zero, as seen in Fig. 1.

Thus to correctly capture the behaviour, we need to incorporate the next term in the expansion in k. This means taking the frequency dependency into account.

To see how the frequency dependency comes into play we consider the second

order formulation of the wave equation with an incoming harmonic wave with am-

plitude A

at the left boundary,

u

− ∂

a(x)

∂

u = 0, x ∈ (0, δ), u

− a(x)u

= 2ikA

, x = 0,

u

+ a(x)u

= 0, x = δ.

For a fixed frequency then we obtain the Helmholtz equation

∂

a(x)

∂

u + k

u = 0, x ∈ (0, δ), iku − a(x)u

= 2ikA

, x = 0,

iku + a(x)u

= 0, x = δ.

Assume we can expand u(x, k) in k, u(x, k) = v

+ kv

+ k

v

/2 + · · · , where we denote

v

(x) = ∂

u(x, 0).

These terms in the k expansion solves a hierarchy of Laplace equations. For v

we simply have

∂

a

(x)∂

v

= 0, x ∈ (0, δ),

∂

v

= 0, x = 0,

∂

v

= 0, x = δ.

For v

then,

∂

a

(x)∂

v

= 0, x ∈ (0, δ), (19)

iv

− a(x)∂

v

= 2iA

, x = 0, (20)

iv

+ a(x)∂

v

= 0, x = δ, (21)

while for ` ≥ 2,

∂

a

(x)∂

v

+ `(` − 1)v

= 0, x ∈ (0, δ) ,

`iv

− a(x)∂

v

= 0, x = 0,

`iv

+ a(x)∂

v

= 0, x = δ.

We note that the solution for v

is a constant, but it is undetermined by the equations since we have only Neumann conditions. On the other hand, the equation for v

only has a solution if

0 = Z

0

∂

a

(x)∂

v

dx

= a(δ)

∂

v

(δ) − a(0)

∂

v

(0)

= i(2a(0)A

− a(0)v

− a(δ)v

).

This gives the solution as

v

= A

2a(0) a(0) + a(δ) . To close the equation for v

we similarly need

−2δv

= Z

0

∂

a

(x)∂

v

dx

= a(δ)

∂

v

(δ) − a(0)

∂

v

(0)

= −i2 (a(δ)v

(δ) + a(0)v

(0)) . (22) The equation (19) gives

∂

v

= C a

(x) , and the boundary conditions (20)–(21) gives

∂

v

(0) = iv

− 2iA

a(0) = −2iA

a(δ) a(0) (a(0) + a(δ)) ,

∂

v

(δ) = − iv

a(δ) = −2iA

a(0) a(δ) (a(0) + a(δ)) , hence

C = −A

2ia(0)a(δ) a(0) + a(δ) . Integrating once more gives

v

(x) = Z

0

C

a(y) dy + D, thus if we insert into (22) we obtain

−2δv

= −i2 a(δ) Z

0

C

a(y) dy + D (a(δ) + a(0))

! ,

and

D = δv

− ia(δ)C R

0

a(y)

dy i (a(0) + a(δ)) . Therefore we have the solution

v

(x) = −A

2ia(0)a(δ)

(a(0) + a(δ))

a(0) Z

0

dy

a(y) − a(δ) Z

x

dy a(y) + δ

!

. (23)

For the case a(0) = a(δ) = a

, a

= 1, we can simplify to v

(x) = −A

ia

2 Z

0

dy a(y) −

Z

dy a(y) + δ

! , and

v

(0) = A

ia

2 Z

0

dy a(y) − δ

! ,

v

(δ) = −A

ia

2 Z

0

dy a(y) + δ

! .

Returning to the problem of computing reflection and transmission coefficients.

The reflection and transmission is given by R(k) = u

u

= u(0, k) − u

(0, k) u

(0, k) , T (k) = u

u

= u(δ, k) u

(δ, k) , which simplifies to

R(k) = (v

+ kv

) − A

A

+ O(k

), T (k) = v

+ kv

A

+ O(k

).

Thus if we identity coefficients in the expansions R(k) = α + ikγ + O(k

), T (k) = β + ikµ + O(k

), we obtain

α = v

− A

A

= 2a(0)

a(0) + a(δ) − 1 = a(0) − a(δ) a(0) + a(δ) , β = v

A

= 2a(0) a(0) + a(δ) , γ = v

(0)

iA

= a

2

Z

dy a(y) − δ

! ,

µ = v

(δ) iA

= − a

2

Z

dy a(y) + δ

! .

We note that to leading order transmission and reflection only depend on the values

of a(x) at the endpoints of the cells, not on the precise shape of the obstacle.

2.3 Numerical examples

To verify that the method works as intended, we present a couple of numerical tests.

Single dominant frequency

When a single frequency is dominating the spectrum, and if it is not too small, it is enough to only include the lowest order reflection/transmission coefficients. From Fig. 1 we see that roughly 10

< kw

< 10

is a valid range.

Let the incoming wave on the macro scale be a right-going modulated gaussian pulse, set as initial conditions

p(x, t) = u(x, t) = exp



− (x − ct)

1/10



exp (i(kx − ωt)),

where c = √

ab. The pulse is started at x = −1/2 and the frequency is specified to k = 20π. The macroscale −1 ≤ x ≤ 1 is discretized using N = 256 cells, giving h = 2/N . The microscale is set to be −h ≤ x ≤ h, and discretized using ˜ N = 512 cells.

We compute two cases, either where the obstacle is defined by a(x) = a

 1 − 7

10 exp



− x

1/1024



, b = b

, (24)

or where the obstacle is defined by b(x) = b

 1 − 7

10 exp



− x

1/1024



, a = a

. (25)

In both cases we have a

= b

= −1.

To compute the effective coefficients α and β, we run the micro problem until we obtain a stationary solution with standard deviation σ < 10

in time. The coefficients are computed by taking the arithmetic mean of the amplitude in the intervals −3h/4 < x < −h/4 and h/4 < x < 3h/4, where the incoming field on the left haft interval is subtracted away before computing the amplitude. As a further check, we verify that |α

+ β

− 1| < 10

.

The resulting reflection and transmission is seen in Fig. 2–3, where we compare against a fully resolved reference solution computed using N = 16384 cells. We see that the resulting effective model captures the behaviour to a large extent.

Broad spectrum reflection

For a broad spectrum pulse we need to take into account the frequency dependency

of the reflection/transmission. This means including the derivative terms specified

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5 1

x

p

EFF HI

−1 −0.5 0 0.5 1

−0.2

−0.15

−0.1

−0.05 0 0.05 0.1 0.15 0.2

x

p error

Fig. 2. The reflection and transmission of a wave packet with a single dominant frequency.

The effective and high-resolution solution is computed using 256 and 16384 cells, respec- tively. The obstacle is defined by (24). The error in the computed energy is εref< 9 · 10⁻³, εtrans< 4 · 10⁻².

by γ and µ. As before, we consider two cases where the obstacle is given either by a(x) = a

 1 − 9

10 exp



− x

1/2048



, b = b

, (26)

or

b(x) = b

 1 − 9

10 exp



− x

1/2048



, a = a

, (27)

with a

= b

= −1. At t = 0 we initialize with a Gaussian pulse centred on x = −1/2 and travelling to the right,

p(x, t) = u(x, t) = exp



− (x − ct)

1/10

 .

By solving the micro problem and finding the stationary solution, we find that α = 0 and β = 1, as expected. The next terms γ and µ are obtained by evaluating the integral (23). Running until t = 1 we then compare with a high resolution computation with N = 16384 to estimate the error, which is shown in Fig. 4. The convergence of the macro scheme as a function of the grid resolution is seen in Fig. 6.

3 Subgrid holes in two and three dimensions

We extend the model to higher dimensions by considering a subgrid hole in a

perfectly reflecting plate. We assume that the plate is along one coordinate axis,

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5 1

x

p

EFF HI

−1 −0.5 0 0.5 1

−0.3

−0.2

−0.1 0 0.1 0.2

x

p error

Fig. 3. Same as Fig. 2, but here the obstacle is defined by (25). The error in the computed energy is εref< 2 · 10⁻², εtrans< 5 · 10⁻².

such that the main transmission is perpendicular to the plate, and along the other coordinate axis. For a plate with a general orientation the method described in this paper can be combined with the boundary condition technique developed in [15].

3.1 Leading order effect

Consider a plate along x = 0, y ∈ R, with a small opening centred at y = −h/2, where h is the grid spacing in the discretization of the macro problem. Assume also that the diameter is smaller than h, and that the grid is aligned such that the u field variable is placed in the middle of the opening. The shape of the hole can be arbitrary. Now consider the update stencil for u given by (5),

u

= u

2

+ ∆tb

2

D

p

1 2

j+¹₂,l+¹₂

.

We observe that only neighbouring points along the x-axis are used in the update.

Thus for the u point centred in the hole we replace the update stencil with the corresponding one-dimensional update formulas (13)–(14). These become

 u

2



= α − 1 α + 1

1 Z p

1 2

−¹₂,−¹₂

+ β α + 1

 1 Z p

1 2 1

2,−¹₂

+  u

2



 ,

 u

2



= − α − 1 α + 1

1 Z p

1 2 1

2,−¹₂

− β α + 1

 1 Z p

1 2

−¹₂,−¹₂

−  u

2



 , where we then use u

in the update formula for p

2,¹₂

, and u

in the update formula for p

2,−¹₂

. Simply put, we are solving the one-dimensional problem perpendicular

to the plate, to compute the transmitted waves.

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5 1

x

p

EFF HI

−1 −0.5 0 0.5 1

−0.1

−0.05 0 0.05 0.1 0.15

x

p error

Fig. 4. The computed reflection using pre-computed coefficients, as well as a fully resolved reference solution. The obstacle is defined by (26). The error in the computed energy for the effective stencil is εref= 2 · 10⁻², εtrans= 2 · 10⁻².

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5 1

x

p

EFF HI

−1 −0.5 0 0.5 1

−0.1

−0.08

−0.06

−0.04

−0.02 0 0.02 0.04 0.06

x

p error

Fig. 5. Same as 4, but the obstacle is defined by (27). The error in the computed energy for the effective stencil is εref= 4 · 10⁻³, εtrans= 2 · 10⁻².

10

10

10

10

10

10

N L

error

Fig. 6. Convergence of the first order approximation of the reflection and transmission.

3.2 Finding the effective coefficients by solving an inverse problem

The two parameters α and β are pre-computed by solving a fully resolved micro problem, only incorporating the hole. The incoming wave is assumed to have a single dominant frequency. Unlike the situation in one dimension, the amplitude of the reflected and transmitted wave does not in a simple way correspond to the values of α and β. Thus to compute these we solve an inverse problem where we find the α and β that best approximate the fully resolved micro problem.

The problem of parameter-matching α and β depends on the objective function.

Another issue that crop up in dimensions higher than one is that since the CFL number must be less than unity, we introduce dispersion errors. Since there are significant dispersion errors when comparing a fully resolved solution with one where a coarse grained discretized is used, we choose to minimize the error in the transmitted and reflected energy, i.e., the L

-norm of the field.

The objective function can be chosen in many other ways, depending what is of interest in the specific case.

3.3 Numerical example

Let the computational domain be [−1, 1] × [−1, 1], discretized on the macroscale

by N = 64 cells in both x and y. Place a plate at x = 0 with a hole of diameter

w = 1/200 centered at y = −h/2, where h = 2/N , giving w/h = 0.16 < 1. This is

done in practice by setting u = 0. Let the incoming wave be generated by a source

x

y

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5 1

−0.2

−0.15

−0.1

−0.05 0 0.05 0.1 0.15 0.2

x

y

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5 1

−0.2

−0.15

−0.1

−0.05 0 0.05 0.1 0.15 0.2

Fig. 7. Illustrating the two-dimensional computation for a hole with relative width w/h = 0.16. The color map is highly compressed to show the transmitted wave. Top: standard Yee scheme computed at a very high resolution (1600² cells). Bottom: effective stencil with α = 0.16, β = 0.84 (64² cells).

1.5 2 2.5 3

−0.2

−0.1 0 0.1 0.2 0.3 0.4

t

u

EFF HI

Fig. 8. The field values at x = 0.156 as a function of time when computed with 64² cells.

We compare against a fully resolved macro problem with 1600²cells.

term F on the x = −1 boundary,

F (x, y, t) = 2 exp (i(k(x − 1) − ωt)) exp



− (t/L − 1)

1/2

 .

The u-field is shown in Fig. 7, first for a computation with very high resolution, and then a low resolution computation where the width of the hole is less then a fifth of the grid resolution (w/h = 0.16). In Fig. 8 we plot the field directly after the hole as a function of time. We see that the error in the transmitted amplitude is similar in size to the dispersion error generated by the coarse grid. To see to which extent the method correctly transmits energy, we plot the error of the computed L

-energy of the transmitted wave as a function of time in Fig. 9. The improvement appears to be roughly one order of magnitude.

All of these confirm that the effective coefficients correctly capture the leading order effect of a small hole. A potential future aspect is to study the sensitivity of the method to the angle of incidence.

3.4 Illustrating the extension to three dimensions

Analogous to the two-dimensional case we can extend the method to three dimen-

sions. Since, unlike in two dimension, there is not a direct substitution that can take

us between the acoustic equations and Maxwell’s equations, we run two separate

computations.

1.5 2 2.5 3 10

10

10

10

t EFF LOW

Fig. 9. The error of the L²value of the transmitted wave as a function of time.

For the three-dimensional acoustic equations, which is expression (1), but with u = (u, v, w) and all variables defined on R

, we run a test case of a small hole in a rigid plate. An incoming harmonic wave is generated at x = −1, with the plate positioned at x = 0, covering the entire zy-plane. The resulting pressure field is given in Fig. 10 for 48

cells, where only the transmitted field is shown. We see that the scheme qualitatively captures the effect of having a hole with diameter smaller than h.

For Maxwell’s equations we consider a perfect electric conductor (PEC) plate in the yz-plane positioned at x = 0. On the plate the tangential electric field is set to zero, i.e., E

= E

= 0. In the middle of the plate there is a small slip along the y-axis at z = 0. An incoming harmonic wave polarized along the y-axis is generated at x = −1. The electric field component E

is shown in Fig. 11 for 48

cells.

While these two test cases illustrate that the numerical scheme on the macro scale is feasible, further work is needed to quantify the accuracy.

4 Conclusions

We have presented a numerical subgrid method for the Yee scheme using precom-

puted effective stencils. The precomputing is performed on a micro grid resolving

the subgrid scale details. The method is general in that no detailed assumptions

are made, and is simple to use in existing FDTD-solvers. Numerical examples are

provided for one, two, and three dimensions. It is shown that the method gives the

Fig. 10. Plot of a three-dimensional computation with 48³ cells. Top: standard Yee. Bot- tom: effective stencil with α = 0.5, β = 0.5.

Fig. 11. Plot of a three-dimensional computation with 48³ cells. Top: standard Yee. Bot- tom: effective stencil with α = 0.1, β = 0.9.

correct energy transmission for a small hole in a plate and does not increase the computational complexity over that of the standard coarse grid scheme.

References

1. M. Berger and P. Colella. Local adaptive mesh refinement for shock hydrodynamics.

J. Comput. Phys., 82(1):64–84, 1989.

2. M. J. Berger and J. Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. J. Comput. Phys., 53(3):484–512, 1984.

3. D. Botteldooren. Finite-difference time-domain simulation of low-frequency room acoustic problems. J. Acoust. Soc. Am., 98(6):3302–3308, 1995.

4. E. B´ecache, P. Joly, and J. Rodriguez. Space–time mesh refinement for elastodynamics.

numerical results. Comput. Method Appl. M., 194(2–5):355–366, 2005.

5. F. Collino, T. Fouquet, and P. Joly. Conservative space-time mesh refinement methods for the FDTD solution of Maxwell’s equations. J. Comput. Phys., 211(1):9–35, 2006.

6. W. E, B. Engquist, X. Li, W. Ren, and E. Vanden-eijnden. Heterogeneous multiscale methods : A review. Commun. Comput. Phys., 2(3):367–450, 2007.

7. J. Fouque, J. Garnier, G. Papanicolaou, and K. Solna. Wave Propagation And Time Reversal in Randomly Layered Media. Applications of Mathematics. Springer, 2007.

8. N. Jansson, J. Hoffman, and M. Nazarov. Adaptive simulation of turbulent flow past a full car model. In State of the Practice Reports, SC ’11, pages 20:1–20:8, New York, NY, USA, 2011. ACM.

9. P. Joly and J. Rodriguez. An error analysis of conservative space-time mesh refinement methods for the one-dimensional wave equation. SIAM J. Numer. Anal., 43(2):825–

859, 2005.

10. J. G. Maloney and K. E. Cummings. Adaptation of FDTD techniques to acoustic modeling. In 11th Annual Review of Progress in Applied Computational Electromag- netics, volume 2, pages 724–731, Monterey, CA, 1995.

11. R. A. Stephen. Modeling sea surface scattering by the time-domain finite-difference method. J. Acoust. Soc. Am., 100(4):2070–2078, 1996.

12. A. Taflove. Application of the finite-difference time-domain method to sinusoidal steady-state electromagnetic-penetration problems. IEEE Trans. Electromagn. C., EMC-22(3):191–202, 1980.

13. A. Taflove and S. C. Hagness. Computational Electrodynamics: The Finite-Difference Time-Domain Method. Artech House Publishers, 3 edition, 2005.

14. A. Taflove, K. Umashankar, B. Beker, F. Harfoush, and K. Yee. Detailed FD-TD analysis of electromagnetic fields penetrating narrow slots and lapped joints in thick conducting screens. IEEE Trans. Antennas Propag., 36(2):247–257, 1988.

15. A.-K. Tornberg and B. Engquist. Consistent boundary conditions for the Yee scheme.

J. Comput. Phys., 227(14):6922–6943, 2008.

16. K. Umashankar, A. Taflove, and B. Beker. Calculation and experimental validation of induced currents on coupled wires in an arbitrary shaped cavity. IEEE Trans.

Antennas Propag., 35(11):1248–1257, 1987.

17. K. S. Yee. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag., 14:302–307, 1966.