An Efficient Positive Definite Method for the Numerical Solution of the Advection Equation

Working Paper

AN EFFICIENT P o s r r m DEFINITE METHOD FOR THE NUMERICAL SOLUTION OF THE ADVECTION EQUATION

J e r z y B a r t n i c k i

August 1986 W P - 8 6 - 3 5

International Institute for Applied Systems Analysis

A-2361 Laxenburg, Austria

NOT FOR QUOTATION WITHOUT PERMISSION OF THE AUTHOR

AN EFFICIENT POSITIVE DEFIHITE METHOD FOR THE NUMERICAL SOLUTION OF THE ADVECTION EQUATION

J e r z y B a r t n i c k i

August 1986 W P - 8 6 - 3 5

Working Papers are interim r e p o r t s on work of t h e I n t e r n a t i o n a l I n s t i t u t e f o r Applied Systems Analysis a n d h a v e r e c e i v e d only lim- i t e d review. Views o r opinions e x p r e s s e d h e r e i n d o n o t neces- s a r i l y r e p r e s e n t t h o s e of t h e I n s t i t u t e o r of i t s National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS 2361 L a x e n b u r g , Austria

PREFACE

Under a c o l l a b o r a t i v e a g r e e m e n t with t h e I n s t i t u t e f o r Meteorology a n d Water Management, W a r s a w (Poland), J e r z y Bartnicki works with IIASA's Acid Rain P r o j e c t on modeling a t m o s p h e r i c t r a n s p o r t a n d deposition of pol- lutants. This p a p e r r e p o r t s o n a new method f o r solution of t h e advection equation of a n a i r pollution t r a n s p o r t model. This new method t u r n e d o u t t o b e a useful tool in t h e u n c e r t a i n t y analysis of t h e RAINS (Regional Acidifi- cation INformation a n d Simulation) model c a r r i e d o u t in t h e p r o j e c t .

Leen Hordijk

L e a d e r , Acid Rain P r o j e c t

I would like t o thank Krzysztof Abert from t h e Institute of Meteorology and Water Management in W a r s a w f o r many valuable discussions about t h e method and especially about t h e filtering p r o c e d u r e . I am a l s o v e r y g r a t e f u l f o r all comments and r e m a r k s of my colleagues from IIASA, and in particu- l a r , I wish t o acknowledge t h e contribution of Maximilian Posch who pro- vided me with t h e three-dimensional g r a p h i c s used in t h i s p a p e r . 1 would a l s o like t o e x p r e s s my g r a t i t u d e t o Dr. J e r z y Pruchnicki from t h e Institute f o r Meteorology and Water Management in W a r s a w who supported p a r t of t h e work done in W a r s a w . Finally I would like t o thank Dr. Leen Hordijk who a g r e e d to review t h e e n t i r e manuscript of t h i s p a p e r .

ABSTRACT

A numerical Positive Definite Pseudo-Spectral (PDPS) method for t h e solution of t h e advection equation is presented. The method consists of t w o p a r t s . F o r e a c h time s t e p f i r s t a solution using a pseudospectral method i s computed. Then t h e solution i s c o r r e c t e d by a filtering p r o c e d u r e which eliminates negative values. The numerical test with t h e r o t a t i o n a l velocity field a n d d i f f e r e n t initial conditions shows t h a t t h e p r e s e n t method h a s t h e a c c u r a c y of t h e pseudospectral o n e without producing negative values. An additional advantage of t h e PDPS method i s t h e elimination of spurious a r t i f i c i a l shortwaves typical for t h e pseudospectral solution.

-

^vii

-

TABLE OF CONTENTS

1. INTRODUCTION 2. NUMERICAL METHOD

2 . 1 Problem Formulation 2.2 Pseudospectral Solution 2 . 3 Filtering Procedure

2.3.1 Method

2.3.2 A One-Dimensional Example 2.4 Stability and Convergence

3 . ADVECTIVE TEST 3 . 1 Basic Equation

3.2 Cone Shape Initial Condition

3 . 3 Rectangular Block Initial Condition 3 . 4 Smooth Initial Condition

3 . 5 Comparison of Different Initial Conditions 4 . CONCLUSIONS

REFERENCES

AN EFFICIENT POSlTIYE DEFINITE METHOD FOR THE NUMERICAL SOLUTION OF THE ADVECTION EQUATION

Jerzy Bartnicki

1. INTRODUCTION

The p a r t i a l differential advection-diffusion equation i s most frequently used f o r t h e mathematical description of t h e long r a n g e t r a n s p o r t of a i r pollutants. This equation i s a l s o a basic one f o r t h e atmospheric p a r t of t h e IIASA RAINS (Regional Acidification INformation and Simulation) model described by Alcamo et. al. (1985) and Hordijk (1985). The atmospheric module of RAINS consists of t h e source-receptor matrices computed by MSC-W (Meteorological Synthesizing Centre-West) in Oslo, using t h e Long Range T r a n s p o r t (LRT) model developed by Eliassen and Saltbones (1983).

This LRT t r a j e c t o r y model i s used by t h e Co-operative Programme f o r Moni- toring and Evaluation of t h e Long Range Transmission of Air Pollutants in Europe (EMEP) f o r r o u t i n e calculations. In o r d e r t o use t h i s o r any o t h e r LRT model within RAINS, i t i s important t o evaluate t h e uncertainty and

credibility of t h e r e s u l t s . Among d i f f e r e n t t y p e s of s o u r c e s of uncertainty in LRT models t h e error introduced by t h e numerical method used t o solve t h e advection4iffusion equation c a n b e a n important o n e , especially for models with nonlinear chemical reactions.

The main goal of t h i s p a p e r i s to p r e s e n t a numerical method which c a n b e used for t h e solution of t h e advection4iffusion equation without produc- ing negative values. T h e r e f o r e t h e method could b e applied to nonlinear problems as w e l l with high a c c u r a c y typical f o r t h e pseudospectral a p p r o a c h and without losing stability (which o c c u r s when negative values a p p e a r ) . When solving t h e advection4iffusion equation, t h e diffusive p a r t i s r e l a t i v e l y l e s s important t h a n t h e advective p a r t concerning numerical problems. A l s o from t h e physical point of view, in t h e synoptic scale of motions, t h e diffusion term i s small compared to t h e advective one a n d i s even neglected in some models (e.g. in t h e MSC-W model). T h e r e f o r e only a n application of t h e method t o t h e advection equation i s p r e s e n t e d in t h i s p a p e r , however, i t c a n b e used f o r t h e advection4iffusion equation as w e l l .

