Citation
Eulerian-Lagrangian localized adjoint methods (ELLAM)

Material Information

Title:
Eulerian-Lagrangian localized adjoint methods (ELLAM) discretization strategies for advection-dominated flows
Creator:
Trujillo, Rick Vance
Publication Date:
Language:
English
Physical Description:
40 leaves : ; 29 cm

Subjects

Subjects / Keywords:
Fluid dynamics ( lcsh )
Lagrange equations ( lcsh )
Algorithms ( lcsh )
Boundary value problems ( lcsh )
Algorithms ( fast )
Boundary value problems ( fast )
Fluid dynamics ( fast )
Lagrange equations ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references.
General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Rick Vance Trujillo.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
25334816 ( OCLC )
ocm25334816
Classification:
LD1190.L622 1990m .T78 ( lcc )

Full Text
EULERIAN-LAGRANGIAN LOCALIZED ADJOINT METHODS
(ELLAM)
DISCRETIZATION STRATEGIES FOR
ADVECTION-DOMINATED
FLOWS
by
Rick Vance Trujillo
B.S., Metropolitan State College of Denver, 1988
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Mathematics


This thesis for the Master of Science degree
by
Rick Vance Trujillo
has been approved for the
Department of Mathematics
by
Thomas F. Russell
Stephen F. McCormick


Contents
1 INTRODUCTION 1
2 BASIC CONCEPTS 6
3 EXACT WEAK FORMULATION 10
4 VANISHING ADJOINT AND TRACKING ALGORITHM 15
5 APPROXIMATE WEAK FORMULATION 20
6 INLET BOUNDARY DISCRETIZATION 23
7 OUTLET BOUNDARY DISCRETIZATION 28
8 MEASURE ADJUSTMENT AND IMPLICIT ELLAM 31
9 FORWARD TRACKING 36
10 SUMMARY 38
References
40


Rick Vance Trujillo (M.S., Applied Mathematics)
Eulerian-Lagrangian Localized Adjoint Methods (ELLAM):
Discretization Strategies for Advection-Dominated Flows
Thesis directed by Professor Thomas F. Russell
Eulerian-Lagrangian Localized Adjoint Methods (ELLAM) are generalizations
of the successful Eulerian-Lagrangian Method (ELM) discretization strategy for
the Finite Element Method (FEM) solution of the equations that govern advection-
dominated flows. In addition to the non-oscillatory and large-Courant-number
properties of ELM solutions, ELLAM offers the features of numerical conserva-
tion and systematic treatment of boundary conditions. The ELLAM discretization
strategy presented here is guided by the analytic solution of advection equations by
the Method of Characteristics (MOC) and is designed to treat advection-reaction-
diffusion equations with variable coefficients and in multiple spatial dimensions. In
this paper, we focus on the theoretical aspects of ELLAM that bear directly upon
numerical implementation. We explore exact and approximate weak formulations,
vanishing adjoints, tracking algorithms, measure adjustments, and boundary dis-
cretizations.
The form and content of this abstract are approved. I recommend its publication.


1 INTRODUCTION
Advection-reaction-diffusion relations arise in many fluid dynamics problems
(e.g. groundwater flows and petroleum reservoir simulations), or in related prob-
lems (e.g. semiconductor modeling). These relations pose stiff tests for numerical
methods, due to the advection or hyperbolic component of the problem. Tradi-
tional Eulerian (static-grid) techniques generally suffer from numerical diffusion
and/or spurious oscillations when confronted by advection-dominated problems;
for a collection of perspectives on this subject, see [5]. These difficulties can be
ameliorated by using very fine computational grids or, preferably, by using opti-
mal test functions (OTF) that satisfy a local adjoint relation. The OTF strategy
( see, e.g., [2] ) allows the test functions to adapt in space in order to minimize
spatial discretization error; however, they are static in time. As a result, the
range of Courant numbers that can treated by OTF methods is usually restricted
to be less than unity.
For advection-dominated flows, physics dictates that numerical methods with
Lagrangian, or moving-grid, properties that can adapt to changing domains of
influence would be more appropriate. In particular, the Eulerian-Lagrangian
Methods (ELM) [1, 4] that track the advecting flow in order to find the proper
domain of influence for the hyperbolic subproblem represent an approach that
is closer to the underlying physics of the problem. On grids whose fineness cor-
responds to an accurate representation of the moving solution front, ELM can


obtain non-oscillatory solutions that are not numerically diffused. However, ELM
suffer from two major difficulties: they fail to conserve mass and to handle general
boundary conditions.
The Eulerian-Lagrangian localized adjoint method (ELLAM) of Celia, Ewing,
Herrera, and Russell [3, 6] combines ideas from OTF and ELM by using space-
time test functions that satisfy a local adjoint relation. In this setting, the test
functions are defined on space-time finite elements that track the evolution of
an advecting flow between static computational and data grids at different time
levels. Furthermore, ELLAM is able to conserve mass and systematically treat
general boundary conditions. Previous papers [3, 6] have presented these ideas
for advection-diffusion equations with constant coefficients in one spatial dimen-
sion. For the case of advection-diffusion equations with variable coefficients and
in multiple spatial dimensions, Russell and the author have recently presented
what we shall refer to as the implicit ELLAM formulation [7]. This new formu-
lation represents a significant break from ELM and earlier ELLAM formulations
in style and substance. The implicit ELLAM formulation views the advection-
diffusion problem in Eulerian coordinates and only applies Lagrangian concepts
to the formal adjoint of the advection term. In this way, the incorporation of
boundary conditions and mass conservation, which distinguishes ELLAM bom
ELM, is straightforward. However, the treatment of the adjoint implictly induces
a Lagrangian coordinate system which is at the heart of the entire method. Thus,
2


