Two dimensional optimization of smoothing properties of multistage schemes applied to hyperbolic equations

van

KARMAN INSTITUTE

FOR FLUID DYNAMICS

Technical Note 173

TWO DIMENSIONAL OPTIMIZATION

OF SMOOTHING PROPERTIES

OF MULTISTAGE SCHEMES

4 M

I

r

1991

APPLIED TO

HYP~RBOLIC

EQUATIONS

L.A. Catalano

&

H. Deconinck

July 1990

~A~

-~O~

RHODE SAINT GENESE BELGIUM

von Karman Institute for Fluid Dynamics

Chaussée de Waterloo, 72

B

-

1640 Rhode Saint Genèse-Belgium

Technical Note 173

TWO DIMENSIONAL OPTIMIZATION

OF SMOOTHING PROPERTIES

OF MULTISTAGE SCHEMES

APLIED TO HYPERBOLIC EQUATIONS

L.A. Catalano & H. Deconinck

Table of Contents

Aeknowledgements ... 111

Abstract ... lV 1. INTRODUCTION ... 1

2. ONE DIMENSIONAL OPTIMIZATION ... 3

2.1 A general von Neumann analysis ... 3

2.1.1 First order upwind scheme ... 3

2.1.2 Second order upwind scheme :. K-schemes ... 4

2.1.3 Forward Euler time integration scheme ... : ... 4

2.1.4 Two-stage time integration seheme ... 4

2.1.5 Multi-stage time integration seheme ... 5

2.2 Optimization strategy ... , ... : 5

2.3 Properties of multi-stage sehemes ... 6

2.4 Two-stage scheme optimization ... 7

2.5 Multi-stage scheme optimization ... 8

3. ONE DIMENSIONAL RESULTS ... ... 10

3.1 Comparison with the optimization strategy of [3] ... ... 10

4. TWO DIMENSION AL OPTIMIZATION ... 12

4.1 2D analysis as an extension of 1D case ... 12

4.2 Properties of 2D multi-stage schemes ... 15

5. TWO DIMENSION AL RESULTS ... 17

5.1 Optimal parameters for a selection of schemes ... 17

5.2 Application to scalar 2D test eases ... . 19

5.2.1 Shear problem ... 19

5.2.2 Rotational flow problem ... . . . .. 20

6. APPROXIMATE OPTIMIZATION STRATEGIES ... 22

6.1 Partial two dimensional optimization ... 23

6.2 One dimensional optimization along the convection direction ... 23