2. NUMERICAL METHOD

Among many d i f f e r e n t methods used f o r t h e numerical solution of t h e advection equation, t h e s p e c t r a l (Orszag, 1971a) and pseudospectral (Gottlieb and Orszag, 1977) a p p r o a c h are relatively efficient and a c c u r a t e . The a c c u r a c y of t h e s e methods i s b e t t e r compared with finite d i f f e r e n c e methods (Orszag, 1971b), and a l s o to o t h e r methods (Long and P e p p e r , 1981;

Chock, 1985). Another advantage of t h e s p e c t r a l methods i s t h e simple mathematical formulation which makes them convenient f o r p r a c t i c a l appli- cations, especially when using numerical Fast F o u r i e r Transform (FFT) (Coo-

ley and Tukey, 1965). S p e c t r a l and pseudospectral methods have been suc- cessfully applied t o t h e a i r pollution t r a n s p o r t models by Christensen and Prahm (1976) and Wangle et al. (1978). Unfortunately, t h e a c c u r a t e pseu- d o s p e c t r a l and s p e c t r a l methods can produce negative values during t h e numerical solution of t h e advection equation. For many p r a c t i c a l problems, like a i r pollutant t r a n s p o r t involving nonlinear chemistry, t h i s phenomenon makes t h e pseudospectral method unstable. T h e r e are o t h e r methods, like t h e flux-corrected t r a n s p o r t (FCT) method (Boris and Book, 1976; Zalesiak, 1979) and a positive definite algorithm developed by Smolarkiewicz (1984) t h a t c a n b e applied in t h i s case. However, t h e s e methods e i t h e r r e q u i r e a long computational time o r are significantly less a c c u r a t e t h a n t h e pseudos- p e c t r a l solution. This p a p e r p r e s e n t s a combined numerical method: The Positive Definite Pseudo-Spectral (PDPS) method, which eliminates com- pletely negative values on one hand, and i s of t h e same o r d e r of a c c u r a c y as a pseudospectral a p p r o a c h , on t h e o t h e r hand.

2.1. Problem Formulation

The multidimensional advection equation t o b e solved h a s t h e following form:

where c

=

c ( z . t ) is t h e concentration (could b e a r b i t r a r i l y s c a l a r ) , assumed t o b e non-negative.

uj

=

uj ( z , t ) i s t h e j-th velocity component

( 2 , t )

=

( z l,... , z N , t ) a r e t h e s p a c e and t i m e coordinates

The numerical method presented in t h i s p a p e r involves two basic s t e p s at e a c h time s t e p when solving equation ( 1 ) :

( 1 ) The pseudospectral method is applied t o equation (1) at time t a n d a solution which contains a l s o negative values of t h e concentra- tion i s achieved.

( 2 ) The filtering p r o c e d u r e , which removes all negative values of t h e concentration, is used t o g e t t h e solution at time t + A t .

Let c m

=

c ( z .m A t ) b e t h e concentration field with periodic boundary conditions at time m At. W e are looking f o r t h e concentration c m + l

=

c ( z ,( m + l ) A t ) at time ( m + l ) A t in t h e uniform mesh of size M l x M 2 ,

...,

^{X M N} where t h e location of t h e mesh points i s given by:

where

f o r any j

=

1,2

,...,

N.

The pseudospectral method c a n b e r e p r e s e n t e d by a n o p e r a t o r

PI

which, applied t o t h e d i s c r e t e concentration field c m at time m At, p r o d u c e s t h e concentration c-+' at time ( m + l ) A t :

c- + I

=

p l ( p ) (4)

The concentration r ? + l c a n still include negative values. The filtering p r o c e d u r e c a n b e r e p r e s e n t e d by t h e o p e r a t o r

f

which transforms r? t o c m containing non-negative values only:

Thus, t h e positively defined pseudospectral method c a n b e defined as:

In principle t h e o p e r a t o r c a n r e p r e s e n t a l s o o t h e r methods, not only t h e pseudospectral method. However, because of i t s simplicity and accu- r a c y , t h e pseudospectral approximation i s a r a t h e r efficient one f o r t h e numerical solution of equation (1).

2.2. Pseudospectral Solution

The pseudospectral a p p r o a c h developed by Gazdag (1973) h a s been chosen as t h e o p e r a t o r P.

-

The principle of Gazdag's method i s to approxi- m a t e t h e time derivatives by a t r u n c a t e d Taylor s e r i e s , and t h e n r e p l a c e t h e time derivatives by t h e s p a c e derivative terms, which are computed using t h e s p e c t r a l method. Mathematically t h e method c a n b e described as follows. Assuming t h a t w e know t h e concentration c m at time m A t , t h e con- c e n t r a t i o n c m

+'

^at t h e n e x t time s t e p (m +1)At c a n b e approximated by t h e t r u n c a t e d Taylor s e r i e s

Following Gazdag (1973), t h e time derivatives of c c a n b e e x p r e s s e d in terms of t h e s p a c e derivatives of c and uj by making use of equation (1):

The s u p e r s c r i p t m h a s been omitted in t h e above equations f o r con- venience. Equations (8-10) show how to compute any o r d e r time d e r i v a t i v e of c from t h e lower o r d e r time d e r i v a t i v e s of u j a n d c . The f i r s t o r d e r time d e r i v a t i v e of c c a n b e computed d i r e c t l y from t h e basic advection equation.

I t remains only t o compute s p a c e derivatives of c which i s done with t h e s p e c t r a l method. Denoting t h e set of all g r i d points (Equations 2-3) by R , t h e finite F o u r i e r t r a n s f o r m C of c c a n b e written as

where i

=

cland k i s t h e wave v e c t o r k

=

( k l

,...,

k j

,...,

k N )

whose components assume i n t e g e r values within t h e limits

From C ( k ,t ) t h e p a r t i a l d e r i v a t i v e s of c ( z ,t ) c a n b e computed as

The numerical computation of t h e s p a c e d e r i v a t i v e s d e s c r i b e d by Equa- tions (11-14) c a n b e c a r r i e d o u t sufficiently f a s t by t h e u s e of t h e numerical Fast F o u r i e r Transform