we are compelled to find the appropriate ways to introduce Lagrangian concepts
for ELLAM in a variable coefficient and multiple spatial dimension setting.
The goal of this thesis is to develop methods for incorporating Lagrangian
ideas in an ELLAM setting. Our approach will explicitly use Lagrangian co-
ordinates and view the entire ELLAM process from the perspective of a La-
grangian reference frame. Finally, we will develop techniques for converting from
Lagrangian to Eulerian coordinates, in order to link our explicit ELLAM formu-
lation to the implicit ELLAM formulation of [7]. The explicit ELLAM discretiza-
tion strategy that will be presented here is designed to treat advection-reaction-
diffusion equations with variable coefficients and in multiple spatial dimensions.
In Section 2: Basic Concepts, we will pose the non-divergence and divergence
forms of the model problems for advection-reaction-diffusion processes. Also,
we will see how the analytic solution of the advection-reaction equation by the
Method of Characteristics (MOC) introduces us to Lagrangian concepts in a very
natural way. Furthermore, we will use our knowledge of the MOC to guide us
through the entire ELLAM formulation.
In Section 3: Exact Weak Formulation, we will discover how to imbed our
knowledge of the analytic solution of the advection-reaction equation into the
weak forms of our model problems. We define a space-time, Lagrangian domain
in terms of backtracking, from a computational grid to a data grid, along the
characteristics of the analytic solution. In this manner, we can view ELLAM
3


as a space-time, Lagrangian variant of the Finite Element Method (FEM) for
non-advective parabolic equations.
In Section 4: Vanishing Adjoint and Tracking Algorithm, we directly en-
counter the inner workings of the Lagrangian coordinate system that is created
by the vanishing adjoint. With an eye towards numerical implementation, we de-
velop a simple, predicted tracking algorithm that is based on the exact solution
of an approximate adjoint relation. And with an eye toward solution accuracy,
we develop a higher-order, corrected tracking algorithm that is based on the
approximate solution of the exact adjoint relation.
In Section 5: Approximate Weak Formulation, we modify the results of Sec-
tion 3, in order to incorporate the developments of Section 4. Although the pre-
dicted tracking algorithm defines an approximate Lagrangian coordinate system
for our equation framework, we can incorporate the corrected tracking algorithm,
which is more faithful to the exact Lagrangian coordinate system, to characterize
the solution.
In Sections 6 ( Inlet Boundary Discretization ) and 7 ( Outlet Boundary Dis-
cretization ), we apply the framework of Section 5 and systematically incorporate
boundary conditions into our explicit ELLAM formulation. Since we are working
in a space-time domain, we are led to incorporating time as a boundary variable.
The resulting measure adjustment ( from space to time ) to the integrals of our
weak formulation is an Eulerian concept and represents a significant development.
4


In Section 8: Measure Adjustment and Implicit ELLAM, we formally de-
velop measure adjustment of the integrals of our weak formulation as a means
to convert from Lagrangian to Eulerian coordinates. With measure adjustment,
the explicit ELLAM formulation conserves mass and, coupled with the results of
Sections 6 and 7, distinguishes itself from its ELM predecessors. Furthermore,
measure adjustment leads us to developing the analogous concept of reference
frame adjustment to Lagrangian coordinates for implicit ELLAM. Thus, we are
able to construct a gateway between implicit and explicit ELLAM formulations.
After performing measure adjustment, we find that the implicit and explicit for-
mulations are equivalent to within O(At)2 discretization error and lead us to the
same discrete equations.
In Section 9: Forward Tracking, we identify the forward tracking relation
that is derived from the predicted backtracking relation of Section 4 as the mech-
anism behind measure and reference frame adjustment. Employing the often-used
piecewise linear Raviart-Thomas representation of the velocity field, we outline
the numerical implementation of forward tracking. As a result, we find that ex-
act forward tracking and mass conservation is possible, whenever we avoid shocks
and the tracking of characteristics is well-defined.
5


2 BASIC CONCEPTS
The ELLAM framework is designed to treat the two canonical forms of
advection-reaction-diffusion equations in m spatial dimensions:
Non-Divergence Form
Cal = Au + Hu + Vu /(x, t) (1)
where
Au = ut + t/(x,i) Vu
Hu = k(x,t)u
Du = V (D(x,f)Vu)
are the advection, reaction, and diffusion components of (1) respectively, and
v(x, t) = (i(x, £),..., um(x, £)), fc(x, £), and D{x, t) are the velocity field, reaction
and diffusion coefficients, respectively, where x = (xi,... ,xTO); and
Divergence Form
Qu = Au + H*u + Du = /(x, £) (2)
where
H*u = Hu + (V u(x, t))u
is a modified reaction term that contains an artificial reaction induced by the di-
vergence of the velocity field, v(x, £), in the Lagrangian reference frame. In a later
section on measure adjustment, we will see that the velocity dependent reaction


is canceled when we convert to Eulerian coordinates. Thus, our formulation for
the Divergence ( or Conservation law ) form will conserve solution mass and is
consistent with physical law.
In relation to previous work [7], we note that the usual statement of the di-
vergence form is (.4* + 71 + V)u = /, where A* = ut + V (u). However,
methods that apply Lagrangian concepts explictly need to work with the La-
grangian derivative .A(-) instead of Therefore, explicit ELLAM splits A*
into A and TV, in order to obtain a consistent formulation for (1) and (2). In
the Lagrangian setting of explicit ELLAM, the adjoint of the advection term for
(1) and (2), A, is A, so Lagrangian ideas can be used throughout the formu-
lation. In the Eulerian setting of implicit ELLAM, the adjoint of the divergence
advection term, A*, is A, so Lagrangian ideas can be applied to the adjoint.
However, in Eulerian coordinates, the adjoint of the non-divergence term, A, is
A* and there is some ambiguity in the treatment of the adjoint.
To complete our specification of the problem, we note that equations (1) and
(2) hold for x £ E C Rm, t J = (t0,T] and with Dirichlet or flux boundary
conditions,
u(z,£) = h(x,t) or (yu DVu)(z,t) n(x,t) = h(x,t) (3)
for x £ 3E, t £ J and an initial condition,
u(x,t0) = w0(s), for x G E.
The ELLAM formulation for (1) and (2) is guided by a knowledge of the
7


analytic solution for the numerically problematic advection component of the
equations. For the time being, we will focus on the non-divergence form (1)
of the advection-reaction-diffusion equation. Taking D(x, t) = 0 and using the
Method of Characteristics (MOC), the solution to the non-diffusive case of (1) is
given by the following system of ODEs:
dt(s,r)
dr
t(s,tn+1) =
dx;(s,T)
-- =
Xi(s,tn+1) = Si, for i = l,..., m (5)
^dr ^ + r) r)U = r)r)
u(s,f*+1) = un+1(s) (6)
where tn+1 is the new time level, s = (si,..., sm) 6 E is a parametrization of
the grid points at the new time level, and un+1(s) is the solution at the new
time level. The idea is to take a point (*i,..., xm, £n+1), called the head of the
characteristic, from the computational grid and track it backwards in time along
the characteristics defined by (5) to its image (sj,... ,a^,£*), called the foot
of characteristic, on a data grid at a previous time level. For the present case,
we assume that t* = tn and our backtracked image remains in the interior of
1,
tn+1 (4)
8


