NASA/TM{2001{210862

An Implicit Characteristic Based

Method for Electromagnetics

John H. Beggs

Langley Research Center, Hampton, Virginia

W. Roger Briley

Mississippi State University, Mississippi State, Mississippi

May 2001
The NASA STI Program Oce ... in Prole

Since its founding, NASA has been dedicated

to the advancement of aeronautics and space

science. The NASA Scientic and Technical

Information (STI) Program Oce plays a key

part in helping NASA maintain this

important role.

The NASA STI Program Oce is operated by

Langley Research Center, the lead center for

NASA's scientic and technical information.

The NASA STI Program Oce provides

access to the NASA STI Database, the

largest collection of aeronautical and space

science STI in the world. The Program Oce

is also NASA's institutional mechanism for

disseminating the results of its research and

development activities. These results are

published by NASA in the NASA STI Report

Series, which includes the following report

types:

 TECHNICAL PUBLICATION. Reports of

completed research or a ma jor signicant

phase of research that present the results

of NASA programs and include extensive

data or theoretical analysis. Includes

compilations of signicant scientic and

technical data and information deemed

to be of continuing reference value. NASA

counterpart of peer-reviewed formal

professional papers, but having less

stringent limitations on manuscript

length and extent of graphic

presentations.

 TECHNICAL MEMORANDUM.

Scientic and technical ndings that are

preliminary or of specialized interest,

e.g., quick release reports, working

papers, and bibliographies that contain

minimal annotation. Does not contain

extensive analysis.

 CONTRACTOR REPORT. Scientic and

technical ndings by NASA-sponsored

contractors and grantees.

 CONFERENCE PUBLICATION.

Collected papers from scientic and

technical conferences, symposia,

seminars, or other meetings sponsored or

co-sponsored by NASA.  SPECIAL PUBLICATION. Scientic,

technical, or historical information from

NASA programs, pro jects, and missions,

often concerned with sub jects having

substantial public interest.  TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientic

and technical material pertinent to

NASA's mission.

Specialized services that complement the STI

Program Oce's diverse oerings include

creating custom thesauri, building customized

databases, organizing and publishing

research results. . . even providing videos.

For more information about the NASA STI

Program Oce, see the following:

 Access the NASA STI Program Home

Page at http://www.sti.nasa.gov

 E-mail your question via the Internet to

help@sti.nasa.gov

 Fax your question to the NASA STI Help

Desk at (301) 621{0134

 Phone the NASA STI Help Desk at (301)

621{0390

 Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076{1320
NASA/TM{2001{210862

An Implicit Characteristic Based

Method for Electromagnetics

John H. Beggs

Langley Research Center, Hampton, Virginia

W. Roger Briley

Mississippi State University, Mississippi State, Mississippi

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23681{2199

May 2001
Available from:

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS)

7121 Standard Drive 5285 Port Royal Road

Hanover, MD 21076{1320 Springeld, VA 22161{2171

(301) 621{0390 (703) 605{6000
Abstract

An implicit characteristic-based approach for numerical solution of Maxwell’s time-dependent curl

equations in flux conservative form is introduced. This method combines a characteristic based finite difference spatial approximation with an implicit lower-upper approximate factorization (LU/AF) time integration scheme. This approach is advantageous for three-dimensional applications because the characteristic

differencing enables a two-factor approximate factorization that retains its unconditional stability in three

space dimensions, and it does not require solution of tridiagonal systems. Results are given both for a Fourier

analysis of stability, damping and dispersion properties, and for one-dimensional model problems involving propagation and scattering for free space and dielectric materials using both uniform and nonuniform

grids. The explicit FDTD algorithm is used as a convenient reference algorithm for comparison. The onedimensional results indicate that for low frequency problems on a highly resolved uniform or nonuniform

grid, this LU/AF algorithm can produce accurate solutions at Courant numbers significantly greater than

one, with a corresponding improvement in efficiency for simulating a given period of time. This approach

appears promising for development of LU/AF schemes for three dimensional applications.

1 Introduction

1.1 Background

In the regime of finite-difference solutions for Maxwell’s time dependent curl equations, explicit methods such as the Finite-Difference Time-Domain (FDTD) method [1]–[4], the Transmission Line Method (TLM), or (more recently) Finite-Volume Time-Domain (FVTD) methods have

been standard techniques for the past 20-30 years. These methods are relatively simple and efficient, and have proven very robust and adequate for many types of problems.

Since these explicit methods are conditionally stable, their maximum time step is limited by

a stability constraint that depends on the local grid spacing. The conditional stability restriction

can increase the cost of analyzing electromagnetic (EM) problems, especially when nonuniform

grids are needed, because a small local grid structure can require a time step smaller than would

otherwise be required to accurately resolve the physical wave propagation. On the other hand,

implicit algorithms can provide unconditional stability which allows the time step to be selected

for accuracy and resolution of the frequency content of the excitation sources, without a stability

restriction.

Implicit algorithms were first proposed for electromagnetics by Holland [5]. They have not received widespread use in electromagnetics, however, perhaps due to increased complexity and in

some cases the need to solve a linear system of equations. For example, centered spatial difference

approximations in conjunction with alternating direction implicit (ADI) schemes require the solution of (for example) tridiagonal systems for each coordinate direction. Furthermore, although this

type of scheme is unconditionally stable for the first-order wave equation in two dimensions, it is

known to be unconditionally unstable in three dimensions [6]. The usual Fourier stability analysis

does not account for nonperiodic boundary conditions, however, and it has been shown using a

matrix stability analysis that this three dimensional centered difference algorithm is conditionally

1
stable when used with upwind boundary conditions [7].

There has been recent work to develop characteristic based explicit and implicit algorithms

for CEM that are related to some solution algorithms from computational fluid dynamics. These

methods have been developed for curvilinear three dimensional grids, using both finite difference

and finite volume approximations. Shankar et al. [8], [9] developed explicit characteristic-based

finite-volume methods, using both flux-vector and flux-difference splittings for electromagnetics.

Shang, et. al. [10]–[19] have developed both implicit and explicit characteristic based methods.

Shang [10]–[11] developed an implicit characteristic based ADI method that leads to simple bidiagonal systems that are easily solved, rather than tridiagonal or pentadiagonal systems. This

characteristic-based scheme is unconditionally stable in two dimensions and conditionally stable

in three dimensions [11]. Shang and Fithen [18] also developed characteristic based finite difference and finite volume schemes using Runge-Kutta time integration, in curvilinear coordinates.

Gaitonde and Shang [20] and Turkel [4] have recently explored high order accurate compact differencing schemes that require tridiagonal solutions for spatial approximations, whether implicit

or explicit time integration is used. Turkel [4] has suggested a two dimensional implicit version

based on an ADI scheme. More recent work includes development of two and three-dimensional

implicit ADI FDTD schemes [21]–[24]. However, these ADI schemes still require solution of a

tridiagonal system of equations.

1.2 Present Work

A characteristic-based approach to spatial differencing has advantages in constructing implicit

algorithms, as well as in utilizing the natural method of incorporating boundary conditions for

hyperbolic systems. The ability of characteristic based methods to provide accurate nonreflecting

solutions at outer boundaries is well known and is discussed for electromagnetics by Shang [10]–

[11] and others.

The present method combines a characteristic based approach to spatial differencing with an

implicit lower-upper approximate factorization (LU/AF) time integration scheme that can be developed as a two-factor scheme that retains its unconditional stability in three space dimensions,

unlike three-factor ADI schemes. The LU/AF scheme also avoids the need to solve tridiagonal

implicit systems. A two-factor approximate factorization for three dimensions was first suggested

by Jameson and Turkel [25]. This type of factorization (two implicit passes in three dimensions)

has been combined with characteristic based differencing by Belk and Whitfield [26]. This general

approach has been further developed for CFD and is discussed, for example, in [27].

In the present report, an implicit characteristic based algorithm is developed for the time dependent Maxwell equations in one space dimension. This (one-dimensional) algorithm is investigated along two lines: A Fourier analysis is used to study both stability and dispersion/damping

properties of the algorithm. Although the Fourier analysis provides useful information and is

widely used to study algorithms, it is limited to periodic solutions on uniform grids. Complex

practical applications often involve nonperiodic solutions in finite regions with complex geometry and nonuniform grids. Grid nonuniformity is an important consideration because conditional

stability criteria typically link the maximum permissible time step to the local grid spacing rather

2
than to the underlying physical time and space ”sampling” requirements of the solution. Stable

implicit methods can benefit from choosing the time step to satisfy an accuracy requirement rather

than a stability constraint, with the consequence that fewer time steps are needed for a fixed simulation time. Accordingly, a one-dimensional model problem with nonuniform grid is developed

to explore the influence of grid nonuniformity on solution accuracy and efficiency. The explicit

FDTD algorithm is used as a convenient reference algorithm for comparison.

In the remainder of this report, Section 2 develops first and second order, characteristic based,

finite difference LU/AF algorithms for Maxwell’s equations in one dimension. The second order scheme includes the conduction current term. Section 3 gives a Fourier analysis of stability,

damping and dispersion properties of the one dimensional algorithms. Section 4 discusses the

characteristic based treatment of boundary conditions. Section 5 gives computed model problem

results related to both accuracy and efficiency for propagation and scattering in free space and dielectric materials, using both uniform and nonuniform grids. Finally, some conclusions are drawn

in Section 6.

2 One-Dimensional LU/AF Algorithm

This section details the theoretical background for development of both first and second-order

accurate LU/AF algorithms. The second order algorithm is used in all comparisons with the

FDTD method, and the first order algorithm is used at points adjacent to the outer boundaries of

the computational grid.

2.1 First-order algorithm

Maxwell’s equations for linear, homogeneous and lossless media in the one-dimensional case

(taking @=@y = @=@z = 0) are

@Ey

@t +

1



@Hz @x = 0 (1)

@Hz @t +

1



@Ey

@x = 0 (2)

These equations can be rewritten in flux conservative form as

@q

@t +

@ f

@x = 0 (3)

where q = [Ey ; Hz ]T

, f = [Hz= ; Ey=]Tand T denotes transpose. To develop the upwind

LU/AF algorithm, the flux conservative form of Maxwell’s equations in (3) is recast in the following form

@q

@t + A @q@x = 0 (4)

where A

is the Jacobian of f and is given by

A =

@ f

@q =

"

0 1=

1= 0 #

(5)

3
The matrix A

has eigenvalues 1;2 = 1=p corresponding to right and left propagating waves

with speeds c. The eigenvalue matrix is given by

 = " 1 0

0 2 #

(6)

The matrix A

can be obtained from the eigenvalue matrix via a similarity transformation given by

A = SS1(7)

where the matrix S

is composed of the eigenvectors of Aand is given by

S =

" p

= p=

1 1 #

; S1 =

1

2 " p

= 1 p

= 1 #

(8)

The eigenvalue matrix 

can be split into two parts, one each for the right and left propagating

waves and is given by

 = + + 

(9)

+ =

" 1 0

0 0 #

;   =

"

0 0

0 2 #

(10)

Substituting (9) into (7) gives

A = S

+ +  S1 = A+ + A

(11)

A = S S1

(12)

The flux vector 

f is then split into two parts given by



f = f + + f (13)



f  = Aq (14)

This flux-vector-splitting method is similar to that developed by Steger and Warming [28] for the

Euler equations governing inviscid fluid flow. To construct the LU/AF algorithm, time and space

are discretized to tn = nt, xi = ix and the term @ 

f =@x is approximated at t = tn + t by

@ f

@x

  



f +

i  

f +

i1

x

+



f 

i+1  

f 

i

x !n+1

+ (1  ) 



f +

i  

f +

i1

x

+



f 

i+1  

f 

i

x !n

(15)

Here,  (0    1) is a parameter for evaluating the spatial approximation at tn + t. This will

give an explicit scheme for  = 0 and an implicit scheme for  > 0. The finite difference equation

for (3) is given by

q

n+1

i  qn

i

t +  



f +

i  

f +

i1

x

+



f 

i+1  

f 

i

x !n+1

+ (1  ) 



f +

i  

f +

i1

x

+



f 

i+1  

f 

i

x !n = 0 (16)

Making the substitution qn

i  q

n+1

i  qn

i

, equation (16) can be rearranged for a uniform grid to

give

h 

I=t + 



1i A+

 + +

1i A

i qn

i = Rn1i

(17)

4
where 

I is the 2 by 2 identity matrix, and



1i() 

()i  ()i1

x

(18)

+

1i() 

()i+1  ()i

x

(19)

are the first-order backward and forward difference operators, respectively. The term Rn

1i

is called

the “residual” and is defined by

Rn

1i  

1i



f +;n

i + +1i



f ;n

i

(20)

Note that equation (17) is an unfactored, implicit scheme of O

x; (  1=2) t; t2

accuracy for

the numerical solution of Maxwell’s equations. The solution of (17) involves solving a block tridiagonal system of equations. The solution of tridiagonal systems can be avoided by an approximate

factorization (AF) of the left side of (17) which results in the following LU/AF algorithm:

h 

I=t + 

1i A+

i q



i = Rn1i

(21)

h 

I=t + +

1i A

i qn

i = 

I=tq



i

(22)

where q



i

is an intermediate solution vector. This LU/AF algorithm is closely related to those

developed previously for Computational Fluid Dynamics (e.g. [27]). It can be shown that these

equations are a consistent and unconditionally stable approximation to (3), but the details of that

analysis are omitted here for the sake of brevity. Note that the solution involves a two-step procedure. In solving (21), a sweep through the grid in the forward direction results in a block lower

bidiagonal matrix and a sweep through the grid in the backward direction results in a block upper

bidiagonal matrix for (22). These equations are then solved by forward and backward substitution, respectively. Thus, a costly inversion of a block tridiagonal matrix is avoided, and each step

in the solution is effectively explicit. Substituting equation (22) into (21) and expanding gives,

h 

I=t + 



1i A+

 + +

1i A

 + 2t 

1i A+

 +

1i A

i qn

i = Rn1i

(23)

Note that equation (17) is recovered except for the third term in brackets, which represents the factorization error. In this particular example, A+

and Aare constant matrices such that A+A  0,

and the factorization error is zero. However, the factorization error is nonzero in the following

example with conduction currents, and more generally in two and three dimensional implementations. This error can be reduced or eliminated by iterative error reduction. The present onedimensional model problem results do not include factorization error except for the case including

conduction currents.

2.2 Second-order algorithm

To develop the O

t2

; x2LU/AF algorithm, the flux conservative form of (4) is rewritten

to include the electric and magnetic conduction currents as

@q

@t + A @q@x = P q (24)

5
where P

is given by

P =

" = 0

0 = #

(25)

and A

is given by equation (5). The time derivative in (24) is approximated by a -weighted,

O(t2) difference equation given by

@q

@t



(2 + 1)qn+1

i  4qn

i + (2  1)qn1

i

2t

(26)

The spatial derivative is again replaced by a -weighted average between time level n + 1 and n

as in (15). This can be rewritten using the operator notation as

@ f

@x

 



2i



f +

i + +2i



f 

i n+1 + (1  ) 2i



f +

i + +2i



f 

i n

(27)

The 

2i

operators are now O(x2) difference operators to be defined shortly. The parameter  can

be used to construct a series of explicit and implicit schemes. For example, if  = 0, this results

in a leapfrog scheme;  = 0:5 results in a Crank-Nicolson scheme and  = 1 results in an Euler

implicit scheme. Using (26) and (27), the finite-difference equation for (24) is

(2 + 1)qn

i  (2  1)qn1

i

2t + 



2i



f +

i + +2i



f 

i n+1 + (1  ) 2i



f +

i + +2i



f 

i n =

Pq

n+1

i  (1  )P

q

n

i

(28)

Using the same flux vector splitting as in (14), this can be rearranged as

h 

I=(2t) + P +  

2iA+

()++

2iA

()i qn

i = Rn2i

(29)

The residual, Rn

2i

, is now defined by

Rn

2i   

Pq

n

i + 

2i



f +;n

i + +2i



f ;n

i 

2  1

2t

qn1

i 

(30)

where   1=(2 + 1). The difference operators in equation (30) are replaced by O(x2

) backward

and forward upwind difference operators on three-point one-sided stencils defined by



2i()  (3()i  4()i1 + ()i2) =(2x) (31)

+

2i()  (3()i + 4()i+1  ()i+2) =(2x) (32)

Substituting (14) into (30), the residual is given by

Rn

2i =  

Pq

n

i + 

2iA+

q

n

i + +

2iA

q

n

i 

2  1

2t

qn1

i 

(33)

Equation (29) is an O(t2

; x2), unfactored, upwind scheme for electromagnetics. The LU/AF

scheme is defined by factoring the left side of (29) into two operators, each designed for a forward

and backward grid sweep as in the first order implementation. The LU/AF scheme is then given

by

h 

I=(2t) + P + 

2iA+

()i q



i = Rn2i

(34)

h 

I=(2t) + P + +

2iA

()i qn

i =



I =(2t) + P

 q



i

(35)

6
3 Fourier Analysis

A Fourier analysis shows that both the first and second-order upwind LU/AF algorithms are

unconditionally stable for   1=2, and that they contain both numerical dissipation (or damping)

and dispersion. The dissipation is present due to even order spatial derivatives in the truncation

error which are a result of the upwind approximation. These schemes are consistent approximations with O((  1=2)t; t2

; x) truncation error for the first-order algorithm and O(t2; x2)

for the second-order algorithm. A Fourier analysis is now given for the first-order LU/AF algorithm; the second-order algorithm analysis follows a similar development.

Assuming an equally spaced mesh with periodic boundary conditions, a trial solution of the

form

q

n

i = q0e

j(ni)

(36)

is substituted into equations (21) and (22) to obtain:

T qn

i = U

q

n

i

(37)

The matrices T

and Uare given by

T =

 

I +



c 

1  e

j A+

  

I +



c 

ej  1 A



(38)

U = 





c 

1  e

j A+ +



c 

ej  1 A

(39)

where c = 1=p

,  is the Courant number defined by  = ct=x and  =  . Again making

the substitution qn

i  q

n+1

i  qn

i

, equation (37) can be rearranged to give

q

n+1

i = G

q

n

i

(40)

where the amplification matrix G

is defined by

G = I + T1U(41)

The stability of this scheme is governed by the spectral radius of the amplification matrix G

, and

its eigenvalues, G1 and G2. Closed form solutions were obtained for G1 and G2, but they are

omitted here for the sake of brevity. A rigorous numerical stability analysis was also performed

on the second-order LU/AF algorithm using the MAP LE 6TM software package. Checking the

magnitude of the eigenvalues G1 and G2 for all values of ,  and  showed that both the first

and second-order schemes are unconditionally stable for   1=2. The stability analysis in this

particular case is identical to the unfactored implicit scheme given in equation (17), because the

factorization error is zero for P = 0.

For the dispersion analysis, the term

q

n+1

i = qn

i e

j

(42)

is substituted in (40) to give

G  e

j 

Iq

n

i = 0 (43)

7
Equation (43) is then solved for 

, the numerical wavenumber. The quantity =m f=g gives

the numerical dispersion, or phase velocity error. Extensive numerical experiments were performed to analyze the dispersive properties of these algorithms. We let  = 2=N, where N is the

grid resolution in cells/. Figure 1 shows the numerical dispersion versus grid resolution with

 = 3=4 for the first and second-order LU/AF algorithms as compared to the FDTD algorithm,

which is O(t2

; x2). The Courant number of  = 3=4 was chosen for these methods simply as

a representative Courant number which could be expected for a problem involving a nonuniform

grid. For a nonuniform mesh with a mesh stretch ratio of Ms : 1 where Ms  xmax=xmin,

the Courant number for FDTD varies pointwise in the grid over the range 1=Ms  i  1, where

i  (ct)=xi . The numerical experiments showed that the first-order LU/AF algorithm has

Yee

First Order

Second Order

–20

0

20

40

60

80

100% error in dispersion

10 20 30 40 50 60 70 80

Grid Resolution (cells/lambda)

Figure 1: Numerical dispersion of FDTD, first and second-order LU/AF algorithms versus grid

resolution for  = 1=2 and  = 3=4.

lower dispersion errors than the second-order LU/AF method for   1. This is not particularly

troublesome, because for   1, explicit schemes are generally more efficient and are preferred.

However, for  > 1, explicit schemes cannot be used and the second-order LU/AF algorithm does

have lower dispersion error than the first-order LU/AF method. For large values of , the dispersive error for the first and second-order LU/AF schemes begin to converge as shown in Figure

2. The results also showed that the lowest numerical dispersion is obtained when  = 1=2 as

shown in Figure 3. Since the LU/AF method uses windward differencing, numerical dissipation

(or damping) is present in the solution. The normalized dissipation is obtained as <e f=g. Figures 4-6 show numerical dissipation for the first and second-order LU/AF schemes versus grid

resolution,  and , respectively. Note that the second-order algorithm has much lower dissipation

than the first order method.

8
First Order

Second Order

0

10

20

30

40

50

60% error in dispersion

2 4 6 8 10

nu

Figure 2: Numerical dispersion of first and second-order LU/AF algorithms versus Courant number  at  = 1=2 and N = 10 cells/ resolution.

First Order

Second Order

10

20

30

40% error in dispersion

0.5 0.6 0.7 0.8 0.9 1

beta

Figure 3: Numerical dispersion of first and second-order LU/AF algorithms versus  at  = 2 and

N = 10 cells/ resolution.

9
First Order

Second Order

0

0.2

0.4

0.6

0.8

Normalized dissipation

10 20 30 40 50 60 70 80

Grid Resolution (cells/lambda)

Figure 4: Normalized dissipation of first and second-order LU/AF algorithms versus grid resolution at  = 0:5 and  = 0:75.

First Order

Second Order

0

0.05

0.1

0.15

0.2

0.25

Normalized dissipation

2 4 6 8 10

nu

Figure 5: Normalized dissipation of first and second-order LU/AF algorithms versus Courant

number  at  = 1=2 and N = 10 cells/ resolution.

10
First Order

Second Order

0.1

0.2

0.3

0.4

0.5

0.6

Normalized dissipation

0.6 0.7 0.8 0.9 1

beta

Figure 6: Normalized dissipation of first and second-order LU/AF algorithms versus  at  = 2

and N = 10 cells/ resolution.

4 Boundary Conditions

4.1 Outer Radiation Boundary Condition

The characteristic based LU/AF method requires no extraneous boundary condition such as

the Liao absorbing boundary condition [29] or the PML [30]. The FDTD method uses a spatial

central difference operator, which for a wave propagating from left to right, eventually requires a

grid point outside the domain. This requirement introduces an additional equation (i.e. boundary

condition) to solve the system and introduces information into the solution that is not required

by Maxwell’s equations. Using an upwind characteristic based approach, the interior point algorithm calculates the left-going characteristic at the left boundary (i.e. at i = 0) and the right-going

characteristic at the right boundary (i.e. at i = imax). Therefore, the only additional information

required is information about waves that are entering the domain. Waves exiting the domain are

naturally handled by the interior point algorithm. Therefore, the characteristic boundary conditions are implemented as follows: at grid point i = 0, equation (22) is used along with a specification of the incoming, right-going flux, f +

0 . At grid point i = imax, equation (21) is used along with

a specification of the incoming, left-going flux, f 

imax. Therefore, the only additional information

introduced at the boundary is nothing more than what is required by the physical system. In multidimensional problems, the local coordinates at the outer boundaries are rotated to align with the

direction of wave propagation defined by E H

. The characteristic equations are developed along

this direction and are appropriately applied at the boundaries. This procedure was discussed and

outlined by Shang [11].

11
4.2 Dielectric Surface Boundary Condition

Since the LU/AF scheme follows the direction of information propagation (i.e. the characteristic), at a material interface, the slope of the characteristic curve (i.e. the speed of light in the

material) changes. Therefore, for the LU/AF method to be widely applicable, a careful treatment

of material interfaces is required. With a material boundary in place, the right-going and left-going

characteristics see a change in characteristic speeds, and therefore, a material interface condition

needs to be implemented to correctly model the physics. To implement this feature, consider the

one-dimensional grid shown in Figure 7 where a dielectric boundary has been inserted at grid

point ib. For the second order LU/AF method, at an arbitrary grid point i, the windward differy

1

i -3 i -1 b i -2 b b

i

1 ε

b

μ1 σ

x

f

-

ε 2 μ 2 σ2

i +1 b i +2 b

i +3 b

f2

+

1

Figure 7: One-dimensional FDTD grid showing a dielectric boundary at grid point ib .

encing uses flux components f +

i2, f +i1, f +i

, f 

i

, f 

i+1, f i+2. This creates a difficulty at a dielectric

interface since the algorithm assumes the material is homogeneous. We know at grid point ib that

the physical boundary conditions require that

Etan1 = Etan2 (44)

Htan1 = Htan2 (45)

which requires that

Ey1 = Ey2 (46)

Hz1 = Hz2 (47)

Imposing these boundary conditions at ib then makes the solution vector q constant between material 1 and 2. Therefore, the same solution vector q can be used for both materials at grid point

12
ib. In general, in regions 1 and 2, we require flux components f +

1 , f 1 , f +2 and f 2 , which are restricted to the respective materials and are dependent on the material properties. However, we

define additional flux components f 

1 located at grid points ib, ib + 1 and ib + 2 in material 2 and

flux components f +

2 at grid points ib and ib  1 in material 1. By defining these components, we

can then apply the second order windward difference formulas at grid points ib  2, ib  1, ib, ib + 1 and ib + 2 to simulate a homogeneous material region. For example, on the forward sweep,

at grid point ib, we use the components f +

1 (ib  2), f +1 (ib  1), f +1 (ib), f 1 (ib), f 1 (ib + 1), f 1 (ib + 2)

to update the solution vector q using the material properties of region 1. Similar analogies can be

drawn at grid points ib  2, ib  1, ib + 1 and ib + 2. These additional fluxes in each region in the

vicinity of the dielectric boundary are computed from (14) based upon the material properties in

each region, and they bridge the solution between regions 1 and 2.

5 Model Problem Results

The model problem results presented in this section demonstrate that the one dimensional implicit characteristic-based LU/AF algorithm can produce accurate results on a nonuniform grid

for time steps significantly larger than the maximum permissible for a typical conditionally stable

scheme. Although the operation count for the implicit scheme is higher, as the mesh variation is

increased, fewer time steps are needed for an accurate solution over a fixed simulation time, and

this savings more than offsets the additional operations (beyond a certain mesh variation).

The second-order LU/AF algorithm was tested by implementing equations (34) and (35) for

interior grid points away from absorbing boundaries and equations (21), (22) at grid points next

to the absorbing boundaries. Characteristic-based boundary conditions were used to terminate

the computational domain and the incoming flux (f 

imax) at the right boundary was set to zero.

The code was initialized by writing a time snapshot of a propagating pulse in the grid. Several

different types of problems were analyzed with the LU/AF algorithm and compared with the

FDTD method in an attempt to assess the characteristics of the LU/AF algorithm as applied to a

one-dimensional problem with a non-uniform grid. These are outlined in the following sections.

5.1 Propagation

To assess the dispersion and dissipation characteristics of the LU/AF algorithm, several propagation problems were analyzed using free space and dielectric materials and were compared with

FDTD. A Gaussian pulse of the form

Eyi = e



(tt0)  2

(48)

was used as an excitation source for both propagation and scattering problems. For propagation,

the code was changed to use periodic boundary conditions to enable simulation of propagation

over large distances compared to the shortest wavelength contained in the frequency spectrum

of the pulse. To illustrate the potential benefit of the LU/AF algorithm over conventional explicit

schemes, the main parameters of interest are the time resolution (i.e. number of time steps/period)

and the grid resolution (i.e. number of cells/wavelength) of the highest frequency of interest in

13
any given problem. These parameters will be designated as Nt and Nx, respectively. The other

parameter of interest is the mesh stretch ratio, Ms. The highest frequency of interest, fmax, is

calculated based upon Nx and the largest cell size, and the time step for the implicit scheme is

calculated based upon fmax and the desired time resolution, Nt . For the FDTD method, the time

step was calculated based upon the Courant stability condition using the smallest grid size.

The first propagation problem involves propagation in free space on a uniform mesh. The

time and grid resolutions are Nt = 30 and Nx = 30, the cell size was 1 cm, the time step was 33.3

ps and  = 1. The maximum frequency of interest was 1 GHz and the problem space was 2000

cells. The pulse was allowed to propagate a distance equal to 100 wavelengths of the minimum

wavelength, which corresponds to a distance of about 30 meters. Figure 8 shows the electric field

versus distance obtained after 3000 time steps computed using the exact solution, FDTD and the

LU/AF algorithm. Note the LU/AF scheme has slightly less accurate results due to the larger dis-

-0.2

0

0.2

0.4

0.6

0.8

1

30 31 32 33 34 35

Electric field (V/m)

x (meters)

Exact

FDTD

LU/AF

Figure 8: Electric field versus distance for free space propagation on a uniform mesh with  = 1=2

and  = 3=4.

persion error; but the general shape of the pulse is still acceptable. Again, this is not troublesome,

because for   1, explicit schemes are more efficient and are preferred. Further numerical simulations confirmed that the numerical dispersion was lower and the results were more accurate

for the second-order LU/AF method for the case when  = 2 as shown in Figure 9. Additional

simulations confirmed that the numerical dispersion error increases with Courant number  as

with the ADI FDTD scheme [24].

To demonstrate the LU/AF scheme on nonuniform meshes, the same problem outlined above

was simulated using mesh stretch ratios of Ms=10:1 and Ms=100:1 and the mesh was periodic every 10 cells. Figure 10 illustrates an expanded section of the nonuniform mesh for a mesh stretch

ratio of Ms =10:1. The smallest grid cells force the FDTD algorithm to take a time step 1/10th of

the previous case, and it must run 10 times as many time steps to propagate the pulse the same

14
-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

30 31 32 33 34 35

Electric field (V/m)

x (meters)

Exact

LU/AF

Figure 9: Electric field versus distance for free space propagation on a uniform mesh with  = 1=2

and  = 2.

0 0.05 0.1 0.15 0.2 0.25 0.3

x (meters)

Figure 10: A section of a nonuniform mesh with a mesh stretch ratio of 10:1 and largest cell size of

1 cm.

15
distance. However, the LU/AF algorithm can take the same time step as before, because the time

resolution, Nt = 30, has been preserved and the time step is set solely upon accuracy requirements, not on the Courant stability condition. The maximum cell size remains constant, and since

Nx is also the same, the maximum frequency of interest remains constant at 1 GHz. The maximum

grid coordinate changes from 20.0 m in the uniform mesh case to 8.17 m for the 10:1 nonuniform

mesh. The Courant numbers in the grid for FDTD are in the range 0:1    1 and the Courant

numbers for the LU/AF algorithm are in the range 1    10. Figure 11 shows the electric

field versus distance obtained with the exact solution and after 30,000 time steps computed using

FDTD and after 3,000 time steps using the LU/AF algorithm. Note the agreement is very good

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

30 31 32 33 34 35

Electric field (V/m)

x (meters)

Exact

FDTD

LU/AF

Figure 11: Electric field versus distance for free space propagation on a nonuniform mesh with

 = 1=2 and  = 1.

for both methods in this case. In fact, Figure 12 shows the error in electric field versus distance

for both methods and the errors are almost equal in magnitude. The LU/AF method still required

about 15 seconds, but the FDTD method required about 22 seconds. Note the improved performance of the LU/AF algorithm away from the area of the main pulse. Since the mesh stretch ratio

was 10 for this problem, the FDTD method is forced by the Courant stability condition to take 10

times the number of time steps to simulate the same amount of physical time. It appears that the

crossover point in efficiency between the LU/AF and FDTD methods is approximately a mesh

stretch ratio of 10:1. This is in contrast to recent results for the ADI FDTD scheme [23] where the

results suggested a crossover point of 30:1 mesh ratio.

For a mesh stretch ratio of Ms =100:1, the Courant number for FDTD was in the range 0:01 

  1 and for LU/AF was in the range 1    100. Execution time for FDTD was 230 seconds

for 300,000 time steps. Execution time for the LU/AF method again required 15 seconds for 3,000

times steps and the error in electric field for this simulation is shown in Figure 13. Note again

that the LU/AF result appears “cleaner” than the FDTD result with less numerical noise. Based

upon these results, we conclude that the LU/AF method is more efficient for either a.) nonuni16
-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

30 31 32 33 34 35

% Error in electric field (V/m)

x (meters)

FDTD

LU/AF

Figure 12: Error in electric field versus distance for free space propagation on a nonuniform mesh

with a mesh stretch ratio of 10:1 and  = 1=2,  = 1 (for FDTD).

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

30 31 32 33 34 35

% Error in electric field (V/m)

x (meters)

FDTD

LU/AF

Figure 13: Error in electric field versus distance for free space propagation on a nonuniform mesh

with a mesh stretch ratio of 100:1 and  = 1=2,  = 1 (for FDTD).

17
form meshes with large mesh stretch ratios or b.) a uniform mesh with a low frequency incident

excitation, where the pulse is highly resolved spatially over its entire frequency spectrum. Further

computational experience on additional free space propagation problems on nonuniform grids reveals that this sweeping method fails for large t (i.e. large Courant number). This can occur for

two reasons, even though the Fourier stability analysis indicates unconditional stability. First, it is

inapropriate to use a time step larger than one which propagates a wave over a distance comparable to the entire solution domain in a single time step, and this can cause instability. Secondly, the

spatial sweeping process associated with each factor of the algorithm can itself become unstable

if it loses diagonal dominance, as observed by Jameson and Turkel [25]. Although a first-order

spatial upwind scheme guarantees diagonal dominance, there is a loss of diagonal dominance for

higher-order upwind schemes as the time step becomes asymptotically large. The upwind scheme

significantly increases the diagonal term compared with a centered scheme, however, and in applying the LU/AF scheme, this problem can be avoided by exercising care in selecting the time

step.

5.2 Lossy Dielectric Materials

To illustrate propagation in dielectric materials, the problem space was filled with a lossy dielectric material with parameters  = 40,  = 0 and  = 0:02. The pulse was allowed to propagate on a uniform mesh for approximately 35 wavelengths at the minimum wavelength. Figure

14 shows the results of FDTD versus LU/AF with excellent agreement between the two methods.

A similar level of agreement is found with the same problem on a nonuniform grid.

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

10 15 20 25 30 35

Electric field (V/m)

x (meters)

FDTD

LU/AF

Figure 14: Electric field versus distance for propagation in a lossy dielectric on a uniform mesh

with  = 1=2 and  = 1.

To illustrate the use of the LU/AF scheme for scattering problems, reflection from a lossy

dielectric half-space was considered as the benchmark problem. The dielectric interface scheme

18
of Section 4.2 was implemented and the results are again compared with the FDTD method. To

solve this problem, the problem space is 2000 cells and it is filled with a lossy dielectric material

from cell numbers 751-2000. The total electric field is recorded versus time at cell number 750. The

incident field is obtained by running the code with free space only and recording the field versus

time at the same location. The incident field is subtracted from the total field to give the scattered

field. A point source located at the left boundary of the grid was used as the excitation source

and the dielectric half-space material parameters are  = 40,  = 0 and  = 0:2 S/m. Figure

15 shows the time-domain scattered field for this case, and the agreement is excellent. Figure 16

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0 10 20 30 40 50 60

Electric field (V/m)

Time (ns)

FDTD

LU/AF

Figure 15: Scattered electric field versus time for scattering from a lossy dielectric half-space on a

uniform mesh with  = 1=2 and  = 1.

shows the reflection coefficient results and again the agreement is excellent. The same problem

was run on the uniform mesh with  = 2, and on a nonuniform mesh with a mesh stretch ratio of

10:1 and  = 0:2 S/m. Figures 17 and 18 show the respective results.

Note that the LU/AF method produces more inaccurate results for moderate conductivity values. This is due to the factorization error terms. As  is increased, the factorization errors become

larger. In principle, this factorization error can be reduced throught iterative error reduction [31].

This has been demonstrated for CFD problems [27], but it is not explored in the present work.

For perfect conductors, we simply set q = 0. Figure 18 shows that the LU/AF method is superior

on nonuniform grids, even for problems involving lossy dielectrics. These results indicate that

the physical dielectric boundary conditions have been properly implemented, thus making the

LU/AF method applicable to a wider class of problems.

6 Conclusions

This report has introduced a LU/AF implicit finite-difference method for computational electromagnetics. A second-order accurate, LU/AF, characteristic based algorithm for electromag19
0.5

0.55

0.6

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1

0 0.1 0.2 0.3 0.4 0.5

Reflection coefficient magnitude

Frequency (GHz)

Exact

FDTD

LU/AF

Figure 16: Reflection coefficient magnitude versus frequency for scattering from a lossy dielectric

half-space on a uniform mesh with  = 1=2 and  = 1.

0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

0.96

0.98

1

0 0.1 0.2 0.3 0.4 0.5

Reflection coefficient magnitude

Frequency (GHz)

Exact

FDTD

LU/AF

Figure 17: Reflection coefficient magnitude versus frequency for scattering from a lossy dielectric

half-space on a uniform mesh with  = 1=2,  = 1 and  = 2 S/m.

20
0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.1 0.2 0.3 0.4 0.5

Reflection coefficient magnitude

Frequency (GHz)

Exact

FDTD

LU/AF

Figure 18: Reflection coefficient magnitude versus frequency for scattering from a lossy dielectric

half-space on a uniform mesh with  = 1=2,  = 1 and  = 0:2 S/m.

netics has been implemented and tested on one-dimensional model problems for uniform and

nonuniform grids. The one-dimensional model problem results for the characteristic based implicit scheme demonstrate that:

1. Accurate solutions for wave propagation can be obtained using Courant number significantly greater than one on nonuniform grids.

2. Stability of the implicit scheme allows a given amount of time to be covered in fewer time

steps than a conditionally stable scheme, and this can significantly improve efficiency when

nonuniform grids are needed. This is also extremely beneficial for low frequency excitation

sources where geometric details may force a highly oversampled mesh (uniform or nonuniform) in certain regions.

3. The method works well with lossy dielectric materials and perfect conductors.

4. The factorization error increases as  increases, but it can be reduced in principle through

iterative error reduction.

5. The lowest dispersion and damping errors occur when  = 1=2.

The present results demonstrate potential advantages in a one dimensional context, and this

approach appears promising for development of stable, accurate and efficient implicit LU/AF

schemes for complex two and three dimensional applications. Extensions to two and three-dimensional

applications have been outlined [31], but details of this work will the subject of future reports and

articles.

21
References

[1] K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp.

302–307, Mar. 1966.

[2] K. S. Kunz and R. J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics,

CRC Press, Boca Raton, FL, 1993.

[3] A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech

House, Boston, MA, 1995.

[4] A. Taflove, Ed., Advances in Computational Electrodynamics: The Finite-Difference Time-Domain

Method, Artech House, Boston, MA, 1998.

[5] R. Holland, L. Simpson, and K. Kunz, “Finite-difference analysis of EMP coupling to lossy

dielectric structures,” IEEE Transactions on Electromagnetic Compatibility, vol. EMC-22, no. 3,

pp. 203–209, Aug. 1980.

[6] R. F. Warming and R. M. Beam, “An extension of a-stability to alternating direction implicit

methods,” BIT, vol. 19, pp. 395–417, 1979.

[7] R. C. Buggeln W. R. Briley and H. McDonald, “Solution of the three-dimensional navier

stokes equations for a steady laminar horseshoe vortex flow,” in Proc. 7th Computational Fluid

Dynamics Conference, Cincinnati, OH, July 1985, AIAA Paper 85-1520-CP.

[8] W. F. Hall V. Shankar and A. H. Mohammadian, “A time-domain differential solver for electromagnetic problems,” Proc. IEEE, vol. 77, no. 5, pp. 709–721, May 1989.

[9] W. F. Hall V. Shankar and A. H. Mohammadian, “A time-domain, finite-volume treatment

for the Maxwell equations,” Electromagnetics, vol. 10, pp. 127, 1990.

[10] J. S. Shang, “Characteristic based methods for the time-domain Maxwell equations,” in AIAA

29th Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1991, vol. AIAA 91-0606.

[11] J. S. Shang, “A characteristic-based algorithm for solving 3-d time-domain Maxwell equations,” in AIAA 30th Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1992, vol. AIAA

92-0452.

[12] J. S. Shang, “A fractional-step method for solving 3-d time-domain Maxwell equations,” in

AIAA 31st Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1993, vol. AIAA 93-0461.

[13] J. S. Shang and D. Gaitonde, “Characteristic-based, time-dependent Maxwell equations

solvers on a general curvilinear frame,” in AIAA 24th Plasmadynamics & Lasers Conference,

Orlando, FL, July 1993, vol. AIAA 93-3178.

[14] K. C. Hill J. S. Shang and D. Calahan, “Performance of a characteristic-based, 3-d timedomain Maxwell equations solvers on a massively parallel computer,” in AIAA 24th Plasmadynamics & Lasers Conference, Orlando, FL, July 1993, vol. AIAA 93-3179.

22
[15] J. S. Shang and R. M. Fithen, “A comparative study of numerical algorithms for computational electromagnetics,” in AIAA 25th Plasmadynamics & Lasers Conference, Colorado Springs,

CO, June 1994, vol. AIAA 94-2410.

[16] J. S. Shang and D. Gaitonde, “Characteristic-based, time-dependent Maxwell equation

solvers on a general curvilinear frame,” AIAA Journal, vol. 33, no. 3, pp. 491–498, March

1995.

[17] J. S. Shang, “A fractional-step method for solving 3d, time-domain Maxwell equations,”

Journal of Comp. Phys., vol. 118, pp. 109–119, 1995.

[18] J. S. Shang and R. M. Fithen, “A comparative study of characteristic-based algorithms for the

Maxwell equations,” Journal of Comp. Phys., vol. 125, pp. 378–394, 1996.

[19] D. C. Blake and J. S. Shang, “A procedure for rapid prediction of electromagnetic scattering

from complex objects,” in AIAA 29th Plasmadynamics & Lasers Conference, Albuquerque, NM,

June 1998, vol. AIAA 98-2925.

[20] D. Gaitonde and J. S. Shang, “High-order finite-volume schemes in wave propagation phenomena,” in AIAA 27th Plasmadynamics & Lasers Conference, New Orleans, LA, June 1996, vol.

AIAA 96-2335.

[21] F. Zheng, Z. Chen, and J. Zhang, “A finite-difference time-domain method without the

Courant stability conditions,” IEEE Microwave Guided Wave Letters, vol. 9, no. 11, pp. 441–

443, Nov. 1999.

[22] F. Zheng, Z. Chen, and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Transactions on Microwave Theory

and Techniques, vol. 48, no. 9, pp. 1550–1558, Sept. 2000.

[23] T. Namiki, “3-D ADI-FDTD method—Unconditionally stable time-domain algorithm for

solving full vector Maxwell’s equations,” IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 10, pp. 1743–1748, Oct. 2000.

[24] A. Taflove and S. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain

Method, 2 ed., Artech House, Boston, MA, 2000.

[25] A. Jameson and E. Turkel, “Implicit schemes and approximate factorizations,” Math. Comp.,

vol. 37, pp. 385–397, 1981.

[26] D. M. Belk and D. L. Whitfield, “Unsteady three-dimesional euler solutions on blocked grids

using an implicit two-pass algorithm,” January 1987, vol. AIAA Paper No. 87-0450.

[27] W. R. Briley, S. S. Neerarambam, and D. L. Whitfield, “Implicit Lower-Upper/ApproximateFactorization algorithms for incompressible flows,” Journal of Comp. Phys., vol. 128, pp. 32–42,

1996.

[28] J. L. Steger and R. F. Warming, “Flux vector splitting of the invisid gas dynamic equations

with applications to finite difference methods,” Journal of Comp. Phys., vol. 4, no. 2, pp. 263–

293, 1981.

23
[29] Z. P. Liao, H. L. Wong, B.-P. Yang, and Y.-F. Yuan, “A transmitting boundary for transient

wave analysis,” Sci. Sin., Ser. A, vol. 27, no. 10, pp. 1063–1076, Oct. 1984.

[30] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,”

Journal of Computational Physics, vol. 114, no. 1, pp. 185–200, 1994.

[31] J. H. Beggs and W. R. Briley, “An implicit characteristic based method for computational

electromagnetics,” Tech. Rep. MSSU-EIRS-ERC-98-11, Miss. State Univ., August 1998.

24
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704–0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis

Highway, Suite 1204, Arlington, VA 22202–4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704–0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE

May 2001

3. REPORT TYPE AND DATES COVERED

Technical Memorandum

4. TITLE AND SUBTITLE

An Implicit Characteristic Based Method for Electromagnetics

5. FUNDING NUMBERS

706–31–41–01

6. AUTHOR(S)

John H. Beggs, W. Roger Briley

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681–2199

8. PERFORMING ORGANIZATION

REPORT NUMBER

L–18051

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546–0001

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA/TM–2001–210862

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Category 33 Distribution: Standard

Availability: NASA CASI (301) 621–0390

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

An implicit characteristic-based approach for numerical solution of Maxwell’s time-dependent curl equations in flux

conservative form is introduced. This method combines a characteristic based finite difference spatial

approximation with an implicit lower-upper approximate factorization (LU/AF) time integration scheme. This

approach is advantageous for three-dimensional applications because the characteristic differencing enables a

two-factor approximate factorization that retains its unconditional stability in three space dimensions, and it does

not require solution of tridiagonal systems. Results are given both for a Fourier analysis of stability, damping and

dispersion properties, and for one-dimensional model problems involving propagation and scattering for free space

and dielectric materials using both uniform and nonuniform grids. The explicit FDTD algorithm is used as a

convenient reference algorithm for comparison. The one-dimensional results indicate that for low frequency

problems on a highly resolved uniform or nonuniform grid, this LU/AF algorithm can produce accurate solutions at

Courant numbers significantly greater than one, with a corresponding improvement in efficiency for simulating a

given period of time. This approach appears promising for development of dispersion optimized LU/AF schemes

for three dimensional applications.

14. SUBJECT TERMS

computational electromagnetics, FDTD methods

15. NUMBER OF PAGES

29

16. PRICE CODE

A03

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

19. SECURITY CLASSIFICATION

OF ABSTRACT

Unclassified

20. LIMITATION OF ABSTRACT