(FFT,

Cooley and Tukey, 1965). According t o Gazdag (1973) i t gives v e r y a c c u r a t e r e s u l t s and t h e r e f o r e h e called i t Accurate S p a c e Derivative (ASD) method.

2.3. Filtering Procedure

The pseudospectral method d e s c r i b e d in t h e previous p a r a g r a p h pro- vide t h e concentration in t h e g r i d system at time (m +l)At , assuming t h a t t h e concentration at time m At i s known (also t h e velocity and i t s time derivatives). Unfortunately, t h e new concentration field may contain nega- tive values. The p r e s e n c e of negative concentrations i s a common phenomenon f o r d i f f e r e n t numerical methods used f o r t h e solution of t h e advection equation. According to Adam (1985), t h i s i s mainly due t o t h e wrong numerical propagation s p e e d of t h e s h o r t e s t waves in t h e spectrum.

He suggests, t h a t t h e situation c a n b e improved by applying digital f i l t e r s . However, t h e most common l i n e a r f i l t e r s d o not completely remove t h e nega- t i v e values. The main f e a t u r e s of a p e r f e c t filtering p r o c e d u r e are: (1) To remove negative values. (2) To c o n s e r v e t o t a l m a s s . (3) To p r e s e r v e t h e s h a p e of t h e function. (4) To p r e s e r v e t h e maxima. (5) To b e f r e e of shortwave noise. Unfortunately, none of t h e existing numerical f i l t e r s s a t i s f y all above requirements.

2.3.1. Method

The multidimensional nonlinear filtering p r o c e d u r e developed in t h i s p a p e r fulfills at least some of t h e conditions mentioned above. I t completely removes negative values and c o n s e r v e s t h e total mass with an a c c u r a c y of 0.001%. Filtered maxima and t h e s h a p e of t h e function a r e relatively close t o t h e original ones. The p r o c e d u r e c a n b e explained as follows. Let c, b e t h e concentration in t h e j - t h point of t h e one-dimensional g r i d system con- sisting of N points ( j

=

1 , .

..

,N). If all c, values are non-negative t h e f i l t e r does not change them. Let u s assume now t h a t t h e concentration field h a s

N 1

positive values (cj >O),

N 2

z e r o values (c,

=

0) and

N 3

negative values (cj

<

0). Obviously

A s assumed under Equation (1) we have:

Ml > M 3

where

is t h e "positive" m a s s and

i s t h e "negative" m a s s . With t h e above assumptions t h e filtering p r o c e d u r e i s defined by t h e following algorithm:

1. Compute t h e negative mass

M 3

and check if i t i s g r e a t e r than zero.

If not, stop.

2. Compute t h e number of positive concentrations

N1.

3. Check t h e sign of t h e concentration c j f o r j

=

1 ,

..., ^N

(a) If c j

>

^0, s u b t r a c t t h e negative mass divided by t h e number of positive concentrations: c j := c j

- - M3

N

¹

(b) If i t is z e r o , d o nothing.

(c) If i t i s negative, s e t i t t o zero: c j := 0.

2.3.2.

A

One-Dimensional Example

The f i l t e r i n g p r o c e d u r e lined o u t in t h e p r e v i o u s p a r a g r a p h i s illus- t r a t e d b y a simple one-dimensional example with a g r i d system consisting of 11 points. The initial d i s t r i b u t i o n shown in Figure l a i s t y p i c a l f o r t h e i n t e r m e d i a t e solution of t h e advection equation with "delta" function (con- c e n t r a t i o n s at a l l points e x c e p t o n e are e q u a l to z e r o ) as initial condition.

Two negative values of t h e c o n c e n t r a t i o n are p r e s e n t in t h e distribution: -4 at point number 4 a n d -5 at point number 8. After t h e f i r s t i t e r a t i o n (Figure l b ) only o n e n e g a t i v e value remains: -0.8 at point number 11. The s e c o n d a n d final i t e r a t i o n (Figure I c ) gives a d i s t r i b u t i o n without negative values.

The maximum i s s l i g h t l y lower: 13 instead of 15 b u t t h e s h a p e of t h e final d i s t r i b u t i o n i s q u i t e c l o s e to t h e initial o n e (Figure I d ) . From F i g u r e I d i t c a n b e also s e e n t h a t t h e s h o r t waves p r e s e n t in t h e initial d i s t r i b u t i o n h a v e b e e n removed from t h e final one.

The b a s i c f e a t u r e of t h e algorithm p r e s e n t e d a b o v e i s t h e c o n s e r v a t i o n of mass, which c a n b e e x p r e s s e d as

M I

-

^{M 3}

⁼

^const. (18)

The algorithm i s c o n v e r g e n t a n d s t a b l e (this will b e p r o v e d in t h e n e x t sec- tion), a n d a l s o simple in i t s numerical realization. Numerical e x p e r i m e n t s with d i f f e r e n t initial d i s t r i b u t i o n s i n d i c a t e t h a t t h e typical number of i t e r a - tions n e c e s s a r y to a c h i e v e a non-negative d i s t r i b u t i o n i s not g r e a t e r t h a n t w o . Also t h e additional computer-time s p e n t f o r f i l t e r i n g i s small ( 5 109.) compared to t h e computer-time r e q u i r e d b y t h e p s e u d o s p e c t r a l method.

INITIAL CONCENTRATION

-6

1 2 3 4 5 6 7

GRID POINT NUMBER

Figure l a . One dimensional t e s t f o r t h e f i l t e r i n g p r o c e d u r e : Initial dis- t r i b u t i o n of t h e c o n c e n t r a t i o n .

CONCENTRAT ION AFTER FIRST ITERATION

CONCENTRATION 6 4

F i g u r e l b .

5 6 7

GRID POINT NUMBER

One dimensional test f o r the f i l t e r i n g p r o c e d u r e : Distribu- t i o n a f t e r f i r s t i t e r a t i o n .

CONCENTRAT ION AFTER SECOND ITERATION

CONCENTRATION 6

4

-6

1 2 3 4 5 6 7 8 9 10 1 1

GRID POINT NUMBER

Figure l c . One dimensional test f o r t h e f i l t e r i n g p r o c e d u r e : Distribu- tion a f t e r second and final i t e r a t i o n .