£ and not on the boundary. With a knowledge of the solution at the old time
level, u(xj,..., x*m,tn) = un(x*), the solution to (6) gives us an expression for
the solution at the new time level, u(zi,... xm,£n+1) = ttn+1(x). Given that (4)
tells us that t = r, equation (6) tells us that
/ -*+1
un+1(x) = un(x*)exp I y fc(x(s,r),r) dr
ftn+1 ( ft**1 \
+ Jtn exp[~JT k{x{s,T,),T,)dr'J /(x(s,r),r)dr. (7)
Now, let At = tn+1 tn. After making 1-point forward and backward Euler
approximations to the first and second integrals of (7) and 0(A£)2 approximations
of the exponential terms, we obtain
un+1(x) = un(x) (1 fcn(x*) A£) + /n+1() At + 0(A£)2 (8)
as a second order approximation to the analytic solution (7). The formulation
for the divergence form (2) proceeds along the same line with k(x, t) + V v(x, t)
replacing k(x,t) in (6), (7), and (8). Symbolically, the key to our discretiza-
tion strategy will be to obtain the weak form of (8) for the advection-reaction
component of (1) and (2).
9


3 EXACT WEAK FORMULATION
For clarity and to establish, a context for our discussion, we will consider
the case of a two-dimensional rectangular computational grid, En+1, that is dis-
cretized into rectangular subpatches, E?*1. As in standard finite element methods
(FEM), we pose the weak form of (1), namely /ni+i w Cu = /n+i / w. However, in
contrast to standard methods, fl?*1 is a space-time finite element that is formed
by backtracking the computational patch, E?*1, along the exact characteristics
given by (5) to the data grid, En. O?*1 defines a Lagrangian coordinate system
where the spatial measure on the /n?+1 integral is tied to the head of the char-
acteristic and the measure of E-1"1. Also, fly1"1 is used to denote the process
rewriting the advection operator from a transport term (ut + v Vu ) in Eulerian
coordinates to a derivative ( utau ) in Lagrangian coordinates. The test function
w refers to those local space-time test functions, Wij, whose values at tn+1 are
non-zero on E-*"1. Additionally, w is standard in space on the computational
domain, E?*1, at tn+1, but is defined in ft-1"1 as the solution of a local adjoint
relation. With the proviso that the subtleties of this Lagrangian, space-time
domain will detailed in later sections, we will concentrate on the mechanics of
creating a weak form of (8).
For simplicity, we assume that we are working on the portion of the com-
putational grid whose backtracked images remain in the interior we consider
boundary interactions in a later section. Considering each component of (1)


separately, we start with, the advection operator and integrate by parts to obtain
/ wAu = / A(uw) / uAw. (9)
Jo"*1 v ynn+i ^ J
J *1
Since .A(-) denotes differentiation along characteristics and integration over fl?*1
can be viewed as integration along these same characteristics, the first integral
on the right hand side of (9) is straightforward. For the second integral on the
right hand side of (9), we define the test function, -u;, to be the solution of the
adjoint relation, Aw = 0. Therefore,
fapi wAu = /En+1 ttn+1(a) w X-+1 w (10)
since the adjoint relation tells us that w = ton+1(z) = wn(x*).
Moving on to the reaction term, we apply a forward Euler discretization to
obtain
J^iw'R,u = At J ^ kn(x*)un(x*)w + O(At)2. (11)
Finally, the diffusion term receives the standard FEM treatment for elliptic
equations as we invoke Greens theorem to obtain
f wVu= [ f V (w DVu) + f DVu'Vw. (12)
Jn$+I v yn~+ v '
Then, we apply the divergence theorem and backward Euler discretization re-
spectively, to the right hand side integrals of (12) to obtain
wDu =
-L
Jn+l
w D Vu n
+At j Dn+1 (z) Vun+l(x) Vw + 0(At)2,
11


where n(x, t) is formed by backtracking the unit outward normal on How-
ever, the dE?*1 boundary integral vanishes, since w = 0 on 5E-1"1. In practice,
E?*1 is subdivided and w is supported on more than one element and the bound-
ary integral does not vanish on those elements. However, the boundary terms
from the sudivisions do not need to be computed or stored, since they will sum
to zero when combined to form the equations for all of E?*1. Therefore, we have
^n+i wVu = At jf ^ Dn+1(x) Vun+1(s) Vw + 0{At)2. (13)
In relation to eariler work in streamline coordinate techniques [5], the treatment
of diffusion, even in advection-dominated flows, is crucial. The backward Euler
integration puts the diffusion on the new time level where there is no distinction
between Eulerian and Lagrangian coordinates. Therefore, the formal incorpora-
tion of diffusion is straightforward. However, if the flow is rapidly converging or
diverging, as determined by the divergence of the velocity field, backward Euler
integration can make large errors unless the time step is small relative to the
time scale of the physical phenomena. Thus, the time evolution of the solution
can be used to determine the time step that is required for the non-destructive
incorporation of diffusion.
To close the system, we apply a backward Euler rule to the source term to
obtain
/E. r+\*)w + 0(Atf. (14)
12


Combining (9), (11), (13), and (14), the weak form of (1) becomes
[ un+1(x)w + At f Dn+1Vun+1-Vw
Jv?+i v ' 7e"+1
U J
= L, "(*') (1 *"(**) Ai) + Af J r" + 0(At). (15)
- Lii Lii
In relation to standard FEM, we can view (15) as a conversion of (1) and (2) into
a'non-advective parabolic problem in the characteristic reference frame,
Setting; D = 0 in (15), we obtain the desired weak form of (8)
/.,(*)= r,K(i<)(l-*:"(x')At) + rA()+0(Ai)!. (16)
"u L*J
The development for the divergence form (2) proceeds along the same line with
fc(,£)+ V v(x,t) taking the place of k(x,t).
To establish a basis for comparision, the implicit ELLAM treatment of (1) (
ignoring boundary terms ) yields
[ un+1(x) w + Atf Dn+1 Vun+1 Vw / / (V v) u
J'Lij JVij Jj+WEy
= [ un(x) (1 kn(x) At)w + Af f /Ti+1 w + O(At)2.
t£7
For the divergence form, the (V v)u integral is missing in the implicit ELLAM
formulation; however, the term reappears in the data integral of the explicit
ELLAM formulation. In Section 8, we will find that the additional (V v) u and
the are equivalent to the £-+1 integration of the data terms on the right hand
side of (15).
Therefore, we have formally captured the essence of the MOC solution of the
non-diffusive case of (1) and (2) in our weak formulation. However, the true
13


essence of the ELLAM discretization strategy can be found in the vanishing of
the adjoint relation, A.w, which in turn invokes the solution of (5), the definition
of the test function, w, and the definition of the Lagrangian reference frame,
Q?.+1.
y
14


4 VANISHING ADJOINT AND
TRACKING ALGORITHM
The key to the ELLAM formulation is the solution of the vanishing adjoint
relation, Aw = 0, which is given by the MOC and the system of ODEs defined
by (4), (5) and
dw(s,r) ^ = tun+1(s) (17)
dr
where twn+1(s) is usually a standard piecewise polynomial test function that is
defined at the mesh points of E?*1. The most familiar examples of iun+1(s) are
bilinear test functions on a 2-D rectangular mesh where u;|j+1(si,s2) = 1 at (Xi,yj)
and is linearly interpolated to 0 at the adjacent mesh points and is identically
zero on the remainder of En+1. Equation (17) says that the test function w is
constant along the characteristics defined by (5) and that the ELLAM scheme is
structurally conservative. That is, when we define the test functions to sum to
1 on E (e.g. = 1, V(x,y) £ En+1), the test functions sum to 1
over the entire space-time domain fin+1. It is clear that the solution to (5) is not
crucial for conservation; that is, the solution of (17) for w = 1 is independent of
ii(s, t). However, the vanishing of the adjoint relation reduces to the solution of
(5)-
Our solution to (5) for the adjoint relation is called the tracking algorithm.
This tracking algorithm is crucial to the entire ELLAM scheme, because it defines


the domain and orientation of integration, Q, it defines w in the fi?*1 space-
time domain, and it defines the backtracked image x* and hence all data integrals.
What we will find is that our tracking algorithm is best described by a Heun
predictor-corrector scheme for (5). In this setting, the space-time elements and
the test functions are defined by a predictor that causes an approximate adjoint
relation to vanish. However, a better approximation, the corrector, is used to
track the solution ( trial functions ), so that our method remains faithful to the
exact adjoint relation and to the spirit of the MOC.
The guidelines for our tracking algorithm are that it be 1) consistent, 2)
simple, and 3) accurate. To see the reason for condition 1), we convert to the
integral equation representation of (5),
rtn+1
?(s) = ?+1(.s) / -Ui((s,r),r)dr for z = l,...,m. (18)
Jtn
Given the way in which the solution appears on both sides of (18), ;(s,t) has
a dual role as a) a solution of (18) and b) defining an integration path for
Vi(z(s,r),r). Thinking in terms of a Picard iteration scheme (Method of Succes-
sive Approximations), it is clear that for general Vt((s,r),r) it is only possible
to acheive an internal consistency in the dual roles of a) and b) in a limiting
sense. However for special v,-, e.g. = constant, it is possible to directly acheive
consistency.
The reason for a simple tracking algorithm is one of practicality. Tracking will
be the numerical workhorse for evaluating the crucial data integrals in (15) and
16