6.3 One dimensional optimizati(;>n along the 450 direct ion ' ... '" .... 23

6.4 A theoretical comparison ... 24

7. MULTIGRID TECHNIQUE ... 25

7.1 Elliptie case ... 25

7.2 Hyperbolic case ... 25

7.3 FMG and FAS cycle [7,13] ... 26

7.4 Saw-tooth cycle ... 28

7.5 Oblique shock test ease-numerical formulation ... 29

7.5.1 Space discretization ... 29

7.5.2 Definition of eell eonvection speed ... 30

8. MULTIGRID RESULTS ... 32

9. CONCLUSIONS ... 35

References ... 36

Tables ... 37

Acknowledgements

The present research was carried out as part of the BritejEuram Contract no.

AERO-0003-C, issued by the European Community.

A multigrid relaxation code developed in the frame wo~k of this contract at CWI, Amsterdam, was used as starti;ng point for the pre!?e~t multigrid computations. We thank Prof. P.W. Hemker for kindly providing this code.

The first author is deeply indebted to Prof. M. Napolitano, from University of Bari,

Italy, and Prof. H. Deconinck for involving him in this BritejEuram Contract.

Abstract

This report provides a numerical technique for optimizing the smoothing properties of multi-stage explicit time integration sche~es for

à

general spátial discretization of the 2D advection equation.

The optimized scheines have been implemented as smoothers in a multigrid code and their efficiency is compared with the classical Runge-Kutta schemes and, when possible, with the Forward Euler scheme, applied as many times as the number of stages of the

optimized one.

The multigrid code has been tested' for different first and sec'ond order (MUS CL ) upwind finite volume space discretization operators, using a FAS or a Saw-Tooth cycle.

Keywords : Multigrids Advection equation Runge-Kutta scheme Hyperbolic equation Smoothing analysis

1. INTRODUCTION

Since many years, Runge-Kutta schemes have been applied as smoóthing operator in multigrid methods for centered difference discretizations of hyperbolic equations such as the Euler equations governing compressible inviscid fluid flow, see e.g., the work of Jameson [1]. Gustafsson and Lötstedt [2] show that two different mechanisms are involved in these multigrid solvers for hyperbolic equations : elimination of low frequencies is due to an increased speed of these waves as a result of the coarse grids. High frequencies, however, are eliminated by damping, since the inevitable artificial dissipation operators intro duce a parabolic type of behaviour. Therefore, in [2] it is concluded that the propertiesof the multigrid method for the hyperbolic problem are similar to the ones for an elliptic system when considering only the high frequencies.

The same conclusion carries over to upwind discretizations, but here the dissipation is implicitly defined by the properties of the upwind discretization which has been chosen. This reasoning has led van Leer, Tai and Powell [3] to the idea of optimizing high frequency damping properties of multi-stage schemes for a given upwind space operator. This is achieved by minimizing the smoothing factor of the multi-stage scheme, defined as the

maximal time amplification factor over the high frequency Fourier domain

[f,

7r].

The

parameters in this optimization are the multi-stage coefficients and the timestep.

The analysis in [3] is done for a scalar advection equation in one space dimension and extended to a system of hyperbolic equations in one space dimension. This naturally leads to optimallocal characteristic timesteps, which are different for each of the decoupled characteristic equations.

von Neumann analysis is used to perform the present optimization, as well as in [3]. Therefore the influence of non linearity and of boundary conditions on the amplification factor is discarded, but the complexity' of the problem does not aUow the use of more complete theories. Anyway it is weU known that von Neumann analysis pro duces good results in the high frequency range of the error distribution, where the smoothing factor is defined.

The extension to a 2D scalar advection equation requires a 2D von Neumann analysis and a 2D optimization of the multi-stage coefficients and timestep over the high frequency domain. As the strategy used in [3] is based on optimizing the position of some frequencies which are completely annihilated by the multi-stage scheme, it is not longer applicable in the case of two space dimensions.

As in the lD case the multi-stage coefficients anq the timestep have to be optimized. An alternative strategy is proposed in the present report, which does not impose any zeroes on the amplification factor and does extend to 2D. Moreover , since the new strategy is more general it allows, when applied in lD, to recover schemes with equal or lower smoothing factor than found in [3].

When applied to the 2D advection equation, the optimal parameters are function of the ratio of the Courant numbers in x- and y-direction. For equal mesh spaci,ng in x- and y-direction this ratio is equivalent to the angle between.x-axis and advection speed vector. The new 2D optimization has been first tested on scalar advection test cases by means of single grid computations (chapter 5). It is shown that even on a single grid strong improvements in asymptotic convergence rate are obtained due to the improved damping properties for high frequencies.

The implementation of the new smoothers in a multigrid code for scalar advection is described in chapter 7, where two different cycles (FAS and SAW-TOOTH) are investi-gated. The multigrid code was tested for three different finite volume space discretization schemes, namely first order upwind and two schemes of the I\:-class of second order upwind (fully upwind or upwind biased) schemes (I\: = -1 and I\: =

t).

The test case studied was a 2D Burgers equation with an oblique shock not aligned with the mesh.

2. ONE DIMENSIONAL OPTIMIZATION

2.1 A general von Neumann analysis

The one dimensional linear convection equation is considered:

au au

-+a-

=0

at ax

a>O

(1)

2.1.1 First order upwind scheme

The space operator is discretized by uSlng first order accurate finite differences for simplicity :

au aót

ót at

= -

Óx [u(x, t) - u(x - Óx, t)]

(2)

with the Courant-Friedrichs-Lewy number classically defined as

aót

v =

-Óx (3)

Similarly to [3], eq. (2) is reduced to an ordinary differential equation, by inserting harmonie data u( x)

=

uoeiwx, the spatial wave number being defined as (3

=

wÓx .

(4)

All the informations about the space differencing operator are contained in the complex function

z({3)

=

e -if3 - 1 = f({3)

+

ig({3) (5)

with f({3) = cos({3) - 1 and g({3)

= -sin({3)

for the 1st order upwind differencing. Therefore eq. (2) can be written as

du

- = v).u ). E C

dt ' (6)

). = ~t being the Fourier transform of the spatial differencing operator.

In this way von Neumann analysis of eq. (6) can be performed independently of the type of space discretization.

2.1.2. Second order upwind scheme :~-schemes

The class of secorid order upwind schemes known as the ~-schemes is defined by

~ l+k .

.6.x 8x

= -4-Ju (x +

6x,

t) -

u(x,

t)] +

(

1

-"2

k)

[u(x,

t) -

u(x - Óx,

t)]-

- 4 -1-

k

[u(x - .6. x ,

t) -

u(x - 2.6. x ,

t)]

and the Fourier transform is proportional to

(7)

2.1.3 Forward Euler time integration scheme

If the ordinary differential equation (6) is solved by a forward Euler time integration

(8)

the superscripts

e

and R.

+

1 denoting respectively the old and the new time level, the amplification factor G = u:~l is a lst degree polynomial function of z :

(9)

2.1.4 Two-stage time integration scheme

A two-stage or a so called predictor-corrector scheme can be applied to eq. (6) to discretize the time derivative, instead of (8):

(IDa)

(lOb)

a being a free parameter.

The amplification factor then becomes

or [3], using

z

= v )..~t :

(llb)

2.1.5 Multi-stage time integration scheme

The use of predictor values can be extended to a n-stage gcheme defined as :

(12a) k = 1, ... ,n . (12b) (12c)

with C_n = 1, v being defined on the fuU time step ~t.

As the amplification factor

(13)

is a n-th degree polynomial function of

z

= vz, it ean also be written as

m

Pn(z) =

rr

(1

+

Vk Z

+

okviz2), n=2m (14a)

k=l

or

n = 2m

+

1 (14b)

Therefore a n-stage scheme can be considered as a sequence of m predictor-corrector schemes and eventually one forward Euler step (if n=2m+1). Equations (14a) and (14b) are more convenient than (13) for our purpose.

In the following the eoeffieient p, equal to 0 or 1 respective1y for an even or an odd number of stages, will be used.

2.2 Optimization strategy

In eqs. (14a-b) n coefficients, namely.m parameters 0 and (m+p) parameters Vk, have to be fixed. Classical Runge-Kutta n-stage schemes are n-th order accurate in time

and the coefficients Ck in eq. (12b) are chosen to satisfy this property. As a consequence

classical Runge-Kutta schemes do not control the magnitude of the amplification factor

(} =

IIPn(z)1I

in the high frequency range

[f,7r].

As the damping properties of the scheme depend on the smoothing factor u, defined

as the maximal value of (} in the high frequency range, an optimization technique is here developed to calculate the coefficients which minimize u, thereby giving up the n-th order time accuracy of the original Runge-Kutta scheme.

The mathematical formulation of the minmax problem is:

Uopt

=

~i!1

[

m~x

IIP

n

[z(!3),ä,VJII]

a,1I PEr 2 ,11"]

(15)

In [3]

Pn(z)

is forced to have m+p zeroes

!3o,k

in the high frequency range and the problem is reduced to

Uopt = min [ max

IIP

n

[z(!3),l1o]lI]

Po

P E[t,1I"] _{, .} (16)

being ä and

ïi

expressed as function of the vector

For the case of the 1st order upwind discretization, the simplified problem (16) is analyt-icaIly shown to be equivalent to (15) and therefore to lead to the optimal solution of the fuIl minmax problem (15). The same cannot be said for a higher order K-scheme, as it

will be shown in the result section. Moreover the teçhnique proposed in, [3] is no longer applicable to a 2D analysis, where a zero of (} is a two-component vector

(!30,x,/30,y)

'

For these reasons a different and more general approach, using the original parameters

a l , ... , am , V}, ... , Vm+p 1 is proposed in the following.

2.3 Properties of multi-stage schemes

Analytical ca.lculations performed on the square of the magnitude of the amplification

local extrema in [0,7r]

*

Two of them are located respectively in (3

=

°

and (3 = 7r for symmetry.

Since {l2 is continuous as wen as its derivatives, a minimum (a maximum) is located between two maxima (two minima). Therefore the greatest possible number of maxima in [0,7r] is

N

=

sup int(

nt

1) = m

+

1.

If only the high frequency range [~, 7r] is considered,

N

is still m+1, because the boundary value {l ((3 = ~) has to be included, while the value (l((3 = 0) has to be excluded. For an odd number of stages (n

=

2m

+

1),

N

=

sup int(nt1)

=

m

+

1. Therefore,

for a n-stage scheme (n

=

2m or n

=

2m + 1) the number of maxima N max is less than

or equal to

N

= m

+

1, depending on the choice of

(a,

iJ).

Acondition for the smoothing factor to be minimized is N max

=

N.

2.4 Two-stage seheme optimization

For the case of a 2-stage scheme

N

=

m + 1

=

2. In [3] the minmax problem (16) is solved instead of (15) : the equal excursion principle dictates that the two maxima must be equal. This allows to write a non linear equation in (30 that can be solved by a simple Newton-Raphson procedure:

(17) where Mb k = 0,1, are the

Fr

maxima arid are only function of (30'

Equation (17) can also be imposed in the solution of the fun minmax problem (15) and represents acondition for expressing a as a function of IJ.

Of course this is not sufficient to minimize the function

0"(

a, IJ). Another relation, required to complete the set of two equations in the two unknowns a and IJ, is used to evaluate the value of IJ for which

0"

is minimum:

dO"[a(IJ), IJ] =

°

dIJ . (18)

*

(~

d(3 = sm

. (

. (3) P2m-l(Cos(3

))

Having imposed eq. (17), the presence of discontinuities in the total derivative du[ot),v)]

is avoided and eq. (18) can be written .

. Therefore the following non linear system has to be solved:

MI(a, v)

=

Mo(a, v) da[a(v), v] -~~---=.

=

0 . dv (19) (18)

Equations (19) and (18) can be analytically solved for the case of a lst order upwind differencing operator. For the more general case of a IC-scheme a numeri cal approach is required:

1) a first solution

(aO, va)

with Nmax =

R

= 2 is fucedj

2) a is evaluated from eq. (19) by a scalar Newton-Raphson technique, keeping v con-stantj a new solution

(al,

vI = vO),satisfying MI = Mo = a, is obtainedj

3) the derivative' ~~ (with a = aCv) calcul~ted from (19) for each value of v) is evaluated numericallyj the value of v ean therefore be inereased or decreased (with (19) always

iUlposed) depending on the sign of ~~, until the eondition (18) is satisfied, the function a being minimjzed.

Intermediate values (ai, vi) are calculated, with a( ai, vi) decreasing during the itera-tive procedure.

More details are given in the next section, for the more general case of n stages.

2.5 Multi-stage scheme optimization

Since

R

= m + 1 both for n=2m and n=2m+1, m eonditions

k

=

1, ... ,m (20)

ean be imposed for ca.lculating the m components of

a

as a function of iJ.

The strategy of solving eqs. (20) in

a

instead of iJ is due to the fact that the number of equations (20) is

R,

just like the number of unknowns ak, both for an even and an odd nu.mber of stages.

In this way Cf = Mo = MI = ... = Mrn is continuous with its derivatives and the

solution is imposed to be mini mal by writing m+p (p=O or 1) conditions 'on the m+p

components of the vector

v:

8Cf[ä(z1),

i1]

--=--...:...,.;:..:..-...:.

=

0,

8

V k k = 1, ... ,m+p (21)

The numeri cal procedure described in the previous section can be extended to the general

case of n stages, ä and

v

having to be considered instead of a and v.

Of course step 2), corresponding to the solution of eq. (20) for ä, requires a vectorial

Newton-Raphson method : the differences

k = 1, ... ,m (22)

are defined as the components of the residual vector

R.

The m eqs. (20) are satisfied wh~n

R=Ö.

Let us suppose that i steps of the iterative procedure have been performed and that

the actual solution is (äi_{, Vi); when eqs. (21) have to be solved, a new value}

_ä

i₊1 _{has to be}

calculated. Starting from the solution äi,o =

ä

i

, the classical vectorial Newton-Raphson

sol ver

allows to linearize eqs. (22). The maxima Mk , k

a~jä-I are evaluated numerically.

(23)

0,1, ... , mand the Jacobian matrix

The linear system (23) is solved in

/).ä

j _{by a simple Gauss-Seidel method and the}

value äi,j is updated:

--i)· --i )·-1

+

A --).

a' = a ' ua (24)

This procedure is repeated until

R

= O.

Step 3) consists in a classical gradient method where the derivatives

%::

are continuous

,

and have to be evaluated numerically too. .

3. ONE DIMENSIONAL RESULTS

3.1 Comparison with the optimization strategy of [3]

A strategy for solving the simplified minmax problem (16) is proposed in [3] and applied to the first order upwind scheme and to the higher order K-schemes, for a number of stages going from 1 to 6.

Here a different and more gener al method has been described, which solves directly the full 1D minmax problem (15), without any simplifications. As it will be shown in the next section, this technique can be extended to the optimization of a multi-stage finite difference scheme applied to the 2D advection equation, which is the main goal of this study.

Some results have been obtained for the 1D optimization, as a check, before extending

the method to the, 2D case.

Analytical calculations show that the full minmax problem (15) is equivalent to the simplified eq. (16) for the case of first order upwind differencing. The same results are in fact obtained by using the two different strategies or an analytical approach, as the one described in [3] : the magnitude of the amplification factor is plotted versus the frequency

f3 in figs. (1-5) for a number of stages going from 2 to 6.

The plots show the, properties characterizing the optimal solution the number of

maxima is N max =

N

= m + 1 and the

N

maxima are equal, in all the cases; the

smoothing 'factor is minimized when m+p frequencies f3o,k are completely annihilated by

the time integration scheme, which is the basis of the method used in [3].

The present strategy has then been applied to the optimization of the 3rd order upwind

biased scheme (K =

t) :

the results, compared again with

[3],

are plotted in figs. (6-10)

for n=2, ... ,6. In all the three cases the smoothing factor obtained by using the present

optimization technique (curve 1) is lower than the one obtained by imposing some zeroes

of {! in the high frequency range (curve 2). The ratio between the two values decreases

new optimization is more important when the number of stages is high~ The new optimal coefficients are proposed in Table 1.

Also the second order upwind scheme

(K.

= -1) has been tested by using boththe two strategies and also in this case the optimal solution is not obtained when m+p roots of

the magnitude of the amplification factor are forced in the high frequency range. Anyway the smoothing factor is only slightly reduced by applying the fuH strategy in this case: the

two curves almost coincide, as shown in fig. (11) for n=4.

It is therefore confirmed that the actual optimization technique does lead to equal

or better results than the strategy proposed in [3], when applied to the 1D case. The

extension to the 2D analysis is proposed in the next section.

Figure (6), as weH as figs. (7-10), are also useful tohave a better understanding of the proposed methodology : the values of a(/3o) and v(/3o) calculated by solving the simplified minmax problem (16), as in [3], ean be used as a first approximation of the fuH minmax problem (15). The smoothing factor is then reduced by the gradient step, always keeping

the two maxima constant: in this way the pattern of the curve is modified from 2 to 1 : this means that the solution calculated in [3] is one of the combinations of (a, v) investigated by the present strategy, which is therefore a generalization of [3]; this explains why equal

or better results than [3] are obtained.

4. TWO DIMENSIONAL OPTIMIZATION 4.1 2D analysis as an extension of lD case

The 2D linear advection equation is written in vectorial form as

au

....

-

+ä·

Vu = 0

at

"

(25)

or, in scalar components:'

au

au

,

au

-+a-+b- =0,

Ot

ax

ay

a>

0, b> 0 (26)

where a and b are the Cartesian components of the convection speed vector

ä.

A

generalized convection angle is defined as

()

=

'

atan

(~~x)

=

atan

(V

y

)

'

.

~y a V_x

where V_x and _{v y} are the CFL numbers in x- and y-direction respectively :

a~t b~t

Vx

= --,

~x v y = - - . ~y

Equal mesh spacing, both in x- and y-direction, will be supposed in all the following

analysis and the definition of convection angle ~hen simplifies to

Of course, this hypothesis is not restrictive since all the following considerations are always valid if the advection speed vector components a and bare scaled respectively by the x-and

y-mesh spacings. Equation (26) is again discretized by using lst order accurate upwind

differencing both in x- and y-direction :

au

~t

at

₌

-vx[u(x,y,t) - u(x -

~x,y,t)]

- vy[u(x,y,t) - u(x,y -

~y,t)]. (27)

Since there is only one free parameter, namely ~t, a, b, ~x and ~y being fixed, a unique CFL number is considered

v:;::: v y , if a ~ b

(28b)

together with the CFL ratio

v y b R

=

~

= - =

tan(8), . ' 'V_x ·a .. if a ~ b

(29a)

V_x a R

= -

= - ::;:

cotan( 8), v y b if a ~ b,

(29b)

being ~x

=

~y.

Again without losing generality and completeness the case a ~ b or 0° ~ 8 ~ 45° will be treated in the following. Therefore eq. (27) is rewritten as

&

.

J

ót iJt = -v{[u(x, y, t) - u(x - ~x, y, t)]

+

R[u(x, y, t) - u(x, y - ó.y,

t)]}

.

(30)

Analogously to lD case, harmonie data u(x, y) ::::; uoei(w"x+wI/Y) are inserted in eq.(30),

the x- and y-direction spatial wave numbers being defined as f3x

=

wx~x and f3y

=

Wy~y: (31)

Equation (31) is equivalent to the ordinary differential equation

du

, dt = v À u " À E C, (6) with all the informations about the space discretization contained in the complex function

F(f3x, {3y)

+

i G(f3x, {3y)

=

[f(f3x)

+

R f(f3y)]

+

i [g(f3x)

+

R g(f3y)] =

[f(f3x)

+

i g(f3x)]

+

R [f(f3y)

+

i g(f3y)] , (32)

where f(f3) = cos(f3) - 1 and g({3)

=

-sin(f3) for the lst order upwind differencing, like in the lD case.

Therefore the amplification factor of a multiTstage scheme applied to the 2D convection

equation is again a polynomial of z, where z is function of the two spatial wave numbers

f3x and f3y and of the CFL ratio R.

The domain of definition of

Pn(z)

is the square

Since the magnitude of amplification factor is point symmetric with respect to the origin of the f3x- and f3y-axis, only half of this domain is considered in the following :

Equations (14a-b) are still valid and can be used for ~inimizing the smoothing factor

u, which is again the maximal value of the magnitude of the amplification factor in the

high frequency Fourier domain (see 'fig. 12a) :

The same lD coefficients ak, k::;:: 1, ... , mand Vk, k = 1, ... , m+p have to be optimized; therefore the lD optimization strategy can be extended to the 2D case if the number of possible maxima is

N

= m

+

1. '

The mathematical formulation of the minmax problem is, similarly to chapter 2:

4.2 Properties of 2D multi-stage schemes

It is very important to underline (see eq. (32» that the complex function

zU3:,;, f3y,JJ.)

is a linear combination of the 1D functions

z(f3x)

and

z(f3y),

R being the coefficient of such

comhination. Since the quadratic funçtion

I • . '

.

( • l

is semi-positive definite, a local minimum will be placed where

{?

= 0, namely where

.

"

Let us look at the other local extrema (F:fO or G:fO): the necessary conditions for

H(f3x,f3y)

to be alocal extremum are

aH

....--- =

0

af3x

which cah be written, according to eq. (32), as

aH

- = 0

af3y

aH

=

aH aF

'

af(f3x)

+

aH aG ag(f3x)

,

0

af3x

aF af(f3x) af3x

8C ag(f3x) af3x

aH

== aH aF af(f3y)

+

aH aq ag(f3y)

=

o.

af3y

aF a f(f3y) af3y

ac ag(f3y) af3y

From eq. (32) it follows:

and

ac =R

ag(f3y)

and eq. (34a-b) ean he simplified to

aH

=

Faf(f3x)

+

Cag(f3x)

=

0

af3x

af3x

af3x

aH

=

R [F 8f(f3y)

+

c ag(f3y)]

=

o.

af3y

af3y

af3y

(34a)

. (34b)

(35a)

(35b)

Since F and Gare function both of

f3x

and

f3y

and are present in both the eqs (35a-b), the

following are neeessary concli tions for satisfying eq s (35a-h ):

af(f3

x

)

~

af(f3y)

af3x -

af3y

·15

og({3x) _ og({3y)

o{3x o{3y' (36b)

Equations (36a-b) show that all the local extrema, except the zeroes of H, lie on the diagonal f3x = f3y.

For the case of lst order upwind differencing, eq. (36) can be very simply written as

sin(f3x) = sin({3y) cos({3x)

=

cos(f3y).

On the diagonal {3x = (3y the function

IIP

n[z({3x,{3y,

R),

v] 11 is equivalent to the 1D function

IIP

n[z({3),(l

+

R)v]!I, where {3

=

{3x or {3

=

{3y indifferently.

Therefore the number of possible local maxima.is the same that for the 1D case :

N=m+l.

Also here the low frequency maximum e({3x .= 0, (3y = 0) has to be substituted by the maximum va.lue of

e

on the boundary between low frequency and high frequency domain, see Fig. 12b.

Therefore the minmax problem (33) can he-solved by applying the same technique proposed in sections 2.4 and 2.5; the low frequency boundary maximum as well as the local maxima are evaluated by a 1D numerical procedure, due to the properties of the 2D multi-stage schemes.

5. TWO DIMENSIONAL RESULTS

. The three finite differencing operators mentioned In chapter 3 have been used to

discretize the 2D advection equation, eq. (25).

5.1 Optima} parameters for a selection of schemes

Theoretical smoothing properties of the 2D multi-stage schemes have been optimized by evaluating the parameters Ck, k = 1, ... , n - 1 and v for 'whlch u is minimized. In eq.

(33), describing the 2D minmax problem, the coefficient R depends on the convection angle

(J. Therefore also the optimal parameters Ck and v are func~ion of (J.

As a first case, the fust order upwind 2-stage seheme has been considered. For such a seheme there is one maximum on the low frequency boundary and an ot her one is located at the point

({3x

=

"Ir,

(3y

=

"Ir)

for symmetry, sirnilar to the lD case. The isolines of the magnitude of the amplification factor (J are shown in figs. 13.a,b,c,d for different values of

the convection angle (J.

In figures 14a-b the optimal paraqleters Cl (equal to a) and v are plotted versus the convection angle, while the smoothing factor distribution is shown in fig. 14c. It also illustrates the wen known fact that for grid aligned advection no damping can be provided for waves aligned with the norm al direction, since no mechanism for cross diffusion is present in this case.

In fact since the surface

e({3x, (3y)

can be simply obtained by shifting a 1D curve

(e

versus

(3x)

along the {3y-axis, the value of

e

on the (3y-axis is

e({3x

=

0,

(3y)

=

e({3x

=

0, (3y

=

0) = 1. The smoothing factor is therefore u

=

1 for (} = 0°, whatever coefficients are used. For this reason and since all the maxima are imposed to be equal,

e({3x, (3y)

is r epre-sented by a constant field

e

=

1 for (J'-+ 0°; fig. 13d shows'in fact that

e({3x,{3y)

is already close to the uni ty in the entire domain: for (} = 1

° .

The 2D equation (25)is reducedto tl~e,ID equatio~ (1), when the advection is grid aligned ((J

=

0°). Therefore the lD optima! parameters are used, instead of the 2D

parameters, for R=O. _':

The same considerations are valid for the optimal parameters and smoothing factor

of a second order fully upwind 2-stage scheme (K = -1), plotted in figs. 15a,b,c. Like in the 1 D case, v is very low for this scheme. The isolines of f! for (0 = 30°) are shown, as an

.,

example, in fig. 15d.

Again, similar plots are presented in figs. 16a,b,c,d for the third order upwind biased 2-stage scheme (K =

k).

Also the 4-stage scheme using the same three different spatial differencing operators mentioned before, has been analyzed. Here the time step, proportional to v, and the three

.'

coefficients Ck of eq. (12b) (cn = 1) have to be calculated as a function of O.

The results are shown in figs. 17a,b,c,d, lSa,b,c,d·and 19a,b,c,d, respectively for the first order upwind, the second order fully upwind (K = -1) and the third order upwind biased (K =

t)

4-stage scheme.

The number of maxima is Nmax. = 3: two of them are local.extrema, while th,e other

one is a boundary maximum, as is shown by figures 17d, 18d and 19d.

, .

In all the cases the dimensionless time step v is of course strongly dependent of 0,

while the coefficients Ck are almost constant. The sensitivity of the smoothing factor to

the coefficients Ck is low, so th at constant values of Ck could be used, almost without

increasing the smoothing factor. Anyway, the practical application of the multi-stage schemes requires the use of the products q(O)v(O), k = 1, n, so that no advantage is obtained by using constant q. In Ta.bles 2-4 an interpolation function is proposed for the

three space differencing operators and a number of stages going from 1 or 2 to 4.

The results presented in the next section are obtained by using a discrete distribution

of the optimal coefficients : the products Ck(O)V(O) are calculated for 0 = k6.0, k=1, ... ,45 and 6.0 = 1 ° and stored in a matrix. For grid aligned advection (0 = 0°) the one di-mensional optimal coefficients are used. In this way few storage is required and only the

convection angle (} (in degrees) is required at each computational cello

5.2 Application to scalar

2D

test cases

The fust order upwind and the second order fully upwind

(K,

= -1) 4-stage optimal

schemes have been first tested by means of single grid scalar computations, in order to check that the real smoothing properties are improved with respect to the classic al Runge-Kutta schemes.

Local time step is generally used in the literature, due to the different speed at which

each wave travels and because the CFL number 1I is normally fixed independentlyon the

convection anglej here not only the time step is different from point to point, but also the CFL number, depending on the value of the convection angle in each cello For the same

reason also different parameters Ck are used for each point.

5.2.1 Shear problem

The first test case is the classical scalar shear problem, where the convection speed

in eq. (25) is constant and equal to

ä=

(a,b)

=

[cos(a-),sin(a-)]. The convection angle

a-is fixed at a- = 22.5°, in tha-is way the advection being as less as possible aligned with the Cartesian grid.

Figures 20a-b respectively show the convergence history of the optimal first order

upwind 4-stage scheme (curve 2), compared with a classic al Runge-Kutta 4-stage scheme

(Ck = n-~+l' and maximum stabie 1I), (curve 1), and the isolines of the solution, for a

80·80 grid.

The flat part of both the convergence histories is due to the convection and strongly

depends on the CFL numberj optimal multi-stage schemes have a greater CFL number, so that this range is reduced. Anyway, when multigrid is applied, the flat part width, which also depends on the initial solution, here taken as u=O everywhere, is normally strongly reduced, due to the increased convection speed for the low frequencies, as a result of the use of coarse grids.

The second part of each curve is instead strongly domi,nated by damping of the high frequency error, introduced by the inevitable parabolic behaviour due to the seheme. The asymptotic convergence rate of the optimal 4-stage seheme, characterizing its smoothing properties, is about twice than the one obtained by using the Runge-Kutta scheme.

The second order

(y;;

=

-1) 4-stage scheme has then been applied to the same test case, without limiters, for a 40·40 mesh. The same conclusions as for the previous scheme are obtained for the convergence histories shown in fig. 21a. The isolines of the solution, of course not monotonic, are plotted in fig. 21 b.

5.2.2 Rotational flow problem

Another linear test case which has been considered is the rotational flow problem, in which the convection speed is no longer constant,. as in the previous case, but is function of both x and y . . The vector ä is shown in figure 22. The solution obtained by using first order upwind differencing is plotted in figure 23b as a y-section: the discontinuous boundary condition on the left side of the figure is propagating along concentric circles (see Fig. 23c for the isolines of the solution) and the step is smoothed because of the artificial viscosity introduced by the scheme (see right side of Fig. 23b).

The convergence history is shown in figure 23a (curve 2) and compared with the one obtained by using classic al Runge-Kutta coefficients (curve 1). A 72·72 grid has been used. The asymptotic convergence rate is strongly increased by using the optimal coefficients.

Also in this case, the second order fully upwind 4-stage scheme of course generates some overshoots, as no limiters have been introduced, see figures 24b and 24c. The con-vergence histories are very similar to those corresponding to first order differencing and are plotted in figure 24a.

5.2.3 Burgers equation

J

The boundary conditions and the analytical solution are shown in fig. 25 : the character-istics are converging and create a normal shock.

The numeri cal solutions respectively obtained by using first and second order fully upwind differencing with a 80·80 mesh, are shown in figures 26b and 27b; the convergence histories are plotted in figures 26a and 27a and denoted by the number 2 for the optimal schemes and by 1 for the classical Runge-Kutta coefficients. The asymptotic convergence rate is still increased, but the gain is lower than for the ot her test cases.

The sol ut ion of the oblique shock Burgers problem is proposed in the chapter 8, this case being tested directly by a multigrid code.

6. APPROXIMATE OPTIMIZATION STRATEGIES

The method proposed in this report solves the problem of ,the minimization',of the smoothing factor of a multi-stage scheme, for any spatial finite differencing operator. In the following part this will be referred to as full 2D optimization.

The numeri cal strategy requires that all the values (ä,

iJ)

investigated during the iterative procedure satisfy the condition '

Nmax

=

iJ ,

( i)

which complicates the practical implementation of the sol ver.

Moreover, as already said in chapter 5, the results show that for grid aligned advection the magnitude of the amplification factor is represented by a constant field {!

=

1. When the convection angle is reducing to B = 0°, the boundary maximum is in fact located closer and closer to the point

(/3x,/3

y ) =

(0,

f)

where the magnitude of the amplification factor

is {! = 1, for B = 0°. Since the maxima are imposed to be equal, the whole field becomes

constant and equal to the maximum value (7.

The problem has been simply overcome by using the 1D optimal coefficients for B = OOi in this way the functions

v(B)

and

Ck(B)

present of course a discontinuity near the value

B = 0° i anyway this has not created any problems during the numeri cal solution of the last two test cases, for which the advection is somewhere grid aligned.

Approximate but faster strategies can be used to calculate coefficients which reduce the smoothing factor, although without being able to have a real control on it.

The procedure described in chapter 4 can in fact always been applied to a part i al, ' optimization, if the condition (i) is satisfied. The basic idea is to apply the same theory on some maxima defined in different domains and then to calculate the resulting smoothing factor, as it will be shown in the following sections.

6.1 Partial two dimensional optimization

The previous considerations about grid alignment suggest a first partial strategy: the boundary maximum is supposed to be located on the the Hnes f3x

=

f

and f3x

=

-f,

see figure 28a, instead of being defined on the entire border between the low and high frequency domain. The local maxima are always defined as 2D extrema and therefore they are situated on the diagonal f3x

=

f3y. In this way the 1D optimal coefficients are automatically obtained when 8 = 00.

6.2 One dimensional optimization along the convection direction

A different approach is to define both the boundary and the local maxima in a 1D domain, represented by the line f3y = qf3x. If q = tan(8), this line is coincident with the convection direction, see fig. 28b; the schemes are therefore optimized with regard to Fourier modes propagating in the direct ion of the convection speed. Of course, every

numerical scheme introduces error modes belonging to the entire 2D spectrum. For this reason a full 2D analysis has been proposed in chapter 4.

The strategy used in [3] cannot be applied, although we are dealing with a 1D domain, because (] is still a 2D function.

The 1D coefficients are again obtained wh en 8 = 00.

6.3 One dimensional optimization along the 450

direction Using the previous approach with q = 1, namely 8 = 450

, both the boundary and the

local maxima are defined on the diagonal f3y = f3x, see fig. 28c. The strategy is therefore

similar to the fuIl 2D optimization, except for the definition of the boundary maximum and it is referred to as the partial optimization along the 450

direction.

A st rong simplification is obtained : since f3x = _{f3 y, (] is a 1D function, as explained}

in chapter 4; therefore only a 1D analysis is required and, moreover, the coefficients Ck are

constant for the entire range of 8 values, and equal to the 1D coefficients. The sum of the

optimal Vx and v y is also constant and allows to write Vx as

· VlD

Vx

=

(1

+

R)

Of course the smoothing factor is still function of (J, as the real low frequency boundary

maximum, coincident with a, is not a lD function and is therefore dependent on

R

or (J.

In this way, the optimal coefficients calculated in [3] have been used by Tai [14] for

2D calculations.

6.4 A theoretical comparison

The smoothing factor distributions for the first order upwind 2-stage scheme, using

the four different strategies, are plotted in fig. 29. Of course partial optimizations always

lead to higher values of a.

A wider investigation has been done for the partialID optimization along the convec-tion direcconvec-tion and has shown that it is not applicable to second order multi-stage schemes:

regions of instabilities are in fact found for a 2-stage scheme, while for a4-stage scheme

there is a wide range in which the smoothing factor is not minimized at all. This is due to

the fact that the smoothing factor is not controlled by the partial optimization. ' ,.'

Therefore the strategy proposed in section 6.2 will not be considered anymore. Also the approximate method proposed in section 6.1 is excluded, because its feature was only to obtain automatically the lD coefficients in case of grid alignment, which is not really needed.

The strategy of section 6.3 can be interesting, as it allows to use constant coefficients

Ck and a simple analytical function for V_{x ,} both coming from the simpler lD analysis,

although the smoothing factor distribution is much worse than the one obtained by using the full 2D optimization, see fig. 29.

7. MULTIGRID TECHNIQUE 7.1 Elliptic case

Multigrid fast solvers have historically been studied first for the solution of

steady-state elliptic problems. For such cases it has been observed that a classical re1axation method rapidly damps the high frequency components of the error distribution, so that

after few iterations only the low frequency modes are present. Due to the bad low frequency damping properties of classical relaxation solvers, the convergenee rate is strongly reduced.

On the other side the absence of high frequeney components implies that the error distribution is smooth and can therefore be transferred from the fine to a coarser grid by

a collection operator

fih,

where hand 2h respectively denote the fine and the coarse grid mesh size.

For an elliptic problem this results in damping the low frequency components at a lower computational cost: these modes are in fact transferred from the fine to the coarser grid, where they look like high frequency oscillations, due to the doubled mesh size, see

figure 30. In this way they are reduced by a cheaper relaxation on the coarser grid. The coarse grid correct ion of the error is then interpolated back to the fine grid by the operator

I;h'

and the solution is then updated. A basic requirement is the use of an efficient smoother, namely a relaxation seheme capable to damp efficiently the high frequency error components.

Different multigrid cycles ean be designed; they will be described in detail in section 7.3.

7.2 Hyperbolic case

A 2D hyperbolic problem, as the set of Euler equations, is eharacterized by the con-vection of some physica.l properties along particular paths called characteristic surfaces.

The same wave propagation is typical of the simpier scalar model, as has been observed in section 5.2; there, it has been underlined that the disturbances are more rapidly expelled

through the outer boundary if a greater time step. is used. Therefore CFL number has to be maximized for a single grid calculation, of course within the stability region.

When a multigrid algorithm is applied to an hyperbolic problem, the speed of low frequency waves is increased because of the doubled mesh size, if CFL number is taken constant over all the grid levels [1,2].

On the contrary, the high frequency modes are eliminated by damping, since artificial viscosity introduced by the scheme implies a parabolic-like behaviour. Therefore most of

the multigrid theory can be. extended to the hyperbolic case, where the use of a good

smoother on the coarser grid is again a basic requirement for the good performance of the multigrid cycle. An efficient smoother is also needed on the top level in order to smooth the fine grid error distribution.

The smoothers can be divided in two classes, the first one including implicit

Gauss-Seidel type schemes (P J, LJ, PGS, LGS, Ze1:>ra, etc.) and the second one being represented

by the explicit schemes. Of course a classical forward Euler time integration scheme is not a good smoother if applied only on ce for each relaxation step. Therefore only the multi-stage time integration schemes and a sequence of applications of the Forward Euler scheme are included in this class.

Classica! Runge-Kutta coefficients are chosen to obtain higher order time accuracYi

for a steady state calculation this is not required and the efficiency of explicit schemes as

smoothers can be improved by using the coefficients calculated in chapter 5.

7.3 FMG and FAS cycle [7,10,13]

The 2D advection equation can be rewritten as

du

dt

+

res(u)

= 0, (37)

where res( u) is the residual, namely the steady-state solution satisfies

res( u)

=

O. Let us

suppose that four levels with equa~ mesh spacing both in x- and y-direction (~x = ~y = h)

The Full Multigrid (FMG) cycle is used to obtain an initial estimate of the solution on the fine grid. Starting from a crude initial guess on the level 1 a FAS cycle is performed from station A to station Band the solution is then bilinearly interpolated up to the finer grid 2, where the scheme is iterated until the finest level is reached.

The FAS cycle at level

e

consists in the following steps:

(1) an estimated first guess or the solution obtained from the FMG cycle is used;

(2) perform p(

e)

pre-relaxations, namely apply the multi-stage scheme p times; (3) calculate the defect dhl = -res( Uht);

(4) interpolate the solution from level

e

to level

e

-1, or use the solution from the previous iteration; in the present tests the latter strategy is shown to be better;

(5) transfer the defect from the fine to the coarser grid by a restriction operator Il h

dht _1

=

Il h dht

=

-Ilh res(uht);

(6) perform a FAS iterations on the grid ht-l solving

where P( Uht_l) = res( Uh t_1) - Ilh res( Uht); the m-th step of the multi-stage scheme results in

U(m) = u(O) _{_ cm6t (res(u(m-l) _ P(u(O) ));}

ht-l ht-l ht_l ht-l

in this way the evolution on the coarse grid is driven by the residual on the fine grid in the first stage. For a = 1 it results in a V -cycle, while for a = 2 it is a W-cycle; in

the present work only the V-cycle is used;

(7) interpolate the correct ion 6Uht_1 = Uh t_1 - Uh~~l and update the solution on level

e:

Uh t = Uh~d

+

I;h 6Uht_l;

(8) improve the solution by means of q( e) post-applications of the multi-stage scheme. Of course the FAS cycle on the coarsest grid does not include steps (3)-(7). In the

present p(e)=q(e)=1 except on the coarsest grid, where the exact solution is required; moreover the post-relaxation on the finest grid has been skipped, resulting in a less

expen-sive cycle.

The restriction and the prolongation operators should satisfy the foll?wing Galerkin relations that guarantee the elimination of fine grid low frequencies [9], [12] :

(38) and

(39) where Arepresents the spatial differencing operator.

Unfortunately, in case of finite volume formulation (cfr section 7.5.1), eqs. (38) and

(39) can be satisfied only by the 1st order scheme, resulting in the following choi~e for the

prolongation and the restriction operators:

and

For the second order scheme this choice leads to the use of a 1st order operator on the

coarser grids, [13]. Therefore, the FAS cycle is used in the present tests both for the first

order upwind and the second order fully upwind scheme.

A FAS Multigrid Euler solver with relaxation smoother, developed at CWI,

Amster-dam, in a joint BritejEuram project

[6],

was adapted for the present multi-stage smoothers

and scalar equation.

7.4 Saw-Toot h cycle

The impossibility of satisfying the second Galerkin relation, eq. (39), for a K-scheme

with K =1= -1 leads to the use of a different cycle, namely the Saw-Tooth cycle, which shows

empirically to perform better than FAS cycle [3].

The Saw-Tooth cycle is shown schematically in figure 32 : a first solution is obtained from an FMG cycle; steps (1-5) of the FAS cycle are here repeated, but the solution on the

coarse grid is collected from the fine grid. A multi-stage relaxation is perforrned at each level, resulting in a distribution of /luht , where f takes all values from the bot torn to the

top level. Af ter the coarsest grid is reached, all the corrections are bilinearly interpolated back to the finest grid, where the solution is updated. Such a definition of the prolongation operator is imposed by the fact that no post-relaxations are performedj of course it implies that also the first Galerkin relation, eq. (38), is not satisfied.

The new operators needed by the Saw-Tooth cycle have been implemented in the existing code.

7.5 Obligue shock test case - numerical formulation

The oblique shock testcase is obtained by solving Burgers' equation on a unit square with boundary conditions u(x=O,y)=1.5, u(x=1,y)=-O.5, u(x,y=O)=1.5-2x, cfr. fig. 33 .

. ,

7.5.1 Space discretization

The discretization of the nonlinear 2D advection equation

~~

+

8~~u)

+

8

g

}yu), = 0 is

the classical flnite volume approach with cell centered unknowns, cfr. figure 34. For cell i,j on the Cartesian mesh it is given by

( au) /lx /ly =

(f~+.1

. -

f~_.1

.) /ly

+

(g* '+.1 -

g~

'_.1)

/lx. (42)

at

"

t 2 ,} ' 2 ,} ',} 2 ',} 2

',}

Here, r+.1 ' and g~ '+.1 are the numerical flux functions in the x-, respectively y-direction,

, 2 ,} ',} 2 given by: f*+l '

=

f+

(u~+.1

.)

+

f- (u-:-+.1 .) , 2 ,} , 2 ,} I 2 ,} (43a) g* '+.1

=

g+ (u+ '+.1)

+

g- (u-:-'+.1) I,} 2 ',} 2 ',} 2 (43b)

where u+ and U- are MUSCL extrapolations allowing to recover higher order schemes:

1+1Ç

1-1Ç

u~+.1 .

=

Ui

+

- 4 - (Ui+l - Ui)

+

- 4 - (Ui - ui-d

, 2')

(44a) (44b)

and similar expressions for u~ ',} '+.1 and 2 U-:-I,} '+.1' 2

For K

=

-1 the scheme is second order accurate; For K

=

i-

third order accuracy can

be obtained only in the lD case; nevertheless it will be referred to as the 3rd order scheme in the following.

j+ and j- are the upwind contributions from the left, respectively right, defined

according to the convective speeds (u in the x-direction, 1 in the y-direction) .:

j+( u)

=

j( u), j+(u) =0, j-(u)

=

Jeu), j-(u)=o, for for for for g+(u) = g(u) g-(u)=O.

7.5.2 Definition of cell convection speed

u~O

u<O

u~O

u>O

The optimal coefficients have been calculated by using a linear von Neumann analysis.

In the case of a non linear equation a proper linearization is required; let us c(:msider only

the flux along x-direction : if the convection speeds on the walls

(i

+

~,j), defined as

f.~~ .-f~~ .

a· ~ . - • 2·J • 2.J , and (i - -21,

j)

have the same sign (suppose a

>

0), the MUSCL 1+2'} - u -u+

i+~,j i+~,j

extrapolation allows to recover the classical fini te difference approximations of ~: by simply writing

au

ax

-+

+

u.+~ .-u. ~ . 1 2'} 1 -2 ,} ~x

where the resulting expression depends on the choice of K.

Therefore the flux difference is linearized as

j + .+~

.

-

j+ ._~ . -

-

ai,j

(+

u.+~ . - u._~

+)

.

I 2 ,J I 2 ,J I 2 ,J I 2 ,J

and the cell convection speed is consequently defined as

jt+~

J' - j'!-i J'

This new definition of convection speed does not require the evaluation of other quantities and it is crucial wh en the solution is still far from convergence.

Equation (45) shows that if the convection speed is nihil, also the flux difference is nihil and no limitation on the time step is required.

For the first order upwind scheme, the cell convection speed coincides with the con-vection speed evaluated at the inlet wall.

For th~ Burgers equation an analytical expression can be given: U-:'-+l _{1 2']} .

+

u-:'-

_{1 - 2 ,]} 1 .

ai,j

=

2

8. MULTIGRID RESULTS .

All the calculations are obtained by using a 64*64

fine

gria and bot tom level grid of

2*2 or 4*4, respectively fOLthe FAS and the Saw-Tooth cyclej one work unit corresponds to one residual calculation on the fine grid, which is equivalent to one Euler explicit time

stepping on the fine mesh. The absolute residual is converged to double precision machine

zero.

The test case is first solved with first order spatial accuracy. The isolines of the solution are shown in figure 35 and a cut crossing the shock at y=0.75 shows that monotonicity is

of course satisfied, cfr. figure 36; this will be no longer true for the higher order schemes

as no TVD limiters have been used. The cut at y=0.25 is presented in figure 37.

The 2- and 3-stage optimized schemes show to perform the best. In figure 38 the con-vergence history lising a 2-stage optimized scheme is compared with ot her explicit schemes. The convergence history is represented by a straight line, showing that multigrid combines weU with the new smoother; the convergence rate is always equal to the asymptotic con-vergence rate of the single grid calculation. The fastest concon-vergence for a sequence of non optimized Forward Euler integration schemes is obtained for p(f)=q(f)=4; the line shows that af ter few iterations the convergence rate slows down, although also in this case few work units are required to have machine zero absolute residual; for many practical Euler calculations, however, Runge-Kutta schemes are preferred instead of a sequence of For-ward Euler schemes, because of stability problems, so that the comparison with a classical Runge-Kutta 4-stage scheme is proposed; the improvement is very clear with a speed up

of a factor 4 for the optimized multi-stage scheme. The comparison with a single grid

calculation shows that the first flat part has been eliminated by the Multigrid-Multi-stage

combination. Finally the comparison with the 2-stage scheme, optimized by using the approximate strategy proposed in 6.3 (optimization along 45° degree direction) with

con-stant vxis shown in the same figure and is in between the Runge-Kutta and the sequence

of Forward Euler schemes.

history for three Multigrid calculations with different number of nodes on the fine grid (32*32, 64*64, 96*96) and same number of stages are compared. The reduction factor,

.1.

defined as p =

(\~::~:~:::D

n is 0.43 for the 32*32 grid, 0.51 for the 64*64 grid and 0.53

for the 96*96 grid. The reduction factor is increased by a factor 1.18 when going from the 32*32 to the 64*64 grid, and only by a factor 1.04 when going from the 64*64 to the 96*96 grid. As the curves are almost superposed, the amount of work needed for converging the solution to machine zero is proportional to the number of grid points, showing correct multigrid performance.

The 32*32 and 64*64 single grid convergence histories are shown only for visual com-parison. In fact it results that the number of work units needed by a single grid to reach convergence is strongly dependent on the number of points, so that a quantitative com-parison has no significative meaning.

For a seconcl order scheme

(K,

= -1) the solution is more accurate but no longer monotonie, as is shown in figs. 40, 41 and 42; the convergence history using the optimized 3-stage scheme, in figure 43, is no more represented by a straight line, but the flat part is strongly reduced and the final slope of the curve is aga.in similar to the one obtained by means of a single grid calculation.

A sequence of three applications of forward Euler scheme, as weU as the 3-stage Runge-Kutta scheme do not perform as good: the lower slope indicates that these schemes introduce new low frequency errors in the fine grid, as the smoothing is not efficient. Again it should be underlined that the comparison with a sequence of Forward Euler schemes is possible only because a scalar equation is solved; for Euler calculations, literature shows that even some stability problems occur in some complex flow situations.

Finally the solution obtained by using the

K,

=

!

scheme is plotted in figs. 44, 45 and 46. The solution in the smooth part is a little more accurate than the second order fuUy upwind solution, but strongly non monotonie in the shock, as shown in figure 45. The performance of the multigrid cyele is probably influenced by the strong overshoots; moreover, the lack in satisfying Galerkin properties influences the multigrid performance.

As a consequence the coarse grid corection is not very efficient. The convergence rate is

almost constant, but the slope of the multigrid curve is lower than the single grid one.

It is believed that the use of an appropriate limiter for eliminating the overshoots would

increase the convergence ratej the implementation of a limiter will be done in the near

future.

Nevertheless, the optimization allows to get an improvement of about 2.5 with respect to the convergence rate obtained by using Runge-Kutta coefficients.

The comparison with a forward Euler scheme is not possible, since Euler explicit is

9. CONCL USIONS

A numeri cal technique for optimizing smoothing properties of an explicit multistage scheme for a general sp ace discretization of a scalar convection equation has been

devel-oped. The new formulation allows to recover optimized multistage schemes both for one

and two dimensionsj in the ID case it produces equal or bet ter results than the Van Leer-Tai-Powell optimization. The extension to three dimensions should be straightforward, although with a substantial increase of the computational cost of the procedure.

The optimized 2D schemes have been tested as smoothers in a multigrid code, sol ving a scalar non linear test case. The results showastrong improvement with respect to the classical Runge-Kutta schemes. Although a Forward Euler scheme (when possible) shows to be more competitive than a classical Runge-Kutta, in the literature it is believed that

multi-stage schemes are more suitable than a single-stage scheme for an Euler calculation,

especially in a multigrid code, because of stability problems.

A further investigation for the set of Euler equations is needed and will be carried out

in the future.

The need of an appropriate limiter for the higher order schemes has been pointed out.