CONCENTRATION AFTER SECOND ITERATION

CONCENTRATION 6 4

Figure l c .

1 2 3 4 5 6 7 B 9 10 1 1

GRID POINT NUMBER

One dimensional t e s t f o r t h e filtering procedure: Distribu- tion a f t e r s e c o n d and final iteration.

I NI TI AL AND FINAL CONCENTRAT ION

CONCENTRATION 6 4

-6

1 2 3 4 5 6 7 8 9 10 1 1

GRID POINT NUMBER

4- CONCENTRATION AFTER SECOND ITERATION

I

I

.!

...

,)

+--$- I

Figure I d . One dimensional t e s t f o r t h e filtering p r o c e d u r e : Comparis- o n of initial and final distribution of t h e c o n c e n t r a t i o n .

2.4. Stability and Convergence

The PDPS method d e s c r i b e d by equation (6) i s a s u p e r p o s i t i o n of t h e

A

o p e r a t o r P

--

t h e p s e u d o s p e c t r a l method

-

a n d

-

t h e f i l t e r i n g pro- c e d u r e . The s t a b i l i t y of t h e p s e u d o s p e c t r a l method i s discussed in d e t a i l by Gazdag (1973). He p r o v e d t h a t t h e s t a b i l i t y condition i s s a t i s f i e d f o r t r u n - c a t e d T a y l o r s e r i e s of o r d e r 3,4,7 a n d 8.

F o r t h e f i l t e r i n g p r o c e d u r e t h e r e are t h r e e possibilities: ( 1 ) The ini- t i a l value of t h e c o n c e n t r a t i o n i s n e g a t i v e (Zj

<

^{0 )} a n d becomes z e r o a f t e r f i l t e r i n g (cj

=

^{0 ) .} ( 2 ) The initial c o n c e n t r a t i o n i s e q u a l to z e r o ( z j

=

^{0 )} a n d r e m a i n s z e r o (cj

=

0 ) . ( 3 ) The initial c o n c e n t r a t i o n i s positive (Ej

>

^{0 ) a n d}

finally r e m a i n s non-negative, b e c a u s e t h e p a r t of t h e n e g a t i v e mass sub- t r a c t e d f r o m i t c a n n o t b e l a r g e r t h a n t h e o r i g i n a l v a l u e (0 S cj

<

^{Z j ) .}

T h e r e f o r e , t h e following condition i s fulfilled by t h e f i l t e r i n g p r o c e d u r e : 0 5 cj 5

I

Ej

1

^j

⁼

^1,

^...,

^N

where

cj

- -

c o n c e n t r a t i o n at point j b e f o r e f i l t e r i n g

cj

-

c o n c e n t r a t i o n at point j a f t e r f i l t e r i n g

A simple implication of r e l a t i o n ( 1 9 ) i s t h a t t h e amplification f a c t o r s are less t h a n one. This means t h a t t h e f i l t e r i n g p r o c e d u r e i s s t a b l e a n d also t h a t t h e PDPS method, as a s u p e r p o s i t i o n of two s t a b l e o p e r a t o r s , s a t i s f i e s t h e s t a b i l i t y condition.

The f i l t e r i n g p r o c e d u r e i s also c o n v e r g e n t . This i s obvious when only non-negative v a l u e s are p r e s e n t in t h e initial d i s t r i b u t i o n . L e t u s assume now a n initial d i s t r i b u t i o n with N 1 positive values, N 2 z e r o v a l u e s a n d N 3 n e g a t i v e v a l u e s of t h e c o n c e n t r a t i o n in t h e initial d i s t r i b u t i o n . Due to

Equation (16) t h e negative m a s s M 3 is smaller than positive m a s s MI.

According t o t h e filtering p r o c e d u r e all negative values become equal t o z e r o and z e r o values are not changed. From e a c h positive value of t h e con- c e n t r a t i o n t h e negative m a s s a v e r a g e d o v e r t h e number of positive values i s s u b t r a c t e d . If e a c h positive value i s g r e a t e r t h a n t h e a v e r a g e negative mass, t h e filtering p r o c e d u r e i s completed a f t e r t h e f i r s t i t e r a t i o n . If not, t h e r e are positive values lower than t h e a v e r a g e negative mass, and they become equal t o z e r o during t h e second iteration. I t means t h a t a f t e r e a c h i t e r a t i o n t h e number of z e r o s i n c r e a s e s at l e a s t by t h e number of negative values. Assuming t h a t t h e filtering p r o c e d u r e i s not convergent, t h e r e will b e z e r o s only a f t e r l e s s t h a n N -N3 i t e r a t i o n s , which i s impossible because of t h e conservation of m a s s (Eq. 18). Thus t h e filtering p r o c e d u r e is con- vergent.

3. ADVECTIVE TEST

In o r d e r t o c h e c k t h e a c c u r a c y of t h e method d e s c r i b e d in t h e previ- ous p a r a g r a p h , a numerical advective test h a s been performed. A s t a n d a r d a r t i f i c i a l velocity field h a s been used with t h e "frozen" initial s h a p e moving around t h e a x i s of rotation. T h r e e d i f f e r e n t initial conditions have been chosen: cone, r e c t a n g u l a r block and smooth s h a p e . The test w a s performed both f o r t h e Positive Definite Pseudo-Spectral (PDPS) method, and Pseudo- S p e c t r a l (PS) a p p r o a c h (Gazdag, 1973).

3.1. Basic Equation

The equation describing t h e r o t a t i o n of t h e "frozen" initial condition h a s been frequently used f o r testing numerical methods (Orszag, 1971a; Gaz- dag, 1973; Long and P e p p e r , 1981; Christensen and Prahm, 1976). I t h a s t h e

following form:

where o is angular velocity

and T i s t h e period of rotation. Equation (1) w a s solved numerically on a g r i d consisting of 32 X 32 points. The time s t e p w a s equal

-

which means

400

t h a t one full revolution r e q u i r e d 400 time s t e p s . The analytical and numeri- c a l solutions were compared a f t e r 1 0 rotations. In addition, s e v e r a l param- eters were computed during e a c h r u n . Namely:

(1) Mass conservation (in X )

- ^M

where c, ( i

,

^j ) is t h e initial concentration

(2) Conservation of t h e s q u a r e of t h e m a s s (in I.)