(16). A complex scheme can grind a numerical implementation to a halt given
that Vi depends on x(s, r), not just x;(s,r), and (18) represents a coupled system
for i(s, r). Also, our knowledge of the velocity field is generally restricted to the
current and previous time levels. Therefore, dynamic time-stepping algorithms,
such as Runge-Kutta, which require knowledge of V{ at intermediate times may
be ill-suited, if not cost-prohibitive, for our task. The simplest tracking rule for
x* in a variable coefficient setting is based on a backward Euler approximation
of (18)
x*(s) = x-^(s) + O(At)2
where
X;^(s) = x"+1(s) v"+1(xn+1(s))A£ for i = (19)
If we use the predicted characteristic given by (19) in a predictor-corrector scheme
for (18), as suggested in (22) below, we find that we can acheive the consistency
criterion when the vt- are constant along the characteristics defined by (19). That
is, the corrected value for x\ will coincide with the predicted value.
The concept of an accurate tracking algorithm needs some clarification. For
our purposes, accuracy is defined with respect to the adjoint relation and not
to the integral equation (18). We need to construct x*(p) so that a modified
form of the original adjoint relation is identically zero; however, x*(p) can be a
coarse approximation to the true backtracked image x*. Given the discussion
of the previous two criteria for a tracking algorithm, we define u(x, t) to be the
17


velocity field defined by a pointwise extrapolation of nTl+1() as a constant along
the approximate characteristics of (19)
v(x,t) = backtrack(vn+1(x)) (20)
i.e. un+1(x) = Un+1(a;) = Then, we can solve (18) exactly with U*
replacing and hence, we can make Aw = 0, where
Aw = wt + v Vw, (21)
by using the coarse sE^hcharacteristics of (19). Although the accuracy of the
approximation of (18) is not crucial for the modified adjoint relation, Aw = 0,
accuracy is necessarily all-important for the solution un+1().
With respect to the solution, we must be able to build a higher-order approxi-
mation, x*(c\ that is invisible with respect to the discretization errors committed
in the construction of (15). That is, x*(c) must direct the solution un+1(x) to
the correct (with respect to the discretization errors) history value un(x*(c)). If
*(c) fails to characterize the correct domain for un+1(x), we expect some Courant
number limitation (Cu < 1) for the accuracy, but not the stability, of the solution
technique. To define x*(c\ we return to (18) in residual velocity form
xt*(s) = a£+1(s) f u,(x*(p)(s, t), r) dr
J
wn+l
f (v u)i(a*^(s,r),T) dr + 0(At)3
where the 0(A<)3 truncation error comes from the expansion of the integrands
about the point x*(p\ Then, we employ a 2-point trapezoidal approximation to
18


obtain
*;w = *:(0)+o(At)J
where
*?c\s) = - xAt(u ti)(i'(p)(s))
(22)
= *?+1(*)-5A£(f+I(a."+*(.)) + .?(x^'>W)). (23)
Despite the appearance of (22) and (23), x*(p) is actually an O(At)2, and not an
0(A£) approximation to x*(c\ To see this fact, consider
(u v)"(*^) = At
u?+1(x) v?(x) v?+l(x) v^(x At
At
= At [0 Anr(**(P)) + 0(At)\
= -AtAiV?(x

) + 0(A£)2,
where
~Zi} = f,+^- (24)
Therefore, we have the more enlightening form of (22) that is given by
*;t0) = *;(',>+i(A +o(a)3. (25)
Since the corrected image x*(c) is an O(Ai)3 approximation, the error in the
corrected tracking scheme is invisible with respect to the 0(A£)2 discretization
error of (15) and is of sufficiently high order for the solution Tin+1(). Therefore,
we can use the x*(p) approximation to define the test functions and an approx-
imate space-time domain, fly+1, with impunity as long as we incorporate x*(c)
into our weak formulation.
19


5 APPROXIMATE WEAK FORMULATION
Using the ideas of the previous section, we will reformulate (15) for the
approximate space-time finite element, flt+1, that is formed by backtracking
E"+1 along the approximate ^(-^-characteristics given by (19). Using (1) as an
example, we consider the approximate weak formulation t+1 w Cu = i fw.
Starting with the advection component, we convert from -4(*) and v(x,t) to
-4.(-) and v(x,t) as defined in (20) and (21) to match our conversion from fl+1
to fit^+1. As a result, we obtain
/ ,, w Au = I ., w Au + I ,tw(v-v)- Vu. (26)
Jn£+1 Jn?1
The first integral on the right hand side of (26) is handled as before (cf. (9))
J_.nA(uw)- f_rnuAw.
(27)
Since ,4(*) denotes differentiation along the approximate characteristics of flt+1
and we know that Aw = 0, we can simplify (27) (cf. (10)) as
(28)
In order to get the sc'^-corrected backtracking scheme of (22) for the solution
tin+1(s), we approximate the residual advection of the second integral of the right
hand side of (26) by a trapezoidal rule as
= At J +i w (v v)n Vun(s*^) + O(Ai)3.
(29)


Combining equations (26), (28), and (29), the approximate weak form of Au = 0
becomes
= / (u"(x*(P)) lAt(v U)n Vttn(**

)) w + 0{At)3. (30)
Now, we recognize that the integrand on the right hand side of (30) is an expan-
sion of un(*(c)) about x*(p\ since
un(x Z
= un(*(p)) hAt)2Avn W*(s*(p)) + O(At)3
Z
= un(*(p)) ^At {v v)n Vun(z,(p>) + O(At)3
z
by equations (22) and (25). As a result, we can rewrite (30) as
/ "+1(*) w = f w + 0( Af)3. (31)
Therefore, the corrected tracking scheme for the solution iin+1(s) does not intro-
duce additional O(At)2 discretization error. Also, the computational cost of (31)
is limited to the evaluation of the corrected tracking relation (23).
For the reaction term, we find that
un(s*(c)) kn(x) + 0(A<)2] [fcn(s*

) + 0(A£)2 ]A£
= un(x*p)) kn(x

) At + O(At)3.
Therefore,
J un(x*

) kn(x

) At = [ n+i un(x) Jfen(s*(c>) At + 0(At)3 (32)
21


and the corrected tracking for the reaction does not introduce more 0(A£)2 error.
Also, equation (32) is acheived with no additional computational effort.
The remaining diffusion and source terms are treated as before (cf. (13), (14))
with ft"+1 replacing fi?*1. After combining (31) and (32) with the modified
versions of (13) and (14), we have the modified form of (15)
/ un+1{x) w + At / +i Dn+1Vun+1(x) Vw
= L, "(*,(c)) f1 *(*'(c))Af) +Ai L, /+I <33)
At this point, we have completed the basic presentation for computational
patches, E-+1, whose backtracked images remain in the interior of E". We will
complete our discussion for the characterization of the solution, n+1(x), by in-
vestigating computational patches, Si"+1, that backtrack to the inlet boundary,
Fn+1 C 3E x Jn+1.
22


6 INLET BOUNDARY DISCRETIZATION
The MOC tells us that points on the computational grid E?*1 backtrack to
the interior of En or to the inlet boundary, rn+1 = dS x Jn+1 where Jn+1 =
[tn, tn+1]. In the previous sections, we have seen the ELLAM formulation for the
former case. In this section, we will investigate the latter case and demonstrate
how the ELLAM formulation systematically treats the influence of inlet flux
and Dirichlet boundary conditions (3) on the solution. Indeed, the treatment
of boundary conditions represents the most dramatic departure of ELLAM from
ELM. The ELM approach uses ad hoc and artificial extensions of the boundary
so that the boundary finite elements can be treated as interior elements. The
ELLAM approach acknowledges the distinction between interior and boundary
in order to incoporate boundary conditions in a manner which is consistent with
the underlying physics of a given problem. For the purposes of notational clarity,
we will consider patches of the computational grid, Ef"+1, that backtrack via
equation (19) to points (s},£*) on the inlet boundary. In the language of the
previous section, the tracked point x} is actually x*t^ and the x\ that will appear
in the integrals is actually x*^C\ Also, we remind the reader that we are working
on the (north, south, east, and west) lateral faces of a cube formed by 3E x Jn+1.
Given that the backtracked image is (x},t*) instead of (s*,tn), we will start
by making the natural adjustments to (33) For the advection term, we find that


in (31) that /E~+i un(x*^)w is replaced by
Ln+iu(xit{c)'t*)w- (34)
Likewise, the reaction term is adjusted by replacing (32) with
f ,n+l At*(x)k(x*I{C\t*)u(x*I{C\t*)w (35)
*'E."
where
At*(x) = tn+1 t*(x).
Also, the source term becomes
/ A4-(*)/"+1M- (36)
Neglecting the diffusion term, we have the weak formulation
J +iun+1(x)w = J +iu(x*!{C),f)(l-k(x](C\f) At*(x))w
+ /' AfWr+'v (37)
JEy
which says that the solution data u(x^C\t*) propagates from the boundary in
a manner that is governed by the physical law (1) over the reduced time step
A£*(). However, we must note that (37) does not contain the influence of the
boundary conditions. For example, if the inlet has a flux boundary condition,
(im DVu) n(x*^C\t*) is known but u(x*^C\t*), which appears on the right
hand side of (37), is unknown. Therefore, our incorporation of the boundary
conditions must come from the treatment of the diffusion term.
24


For the caise of a flux boundary condition, we treat the diffusion term as in
(12), (13). However, the element boundary, is defined as the union of
the surface formed by backtracking the boundary of the computational patch,
5S;"+1, and the surface ri"+1, which is formed by by backtracking Ei+1 to the
rn+1 boundary. Since w = Wij vanishes on 5St-+1, we can write the element
boundary integral as
/ w DVu -n = f, w DVu n. (38)
Jsn>+1
Equation (38) represents one term of the flux, (yu DVu) n, on Tn+1. In order
to form the remaining flux term, we must perform a measure change for (34) from
E;"+1 to the domain of definition for the integrand, r|+1. Using the relation (19)
between x and t in we can modify (34) as
bn+iU(x*i(C)t*)w = ~ (39)
1 j
where dEi"+1 = |U*7i| dt dS and 5 is a spatial boundary variable ( see Section 8 for
the details of measure adjustment ). Since v points inward at an inlet boundary
and n is the outward normal, v n < 0 and v n = |iJ n\. For future reference,
we must note that equation (39) represents a simple, but significant, break from
the explicit tracking formulation. Indeed, we will see in a later section how
the concept of measure adjustment links the current formulation to the implicit
tracking formulation of [7].
Before incorporating (38) and (39), we must apply corrected tracking to the
integrands and replace v with v, since the flux boundary condition is phrased is
25


terms ofv and /von rn+1. Using the arguments developed in Sections 5 and
6 for the inclusion of corrected tracking and the treatment of residual velocity
terms, we can make the aforementioned modifications to (38) and (39) at the
(necessary) expense of an O(At)2 error term.
After combining the modified forms of (38) and (39) into the framework of
(37), we can bring the modified form of (38) to the right hand side and replace
the flux integrand with h(,£), by the definition (3), to form the modified version
of (33),
L *"(*) + d+i v
= ~L'n+i w /,+1 ku{x]{c\t*) At\x)w
+ 1' Af(*)r+'w. (40)
We see that in the limiting case of D = 0, h(x,t) = vu to and we obtain the
desired form (37), after making the inverse measure change from TiJl+1 to E;"+1.
In the Dirichlet case, we handle the diffusion term by inverting the order of
operations used in (12) and (13). That is, we start by performing a backward Eu-
ler time discretization and end by invoking Greens and the Divergence Theorem.
In this way, we obtain
wT>u =

+ /, > Dn+1Vun+1(x) (V(At*(x)); + A£*(*) Vw).
26


Since w = 0 on 3E-l+1, the 5Et-jl+1 boundary integral vanishes and
w Dn+lVua+l(x) V(Af(x))
+ f,nl At*{x) Dn+1 Vun+\x) Vw. (41)
The first integral on the right hand side of (41) takes the place of the T,-"*"1
boundary integral of the flux case (40). Now, we modify (34) and (35) to fit the
Dirichlet context by replacing u(x*^C\t*) with h(x^C\t*) from (3). As a result,
we have formed the Dirichlet version of (33)
L + L. n+I W+1(x) V(Af(x))
+ ^,n+i At*(x) Dn+1Vun+1(x) Vw
= / m*;( "ij
(42)
For equation (42), we obtain (37) in the limiting case, D = 0. As a basis for
comparison, implicit ELLAM yields inlet boundary equations that have the same
form as (40) and (42).
Equations (40) and (42) conclude our discussion of inlet discretization and
the specification of ttn+1(x) on E"+1. However, En+1 backtracks to a subset, E*,
of En and we have not used all of the information at our disposal. Physically, the
information residing on En E* flows to the outlet boundary, 0n+1 C 5E x Jn+1.
In the next section, we will characterize the solution on 0n+1 by discretizing 0n+1
in both time and space and then backtracking to En E*.
27


7 OUTLET BOUNDARY DISCRETIZATION
In this section, we will see how to use the wealth of information on En E*
that has not been used by un+1(s). Since this information tracks to a known
outlet boundary condition, we will solve for information on the outlet boundary
that is dual to the known boundary condition. That is, if the outlet boundary
condition is flux ( (vu DVu) n ), we will solve for the value of the solution, u,
on the outlet boundary. Also, if the outlet boundary condition is Dirichlet, we
will solve for the value of the flux, Vu ra, on the outlet boundary.
To start our outlet formulation, we assume that the outlet boundary 0n+1 C
dYt x Jn+1 is known. Alternately, we could define 0n+1 by forward tracking
- which is discussed in the next section from En E*. We will define the
time discretization of 0n+1 by 0"fc+1 with k denoting time and replacing j, which
denotes the loss of spatial dimension on 0n+1. In order to maintain the equation
framework that we have established, we can formally think of integration with
respect to a computational patch, Ej^"1, outside of En+1 that backtracks to 0^x
on the outlet boundary. Also, we can think of test functions, to;*., that are defined
on E"fc+1 such that they have a standard representation with respect to 0^x. As
was seen in the inlet boundary formulation, we will perform measure changes in
order to convert from /E+i to /0?+i. Therefore, the precise nature of E^fe+1 and
u>2b+1 do not need to be specified.
Given that the backtracked images from E^1 to Sfe1"1 and from EJ*-1 to
28


E71 E* are, respectively, (xo,t*) and (x*(c\tn), we will start by making tbe
natural adjustments to (33) For the advection term, we find that
f A(uw) = f (u(x*0,t*) un(x - J^+i(vu n(x0,t) vun -n(x^C)))w,
(43)
since v n is positive on the outlet boundary. Also, the XqC^ in /0n+1 denotes
tracking from standard points on 0^+1 and not on E"fc+1. Likewise, the reaction
term is adjusted as
/E?+1 fc" (f <) = jf ^ V n *" u"(xi) (t-e)w. (44)
Also, the source term is converted as
X-+1 /(*0n (t* -tn)w = JQn+i V n f(xo ,t){t- tn)
Finally, the diffusion is treated as in (38) and find that
w.
(45)
f ..wDVu-n+f DVu Vw (t* tn)
Jan?*1 v 1
= Jn+iwDVu-n Jn+iv-nDVu-Vw(t tlx). (46)
Now, we combine terms ((43), (44), (45), and (46)) to obtain the outlet version
of (33) given by
J^+i(vu D^7u) nw + /n+1 ^ tiDVu Viw(t tn)
fc ik
= fQn+1 V n (un(xic)) (1 kn(x$C)) (* **) )
+ ^B+1 v-nf{xQ,t)(t- tn)w.
(47)
29


For a flux outlet boundary condition, we modify (47) by noting that (cf. (3))
DVu) nw =
(48)
since v = v on 0^+1. Then, we solve for u in terms of the Vu of the second
integral on the left hand side of (47). For a Dirichlet outlet boundary condition,
we modify (47) by noting that (cf. (3))
I vu-n= I v h(xn,t) *7i (49)
/* y?+ v ' v
tk ifc
and solve for Vu n after decomposing the Vu, in the second integral on the left
hand side of (47), into tangential, Vh(x,t), and normal, Vu n, components. As
a basis for comparsion, implicit ELLAM yields outlet boundary equations that
have the same form as (48) and (49).
With equations (33), (40), (42), (47), (48), and (49), we have completed the
explicit tracking ELLAM formulation with data integrals defined with respect to
the measure at the head of the characteristic on E"+1, E"*+1, and Q^T*"1. However,
equation (39) represents an important break from this framework. In order to
incorporate a flux inlet boundary condition, we were able to perform a simple
measure change in (39) from E-"+1 to the domain of definition of the integrand,
ri"+1. In the next section, we will see how we can establish an equivalence ( up
to 0(At)2 discretization error ) between the current formulation and the implicit
tracking formulation of [7], where all data integration is performed with respect
to the measure on the domain of definition.
30


8 MEASURE ADJUSTMENT AND
IMPLICIT ELLAM
So far, our formulation has successfully incorporated many desirable prop-
erties: bridging the gap between exact and weak formulations, simple tracking for
the vanishing adjoint, and systematic treatment of general boundary conditions.
The final step is to demonstrate that this explicit ELLAM formulation conserves
mass, in order to distinguish it from its ELM predecessors. Since the test func-
tions are defined so that they sum to unity on numerical conservation of
mass corresponds to the exact evaluation of the /E+i data integral. Con-
servation can be achieved when quadrature points selected on E-+1 backtrack
to standard quadrature points for the representation of un(x) on En, since the
integration is performed with respect to E?*1 and the integrand is defined on
En. Given the nature of our tracking relation (19), it appears that conservation
would be difficult to achieve with backtracking. In general, our only alternative
would be to saturate En with backtracked quadrature points, until a desired error
tolerance is achieved. The best case scenario would be to change the measure
of the data integral from E-+1 to the measure of the local domain of definition,
E,without additional truncation error.
Given the simple nature of (19), we can define local measure adjustment from
to Ej in our approximate weak formulation. First, we consider the interior


case. Since vn+1(x) = we can use the alternate form of (19)
*(P)(s) = ?+1(s)
to find that
<*r* ., .,
dx:F>
Therefore,
(P)
and
HSp*) **
ir =I (*) n (i+^rA()
(50)
(51)
(52)
(53)
represents the conversion to Eulerian coordinates for the interior data integral.
Furthermore, notice the (53) is achieved without additional O(At)2 error.
For measure adjustment on the boundary, we use (51) in a chain rule argu-
ment. That is, the boundary version of (50)
(Sfl)i = *."+1 -Vi(*B,f)(tn+I-f) (54)
tells us that
(1+^fer(+1
d(xB)j
dt*
Vi(xB,t*).
(55)
Therefore, the boundary version of (51) ( which is represented by the first term
on the right hand side of (55) ) and the chain rule tells us that
dx?+l d{xB)i
d(xB)i dt*
= Vi(xB,t)
(56)
32


represents the measure adjustment for the boundary for the spatial variables that
are replaced by time when we track to the space-time boundaries. Since v-n tells
us which velocity components are responsible for tracking to the boundary, |tJ n\
gives us the correct ( | | for area preservation ) measure adjustment for the inlet
and outlet boundary integrals of Sections 6 and 7. For the remaining spatial
variables, we use the boundary version of (51).
Before moving on to the implementation of measure adjustment, we will es-
tablish the equivalence between the current, explicit ELLAM formulation and
the implicit ELLAM formulation of [7]. The implicit formulation assumes that
some form of backtracking has been performed in order to define space-time finite
elements with respect to Eulerian, and not Lagrangian, coordinates. Then, the
weak form of (1) and (2) is posed with respect to an Eulerian grid and we perform
standard integration by parts to form an adjoint relation. Since the vanishing
adjoint induces Lagrangian concepts, we can apply the tracking and approximate
weak formulation ideas from the explicit ELLAM formulation to the implicit
ELLAM formulation. After making the appropriate adjustments, the implicit
ELLAM formulation ( ignoring boundary terms ) treats the advection operator
as
/ f (ut + V''Vu)w= f un+1(x)w f un(x)
Jjn+lJ £;/ JVij JZJ,
/ (wt + V -Vw)u I I (V-v)uw
Jj"+1 JUii Jj+1 J-Sii
w
33


Since the approximate adjoint vanishes, the weak form of ut + v Vu = 0 becomes
f un+1(x)w f f (V -v)uw = f un(x)w. (57)
7jti+1 JEjj
Because the sense of the integration in (57) is Eulerian, and not Lagrangian,
the measure of the /s. data integral has the desired form for mass conservation.
However, the implicit ELL AM formulation does not tell us how to treat the V-v
integral on the left hand side of (57). Although much of the external trappings
of the explicit ELLAM formulation are missing, the Lagrangian reference frame
that is induced by the vanishing adjoint is at the core of the implicit ELLAM
formulation. From our discussion of measure adjustment for the explicit ELLAM
formulation, we recognize the V v integral as a reference frame adjustment
from Eulerian coordinates to the Lagrangian coordinates that are induced by
backtracking. Using forward Euler discretization, we see that (57) becomes
f ttn+1() w = f u"(x) (1 + V v"(x)) w + 0(A£)2. (58)
J'Zij Jjzr.
Reverting back to the explicit ELLAM counterpart, (53), of (58), we see that
ft (l + A<) = 1 + V vn + O(At)2 (59)
and conclude that (53) and (58) are equivalent to within O(At)2 discretization
error. For the divergence form (2), explicit ELLAM with measure adjustment
treats the weak form of Ati + uS7 v = 0 as
/ n+i un+1(x)w = /. "(*) w (1 V vn(x)At) (1 + V vn()At) + O(At)2
EJ
= J -un(s) w + 0(At)2,
34


which is identical to the result achieved by implicit ELLAM. Therefore, the inte-
rior implicit and explicit ELLAM formulations are equivalent to within O(Ai)2
discretization error. Finally, we note that similar results follow for inlet and
outlet elements of the implicit and explicit ELLAM formulations.
Measure and reference frame adjustment establish a formal equivalance be-
tween implicit and explicit ELLAM. More importantly, these formal adjustments
enable both ELLAM formulations to conserve mass. However, we must be able
to implement measure and reference frame adjustment. In the next section, we
will identify the mechanism behind measure and reference frame adjustment and
outline a numerical implementation.
35


9 FORWARD TRACKING
Although (53) gives us the formal measure adjustment needed for mass con-
servation and (58) gives the formal reference frame adjustment to Lagrangian
coordinates, we need to evaluate the integrands of (53) and (58) in a numerical
implementation. Since measure and reference frame adjustment depends on a
knowledge of U, we see that the mechanism behind measure and reference frame
adjustment is the forward tracking relation (cf. (50))
?+1 = * + Vi(x*)At, (60)
where xf is known and a:?*1 is the head of the characteristic under forward
tracking. In (53), equation (60) would be used to evaluate w on its domain of
standard definition, E?*1.
As an illustration of a forward tracking scheme, consider the case of a piecewise
linear Raviart-Thomas representation of the velocity field. That is, u+1() is
piecewise linear in Xi and piecewise constant in Xj, for i ^ j. Then, (19)
becomes
= *rtl-(i*r+1+6)^ (6i)
and solving for a+1 we find that
s?+1 = + [(aiAi)x; + &iAi(l + aiAt)]-|-0(At)2
= *? + (a'.*? + b^At + OiAt)2. (62)
Therefore, (62) says the velocity field remains approximately piecewise linear


Raviart-Thomas under backtracking. Since (62) is an approximation to the for-
ward tracking relation (60), we can use (62) to determine a provisional value for
the head of the characteristic, "+1- Then, equation (61), which defines a:+1
implicitly, can be used in an iteration scheme to determine the true value of x+1
and U", since ^(e*) = v7l+1(n+1). By the Implicit Function and Fixed Point
Theorems, we can find z+1 via (61), if
du?+1(s)
dxi
At
< 1.
(63)
Since Vi has a piecewise definition, (63) can be restated as
|a;At| =
Au?+1
A Xi
At
= |A(C7u)i| < 1
(64)
and can be interpreted as a local shock condition. That is, whenever we avoid
shocks and the characteristics are well-defined, equations (63) and (64) are satis-
fied and exact forward tracking is possible. Therefore, conservation of mass and
reference frame adjustment, which depend on forward tracking, is possible.
37


10 SUMMARY
ELLAM represents a natural strategy for the treatment of advection in
advection-reaction-diffusion equations. In Sections 2 and 3, we were given equa-
tions (1) and (2) that contain advection, v Vu, terms in Eulerian coordinates.
Guided by the MOC, we rewrote our equations in Lagrangian coordinates, in
order to remove the advection terms. Numerically, the solutions of our origi-
nal, Eulerian, equations have undesirable, oscillatory properties; however, the
solutions to our transformed, Lagrangian, equations have desirable, smoothing
properties. Having seen the promise of ELLAM, we are faced with the choice
of explicit and implicit ELLAM formulations. In Sections 4 and 5, we saw that
the explicit ELLAM formulation is best suited for incorporating the inherently
Lagrangian concepts of vanishing adjoints, tracking algorithms, and approximate
weak formulations. On the other hand, the implicit ELLAM formulation is best
suited for incorporating the inherently Eulerian concepts of mass conservation
and boundary conditions. However, in Sections 6, 7 and 8, we saw that these
l
crucial concepts can be incorporated into the explicit ELLAM formulation. Addi-
tionally, the explicit ELLAM concept of measure adjustment allowed us to bridge
the gap between implicit and explicit ELLAM and to freely exchange ideas be-
tween the two formulations. In Section 9, we saw that forward tracking was a
crucial component of and a crucial link between both formulations. Armed with
a knowledge of implicit and explicit formulations, the potential ELLAM user has


the freedom to create hybrid schemes that are tailored to specific problems. In-
deed, consistency and flexibility axe the hallmark of implicit and explicit ELLAM
and represent ELLAMs greatest strength.
39


References
[1] A. M. Baptista, Solution of advection-dominated transport by Eulerian-
Lagrangian methods using the backwards method of characteristics, Ph. D.
Thesis, Massachusetts Institute of Technology, 1987.
[2] M. A. Celia, I. Herrera, and E. T. Bouloutas, Adjoint Petrov-Galerkin
methods for multi-dimensional flow problems, Proceedings of the Seventh
International Conference on Finite Element Methods in Flow Problems,
Huntsville, AL, April 3-7, 1989, pp. 953-958.
[3] M. A. Celia, T. F. Russell, I. Herrera, and R. E. Ewing, An Eulerian-
Lagrangian localized adjoint method for the advection-diffusion equation,
Adv. Water Resour., to appear.
[4] J. Douglas, Jr., and T. F. Russell, Numerical methods for convection-
dominated diffusion problems based on combining the method of character-
istics with finite element or finite difference procedures, SIAM J. Numer.
Anal, 19 (1982), 871-885.
[5] T. J. R. Hughes, ed., Finite Element Methods for Convection Dominated
Flows, AMD Vol. 34, ASME, 1979.
[6] T. F. Russell, Eulerian-Lagrangian localized adjoint methods for advection-
dominated problems, Proceedings of the 13th Biennial Conference on Nu-
merical Analysis, University of Dundee, Pitman, 1990, to appear.
[7] T. F. Russell and R. V. Trujillo, Eulerian-Lagrangian localized adjoint
methods with variable coefficients in multiple dimensions, Proceedings of
the 8th International Conference on Computational Methods in Water Re-
sources, Venice, Italy, June 11-15, 1990, pp. 357-363.
40