-

^SU

(3) Minimum of c ( i , j )

-

/UI?Y

(4) Maximumof c ( i , j ) - M A X (5) Maximum a b s o l u t e e r r o r

-

^MER

MER

=

m a x ( I c ( i , j )

-

c , ( i , j ) l ) i lj

(6) Average a b s o l u t e error

-

AER 1 32 32

AER

=- C C

I c ( i l j >

-

c , ( i , j ) l 3 2 x 3 2

,=,

^{j = ,}

All a b o v e p a r a m e t e r s a r e functions of time and a r e d i f f e r e n t f o r e a c h initial condition. The maximum of e a c h t e s t e d initial condition w a s k e p t c o n s t a n t a n d e q u a l 100.

3.2. Cone Shape Initial Condition

The "cone" s h a p e initial condition (Figure 2 a ) i s a s t a n d a r d o n e a n d was applied as a t e s t c a s e t o almost a l l numerical methods used f o r solving a n advection equation. In t h e g r i d system t h e "cone" s h a p e i s defined as:

In F i g u r e s (2b) a n d (2c) t h e numerical solutions a f t e r t e n r o t a t i o n s a r e shown f o r t h e PDPS a n d P S methods, r e s p e c t i v e l y . The d i f f e r e n c e in s h a p e s i s small a n d both numerical solutions a r e q u i t e c l o s e t o t h e a n a l y t i c a l one.

However, negative values a p p e a r in t h e P S solution.

The m a s s c o n s e r v a t i o n M , defined by Equation (22), i s equal t o 100%

d u r i n g t h e e n t i r e r u n f o r both PDPS a n d P S with a c c u r a c y b e t t e r t h a n 0.001%. The s q u a r e mass c o n s e r v a t i o n S U , defined by Equation (23), i s

shown in F i g u r e 3. The s q u a r e mass i s well c o n s e r v e d b y t h e P S method (99.7% a f t e r 1 0 r o t a t i o n s ) and slightly worse b y t h e PDPS method (92.6%

a f t e r 1 0 r o t a t i o n s ) . In t h e l a t t e r case t h e s q u a r e mass d e c r e a s e s r a p i d l y d u r i n g t h e f i r s t r o t a t i o n and t h e n s t a y s almost at t h e same level.

The minimum values MJN are shown in F i g u r e 4 f o r b o t h methods. In case of PDPS n e g a t i v e v a l u e s are n o t c r e a t e d and t h e numerical minimum i s equal to t h e a n a l y t i c a l o n e , which i s z e r o . In case of P S n e g a t i v e values are c r e a t e d , r e a c h i n g -1.81 a f t e r t e n r o t a t i o n s .

The a n a l y t i c a l maximum

MAX

i s e q u a l to 1 0 0 a n d i s s l i g h t l y a b o v e t h e numerical o n e s (Figure 5 ) . A f t e r t e n r o t a t i o n s t h e maximum f o r t h e PDPS method i s e q u a l to 91.75 w h e r e a s i t i s 94.02 f o r t h e P S method. F o r b o t h methods t h e maximum d e c r e a s e s mainly d u r i n g t h e f i r s t r o t a t i o n and t h e n s t a y s at t h e same level.

F o r both PDPS a n d P S t h e maximum a b s o l u t e error MER, defined by Equation (24), o c c u r s at t h e t o p of t h e c o n e . I t i s slightly h i g h e r f o r PDPS t h a n f o r P S (Figure 6 ) , a n d i s l e s s t h a n 1 0 a f t e r t e n r o t a t i o n s .

F o r t h e P S method t h e a v e r a g e a b s o l u t e e r r o r AER, defined by Equa- tion (25), i n c r e a s e s r a p i d l y d u r i n g t h e f i r s t r o t a t i o n a n d t h e n , with some fluctuations, r e m a i n s at t h e same l e v e l of 0.14 (Figure 7 ) . In case of t h e PDPS method, AER i n c r e a s e s slowly, r e a c h i n g 0.172 a f t e r t e n r o t a t i o n s .

9.9. Rectangular Block Initial Condition

The ' R e c t a n g u l a r Block" initial condition i s shown in F i g u r e Ba. I t i s defined on t h e g r i d as:

I

1 0 0 if 5 5 i 511 a n d 135 j 5 1 9

c ( i , j )

=

₀ o t h e r w i s e i , j

=

1,...,32

Figure 2 a . S h a p e of t h e c o n e a f t e r 10 r o t a t i o n s : a n a l y t i c a l solution.

Figure 2b. S h a p e of t h e cone a f t e r 10 rotations: PDPS method.

Figure 2 c . Shape of the cone a f t e r 10 rotations: PS method.

CONE - SQUARE OF MASS CONSERVATION

0

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Figure 3. S q u a r e of mass conservation, with t h e cone s h a p e a s initial condition.

CONE - MI NlMLlM

0

.o

-0.2 -0.4 -0.6 -0.8 M/N

-

1

.o

-1.2 -1.4 -1.6 -1.8 -2

.o

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

F i g u r e 4 . Minimum v a l u e s f o r t h e PDPS a n d PS methods with t h e c o n e s h a p e i n i t i a l condition.

CONE - MAXIMUM

30 20

.

10

.

0

-

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Figure 5 . Maximum v a l u e s f o r t h e PDPS a n d PS methods with t h e c o n e s h a p e initial condition.

CONE - ^MAXIMUM ABSOLUTE ERROR

Figure 6. Maximum absolute e r r o r f o r t h e PDPS a n d P S methods with t h e c o n e s h a p e initial condition.

F i g u r e 7.

CONE - AVERAGE ABSOLUTE ERROR

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

A v e r a g e a b s o l u t e e r r o r f o r t h e P D P S a n d PS m e t h o d s with t h e c o n e s h a p e i n i t i a l condition.

The numerical solution a f t e r t e n r o t a t i o n s is shown in Figure 8 b f o r t h e PDPS method, a n d in F i g u r e 8c f o r t h e P S method. The d i f f e r e n c e between t h e a n a l y t i c a l a n d t h e numerical solutions i s g r e a t e r t h a n f o r t h e c o n e s h a p e initial condition discussed e a r l i e r , b u t t h e initial s h a p e i s k e p t q u i t e well. An i m p o r t a n t a d v a n t a g e of t h e PDPS method i s t h e a b s e n c e of s h o r t w a v e noise, p r e s e n t in t h e solution given by t h e P S method.

The s q u a r e mass i s b e t t e r c o n s e r v e d by t h e P S method (95% a f t e r t e n r o t a t i o n s ) t h a n by t h e PDPS method (69.36% a f t e r t e n r o t a t i o n s ) . Again, l i k e in t h e case of t h e c o n e s h a p e , t h e s q u a r e mass d e c r e a s e s mainly d u r i n g t h e f i r s t r o t a t i o n (71.42%) f o r PDPS, and t h e n remains a t t h e same level 70%

(Figure 9).

The minimum f o r t h e PIIPS method equals t o z e r o d u r i n g t h e e n t i r e r u n . F o r P S i t v a r i e s from -9.28 a f t e r t h e t h i r d r o t a t i o n t o -13.07 a f t e r t e n r o t a - tions (Figure 10). The minimum value p r o d u c e d by t h e PS method is h i g h e r (in a b s o l u t e value) in case of t h e r e c t a n g u l a r block initial condition t h a n in case of t h e c o n e s h a p e initial condition.

The maximum f o r t h e PDPS method (Figure 1 1 ) i n c r e a s e s t o 105.6 a f t e r t h e f i r s t r o t a t i o n a n d t h e n continuously d e c r e a s e s t o 101.0 a f t e r t e n r o t a - tions. The maximum f o r t h e P S method (Figure 1 1 ) i s r e l a t i v e l y high, r e a c h e s 120.16 a f t e r t h e f o u r t h r o t a t i o n a n d remains lower a f t e r w a r d s .

The maximum a b s o l u t e e r r o r i s much h i g h e r compared t o t h e c o n e s h a p e initial condition, both f o r PDPS and P S (Figure 1 2 , c f . Figure 6). In case of PDPS t h e maximum a b s o l u t e error i n c r e a s e s r a p i d l y t o 38.91 a f t e r t h e f i r s t r o t a t i o n a n d t h e n slowly g o e s t o 47.08 a f t e r t e n r o t a t i o n s . The maximum a b s o l u t e e r r o r i n c r e a s e s a l s o in case of t h e P S method (Figure 1 2 )

Figure 8 b . Shape of the rectangular block a f t e r 10 rotations: PDPS method.

Figure 8 c . Shape of t h e r e c t a n g u l a r block a f t e r 10 rotations: P S method.

Figure 9.

RECT. BLOCK - SQUARE OF MASS CONSERVATION

0 1 2 3 4 5 6 7 8 9 10

ROTATION

S q u a r e of mass c o n s e r v a t i o n , with t h e r e c t a n g u l a r block in- i t i a l condition.

t o 38.8 a f t e r s e v e n r o t a t i o n s .

The a v e r a g e a b s o l u t e error (Figure 1 3 ) i s a b o u t 1.5 times s m a l l e r f o r t h e PDPS (0.161-0.181) t h a n f o r t h e PS method (0.232-0.314).

Compared t o t h e c o n e initial condition, t h e r e c t a n g u l a r block initial condition i s a more c r i t i c a l test f o r t h e numerical methods. Differences between numerical a n d a n a l y t i c a l solution are h i g h e r a n d n e g a t i v e numbers are b i g g e r . In t h e P S solution, t h e r e are a l s o s h o r t r a n g e waves p r e s e n t , which did n o t o c c u r b e f o r e . In t h i s case t h e PDPS method p a s s e d t h e test q u i t e well, a n d e s p e c i a l l y , i t p r e s e r v e d t h e numerical maximum c l o s e t o a n a l y t i c a l o n e a n d did n o t p r o d u c e t h e s h o r t w a v e noise, p r e s e n t in PS solu- tion.

3.4. Smooth Initial Condition

C o n t r a r y to t h e two p r e v i o u s cases t h e last numerical test was p e r - formed with t h e following smooth initial condition:

The s h a p e of t h e d i s t r i b u t i o n defined by Equation (28) i s shown in Fig- u r e 1 4 a .

The numerical solution a f t e r t e n r o t a t i o n s i s shown in Figure 1 4 b f o r t h e PDPS method, a n d in F i g u r e 1 4 c f o r t h e PS method. In b o t h cases t h e d i f f e r e n c e s between t h e a n a l y t i c a l a n d t h e numerical solutions are s m a l l . In case of t h e P S method n e g a t i v e values a p p e a r e d a g a i n b u t s h o r t w a v e s c a n n o t b e s e e n o n t h e g r i d . The s h o r t w a v e s are also n o t p r e s e n t in t h e PDPS solution.

RECT. BLOCK - MINIMUM

-

¹⁴

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Figure 10. Minimum values f o r t h e PDPS and P S methods with t h e r e c - t a n g u l a r block initial condition.

RECT. BLOCK - MAXIMUM

112 MAX

108

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

F i g u r e 11. Maximum v a l u e s for t h e PDPS a n d P S methods with t h e r e c - t a n g u l a r block initial condition.

RECT. BLOCK - MAXIMUM ABSOLUTE ERROR

Figure 12.

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Maximum a b s o l u t e e r r o r f o r t h e PDPS and P S methods with t h e r e c t a n g u l a r block initial condition.

RECT. BLOCK - AVERAGE ABSOLUTE ERROR

Figure 13. Average absolute e r r o r f o r t h e

PDPS

and PS methods with t h e r e c t a n g u l a r block initial condition.

Figure 1 4 a . Smooth s h a p e a f t e r t e n r o t a t i o n s : analytical solution.

Figure 14b. Smooth s h a p e a f t e r ten rotations: PDPS method.

Figure 1 4 c . Smooth s h a p e a f t e r t e n rotations: P S method.

The s q u a r e mass, shown in Figure 15, i s well c o n s e r v e d by both t h e PDPS method (98.91% a f t e r t e n r o t a t i o n s ) and t h e PS method (99.76% a f t e r t e n r o t a t i o n s ) . Again, like in both p r e v i o u s c a s e s , t h e s q u a r e mass d e c r e a s e s mainly d u r i n g t h e f i r s t r o t a t i o n s f o r t h e PDPS method.

The minimum i s z e r o f o r t h e PDPS method d u r i n g e n t i r e r u n (Figure 16). F o r P S , i t slowly d e c r e a s e s from -0.19 a f t e r t h e f i r s t r o t a t i o n t o -0.43 a f t e r t e n r o t a t i o n s . However, t h e a b s o l u t e values of t h e minimum are r a t h e r low compared t o t h e c o n e s h a p e a n d r e c t a n g u l a r block initial conditions.

The maximum, shown in Figure 1 6 , i s v e r y c l o s e t o t h e a n a l y t i c a l value 1 0 0 f o r both methods. A f t e r t e n r o t a t i o n s t h e maximum i s e q u a l t o 99.32 f o r t h e PDPS method a n d 99.70 f o r t h e PS method.

The maximum a b s o l u t e e r r o r (Figure 1 7 ) is 0.70 f o r t h e PDPS method a n d 0.69 f o r t h e P S , a f t e r t e n r o t a t i o n s . I t i n c r e a s e s f a s t e r f o r t h e P S method (0.27 a f t e r t h e f i r s t r o t a t i o n ) t h a n f o r PDPS (0.53 a f t e r t h e f i r s t r o t a t o n ) .

The a v e r a g e a b s o l u t e e r r o r is p r a c t i c a l l y t h e same f o r both methods a n d i t slowly i n c r e a s e s from 0.04 a f t e r t h e f i r s t r o t a t i o n t o 0.05 a f t e r t e n r o t a t i o n s .

3.5. C o m p a r i s o n o f D i f f e r e n t Initial C o n d i t i o n s

The r e s u l t s of t h e a d v e c t i v e test depend both on t h e numerical method a p p l i e d t o t h e advection equation a n d on t h e s h a p e of t h e initial condition.

F o r e a c h of t h e t h r e e d i f f e r e n t initial conditions t h e PDPS a n d PS methods c o n s e r v e initial m a s s with a n a c c u r a c y b e t t e r t h a n 0.001%. Also, f o r a l l of them, t h e PDPS method d o e s n o t p r o d u c e n e g a t i v e v a l u e s a n d t h i s i s t h e most important f e a t u r e of t h e method. However, t h e values of t h e o t h e r

SMOOTH SHAPE - SQUARE OF MASS CONSERVATION

ROTATION

Figure 15. S q u a r e of mass c o n s e r v a t i o n : smooth initial condition.

SMOOTH SHAPE - MI NlMUM

-0.45

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Figure 16. Minimum values f o r t h e PDPS a n d P S methods with t h e smooth initial condition.

SMOOTH SHAPE - MAXIMUM

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Figure 17. Maximum v a l u e s f o r t h e PDPS a n d P S methods with t h e smooth initial condition.

m e a s u r e s (SM, MlN f o r t h e P S method, MAX, MER and AER), defined in Sec- tion 3.1, depend o n t h e s h a p e of initial condition.

The s q u a r e mass i s b e t t e r c o n s e r v e d (Figure 15) f o r t h e smooth s h a p e t h a n f o r t h e c o n e s h a p e (Figure 3 ) a n d t h e r e c t a n g u l a r block (Figure 9).

A f t e r t e n r o t a t i o n s , only 0.46% of t h e s q u a r e mass i s l o s t in t h e P S solution a n d 1.09% in t h e PDPS solution. F o r t h e r e c t a n g u l a r block t h e s q u a r e mass d e c r e a s e a f t e r t e n r o t a t i o n s is: 4.28% f o r P S a n d 30.64% f o r PDPS. The c o r r e s p o n d i n g numbers f o r t h e c o n e initial conditions are: 0.93% f o r P S a n d 7.33% f o r PDPS.

The minimum v a l u e s g e n e r a t e d by t h e P S method are smaller (Figure 1 6 ) ,

--

in a b s o l u t e units

-

f o r t h e smooth s h a p e (-0.43 a f t e r s e v e n r o t a t i o n s ) t h a n f o r t h e r e c t a n g u l a r block (-13.07 a f t e r t e n r o t a t i o n s ) , a n d f o r t h e c o n e

(-0.181 a f t e r t e n r o t a t i o n s ) .

Similarly, t h e maximum f o r both methods i s closer to t h e a n a l y t i c a l solution in c a s e of t h e smooth s h a p e (99.32 f o r PDPS, 99.70 f o r PS; F i g u r e 1 7 ) t h a n in case of t h e r e c t a n g u l a r block (101.00 f o r PDPS, 114.04 f o r PS), a n d t h e c o n e s h a p e (91.45 f o r PDPS, 94.02 f o r PS).

The maximum a b s o l u t e error a f t e r t e n r o t a t i o n s (Figure 1 8 ) i s a l s o s m a l l e r f o r t h e smooth s h a p e (0.70 f o r PDPS, 0.69 f o r PS) t h a n f o r t h e rec- t a n g u l a r block (47.08 f o r PDPS, 35.08 f o r P S ) a n d t h e c o n e (8.55 f o r PDPS, 5.98 f o r PS).

Finally, t h e a v e r a g e a b s o l u t e error (Figure 1 9 ) i s slightly h i g h e r f o r t h e smooth s h a p e (0.05 f o r PDPS a n d PS, a f t e r t e n r o t a t i o n s ) t h a n f o r t h e c o n e s h a p e (0.172 f o r PDPS, 0.146 f o r PS) a n d t h e r e c t a n g u l a r block (0.181 f o r PDPS, 0.301 f o r PS).

Figure 18.

SMOOTH SHAPE - MAXIMUM ABSOLUTE ERROR

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

Maximum a b s o l u t e e r r o r Tor t h e

PDPS

a n d PS methods with t h e smooth initial condition.

SMOOTH SHAPE - AVERAGE ABSOLUTE ERROR

0 1 2 3 4 5 6 7 8 9 1 0

ROTATION

F i g u r e 19. A v e r a g e a b s o l u t e error f o r t h e PDPS a n d P S methods with t h e smooth i n i t i a l condition.

Comparing d i f f e r e n t intial conditions, i t seems t h a t t h e r e c t a n g u l a r block s h a p e is t h e most c r i t i c a l test f o r t h e numerical methods. I t i s a l s o confirmed by t h e generation of t h e a r t i f i c i a l shortwave noise, with t h e high amplitude f o r t h e r e c t a n g u l a r block initial condition (Figure 8b-c), smaller f o r t h e cone initial condition (Figure 2b-c) and practically invisible for t h e smooth initial condition (Figure 14b-c).

4. CONCLUSIONS

The PDPS method p r e s e n t e d in t h i s p a p e r is simple and comprehensive both in mathematical formulation and in p r a c t i c a l application. I t does not produce negative values and c o n s e r v e s initial m a s s with 100% a c c u r a c y . The method consists of two basic p a r t s : (1) The pseudospectral solution, and (2) t h e filtering p r o c e d u r e . Compared t o t h e pseudospectral a p p r o a c h t h e additional computer-time f o r t h e PDPS method i s only about 10% higher. The multidimensional filtering p r o c e d u r e is g e n e r a l enough to b e combined with methods o t h e r t h a n PS, and especially with explicit time integration algo- rithms. The PDPS method can a l s o b e applied t o t h e advection-diffusion equation in t h e same way as t o t h e advection equation.

From t h e numerical t e s t s , performed with t h e PDPS and P S methods in t h e rotational velocity field, i t seems t h a t

-

f o r both methods

-

r e s u l t s depend on t h e initial condition. The most commonly used cone s h a p e gives relatively good r e s u l t s concerning a c c u r a c y and shortwave noise. From t h e t h r e e d i f f e r e n t conditions t e s t e d , t h e most c r i t i c a l one i s t h e r e c t a n g u l a r block s h a p e with v e r y s t e e p gradients. In t h i s c a s e , t h e magnitude of t h e negative values g e n e r a t e d by t h e PS method i s l a r g e r (Figure 10) t h a n f o r t h e o t h e r initial conditions (Figures 8c and 14c). Also t h e amplitude of t h e

s h o r t w a v e s o n t h e e n t i r e g r i d system i s l a r g e r (Figure 8 c ) . Both, negative values a n d s h o r t w a v e s noise ( e x c e p t some small d i s t u r b a n c e s c l o s e t o t h e r e c t a n g u l a r block) a r e n o t p r e s e n t in t h e PDPS solution.

Also t h e a v e r a g e a b s o l u t e e r r o r i s smaller f o r t h e PDPS method, e x c e p t in t h e c o n e c a s e , when i t i s slightly h i g h e r compared t o t h e P S solu- tion.

Summarizing, t h e a c c u r a c y of t h e PDPS method i s v e r y c l o s e t o t h e P S one. The a d v a n t a g e of PDPS i s a complete elimination of n e g a t i v e values from t h e solution and t h e r e f o r e i t s possible application t o non-linear p r o b - lems (e.g. chemical r e a c t i o n s d u r i n g t h e t r a n s p o r t ) . An additional advan- t a g e of PDPS i s t h e a b s e n c e of shortwaves, t y p i c a l f o r t h e P S solution, in case of s t e e p g r a d i e n t s in t h e c o n c e n t r a t i o n field.

REFERENCES

Adam, Y. 1985. Nonlinear instability in advection-diffusion numerical models. Appl. Math. Modeling 9 , 434-440.

Alcamo, J., L. Hordijk, J. Kamari, P . Kauppi, M. P o s c h and E. Runca. 1985.

I n t e g r a t e d a n a l y s i s of acidification in E u r o p e . J . Env. Management 21, 47-61.

Boris, J.P. a n d D.L. Book. 1976. F l u x - c o r r e c t e d t r a n s p o r t 111. Minimal error FCT algorithms. J . Comp. Phys. 20, 397-431.

Chock, P.P. 1985. A comparison of numerical methods f o r solving t h e a d v e c - tion equation 11. A t m o s . Environ. 1 9 , 571-586.

C h r i s t e n s e n , 0. a n d L.P. P r a h m . 1976. A p s e u d o s p e c t r a l model f o r d i s p e r - sion of a t m o s p h e r i c pollutants. J. Appl. Meteor. 1 5 , 1284-1294.

Cooley, J.W. a n d J.W. Tukey. 1965. An algorithm f o r t h e machine c a l c u l a t i o n of complex F o u r i e r s e r i e s . Math. Comp. 1 9 , 297-301.

Eliassen, A., a n d J. S a l t b o n e s . 1983. Modeling of long r a n g e t r a n s p o r t of sul- f u r o v e r E u r o p e : a two-year model r u n a n d some model e x p e r i m e n t s . Atmos. Environ. 1 7 , 1457-1473.

Gazdag, J. 1973. Numerical c o n v e c t i v e s c h e m e s b a s e d on a c c u r a t e computa- tion of s p a c e d e r i v a t i v e s . J. Comp. Phys. 1 3 , 100-113.

Gottlieb, D. and S.A. Orszag. 1977. Numerical analysis of s p e c t r a l methods:

Theory and Application Society f o r Industrial and Applied Mathematics.

Philadelphia, Pennsylvania 19103.

Hordijk, L., 1985. A model f o r evaluation of a c i d deposition in Europe. In:

A. Sydow, M. Thoma a n d R. Vichnevetsky (Eds.), Systems Analysis and Simulation 1985, Vol. 11, Akademie-Verlag Berlin, 1985, pp. 30-39.

Long, P.E. and D.W. P e p p e r . 1981. An examination of simple numerical schemes f o r calculating s c a l a r advection. J . Appl. Meteor. 20, 146-156.

Orszag, S.A. 1971a. Numerical simulation of incompressible flows within sim- p l e boundaries. I. Galerkin ( S p e c t r a l ) R e p r e s e n t a t i o n s Study in Appl.

Math. 5 0 , 293-327.

Orszag, S.A. 1971b. Numerical simulation of incompressible flows within sim- ple boundaries: a c c u r a c y . J. Fluid. Mech. 17, 75-112.

Smolarkiewicz, P.K. 1984. A fully multidimensional positive definite advec- tion t r a n s p o r t algorithm with small implicit diffusion. J. Comp. Phys.

5 4 , 325-362.

Wangle, H., B. Van den Bosk and J.H. Seinfeld. 1978. Solution of a t m o s p h e r i c diffusion problems by p s e u d o s p e c t r a l and o r t h o g o n a l collocation methods. Atmos. Environ. 1 2 , 1201-1032.

Zalesiak, S.T. 1979. Fully multidimensional flux c o r r e c t e d t r a n s p o r t algo- r i t h m s f o r fluids. J. Comp. Phys. 31, 335-362.