- Email: [email protected]
Finite element or finite difference methods for the two-dimensional shallow water equations?
COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 7 (1976) 351-357 0 NORTH-HOLLAND PUBLISHING COMPANY
FINITE ELEMENT OR FINITE DIFFERENCE METHODS FOR THE TWO-DIMENSIONAL SHALLOW WATER EQUATIONS? T.J. WEARE Hydraulics Research Station, Oxon OX1 0 8BA, Wallingford, UK
Received 21 April 1975 Revised manuscript received 30 June 1975
The computational costs of some finite element and finite difference methods currently used to solve the shallow water equations are compared on theoretical grounds. It is shown that because band algorithms are employed, the finite element methods considered are not economically attractive for practical calculations. This weakness is not intrinsic to the finite element approach and ways in which it might be avoided are briefly discussed.
1. Introduction In recent years considerable effort [ l-31 has been devoted to the problem of applying the method of finite elements [4] to the solution of the two-dimensional shallow water equations. Hitherto most workers have been concerned with those matters of a fundamental nature which arise when adapting a well established method to a completely new class of problems. Less attention has been paid to the question of the relative economics of the finite element method compared with the widely used method of finite differences. The following theoretical considerations indicate that the finite element methods currently proposed in the literature [ l-31 are uneconomic for practical calculations because band algorithms are employed in the solution. The usual starting point for such models is the two-dimensional shallow water equations which are applicable to nearly horizontal flow under well mixed conditions. These may be written in terms of the surface elevation f and the x, y components of the depth-averaged velocity U, u as follows:
&.$pd)-$(ud)-f,, (1)
where i - at/at etc., f is a friction factor, ter. The problem is to solve these equations
d is the total water depth and a is the Coriolis paramewhen the initial conditions
of the variables {, u and u
352
T.J. Weare,Finite element or finite difference methods
and appropriate boundary conditions are given. To facilitate consideration of the relative economics of finite difference and finite element methods, a brief description of the salient features of each is given. In both methods the continuous variables {, U, u etc. are replaced by a set of values on a discrete space-time mesh.
2. The finite element method In the finite element method each variable is expanded in terms of a discrete set of basis functions Nj(x, y) and nodal values ci(t), z+(t) and vi(t). The system of coupled equations which results is assembled into a single matrix equation of the form Mir-F(U)=o,
(2)
where U is the vector of nodal values n = total number of nodes . The vector F is a nonlinear function of the nodal values. The 3n X 3n global matrix M is a function of the element geometry only. in general M is symmetric and banded with a semiband-width 3m, where m is the largest difference in node numbers that occurs in any element. The central task in the finite element method is thus that of integrating in time a set of coupled equations (2) for the nodal values from given initial conditions and subject to specified boundary constraints. The question of accommodating the boundary constraints is not relevant to the present discussion and will not be considered here, Many numerical integration schemes have been developed to handle this task. They can be divided into two basically different approaches - iterative and direct methods. Since the equations to be integrated are nonlinear, care must be exercised in order to ensure that the integration is stable. Iterative integration
schemes
The models described in [ 1] and [ 21 handle the nonlinearity of equation (2) by an iterative procedure. For example, in the trapezoidal scheme the first approximation for the unknown values at time t + At, denoted UT, is obtained from the known values of U at time t through the approximation
and successive approximations Mz
1
(u;l
- u) - ;(F(U;)
are obtained by + F(U)) = o .
(4)
Convergence of this iterative procedure is usually accelerated by a relaxation technique. Both approximations (3) and (4) can be written in the form
T.J. Weare,Finite element or finite difference methods
353
M(Ui+I - u> = Fi 3
(5)
where Fi = i At(F(Ui’> + F(U)) so that the vectors U and Fi are known. The standard method of solving the matrix equation (5) is by Gaussian elimination. Since the matrix to be inverted in order to solve the coupled equations in (5) is constant, the inversion has to be carried out once only. However, at each iteration the vector Fi has to be used in a forward substitution, and a back substitution must be performed (that is the matrix multiplication M-‘Fi has to be repeated). If i iterations are necessary per time step, the number of operations involved is i X 2 X (3~) X (3m)/n
= 18 mi /node/time
step .
(6)
The essential point to note is that the cost of generating to the band-width, i.e. to m.
each new nodal value is proportional
Direct integration schemes In the second approach [ 31 the nonlinear term F(U) is written and equation (2) is represented by the Crank-Nicolsin form M it
(U+-
U) - $U);(U++
in a quasilinear
form E(u)U,
U) = o ,
i.e. M+U+=M-U,
(7)
where M’ =M T ;At&J)
.
Since the matrix M+ is time dependent, the matrix inversion involved in solving equation (7) has to be repeated at each time step. If we assume that M+ is a 3n X 3n symmetric banded matrix with semiband-width 3m, then the number of operations involved in the solution of (7) is 3n($(3m)2 + 2 X 3m)/n = ym2 + 18m /node/time Thus this form of implicit finite element depends quadratically on the band-width.
method
step . leads to a computational
(8)
cost per node that
3. The finite difference method There are many finite difference methods in current use. The usual technique is to approximate the partial derivatives in equations (1) directly by differences of values on a uniform grid. A continuous variable such as {(x, y, t) is replaced by a discrete set of variables {(iAx, kAy, ZAt) = cjlk on a uniform space-time mesh with spacings of Ax, Au and At in the x, y and t directions respectively. Partial derivatives can then be replaced by second order approximations such as
3.54
T. J. Weare, Finite element or fiizite difference methods
i-;_s,k)+ @Ax2). As with the finite element method the finite difference approximation or implicit numerical integration schemes.
may lead to either explicit
Explicit finite differences rnet~~ds
Matrix techniques are not generally required in explicit finite difference methods, and a stable solution is possible without resorting to an iterative procedure [ 5,6]. The essential feature of this approach is that the right hand side of equations (1) are represented in terms of known time levels. For example the first equation might be represented by the finite difference equation
(10) Thus if the variables at time levels 1 and 1 + fr are known, equation (10) can be used to predict the elevation {k’ at time level E+l directly. In this way the computational effort required to advance a solution is simply proportional to the number of grid points involved. In order to draw a comparison with the finite element methods discussed in the previous section, it is necessary to define a unit of work. In any FORTRAN program set up to handle the triangularisation and back substitution involved in the Gaussian elimination the basic operation is of the form DO 10 I = ...
10
B (I) = B (I) - A(1, J) * B(J) .
For an ICL 1904s computer this is equivalent to approximately 7 multiplications. A basic unit of work of 5 multiplications is assumed for simplicity. Obviously this unit is machine dependent, but since the same unit is used throughout, the comparisons should be generally valid. An explicit finite difference scheme used by the author [ 51 requires approximately 40 units of work/grid point/time The computational
(11)
step .
cost per grid point is independent
of the overall size of the grid.
Implicit finite difference methods
To avoid the sometimes severely restrictive stability constraint of explicit schemes, it is often desirable to use an implicit integration scheme. This can be achieved by representing the main spatial derivative terms in equations (1) by a Crank-Nicolson form, that is by averaging between the known and unknown time levels. However, using this representation simultaneously for both the x and y derivatives leads to a two-dimensional coupling of the equations for the unknowns. That is, the finite difference equations will be similar to the matrix equation (7) arising in the implicit finite element method, and consequently the process of matrix inversion and multiplication
T. J. Weare, Finite element or finite difference
methods
355
will be required at each time step. Thus the computational cost of this scheme is the same as given in (8) where m now corresponds to the number of grid points across the width of the finite difference mesh. The usual way of circumventing this two-dimensional coupling is to split the implicit representation of the spatial derivative terms into two steps by the alternating direction method [ 5,7,8]. In the first step the appropriate x derivatives (say) are treated implicitly. The integration is completed in the second step when the y derivatives are treated implicitly. The coupling is thereby reduced to one dimension in each step and the matrices appearing in equation (7) become tridiagonal. Since the band-width is then fixed, the computational effort required with the alternating direction implicit finite difference method is simply proportional to the number of grid points, as with the explicit finite difference scheme. An implicit scheme used by the author [ 51 requires approximately 50 units of work/grid
4. Relative economics
point/time
step .
of finite elements
(12)
and finite differences
The most important feature of the economics of both types of finite element methods described in (2) is that the unit cost of computation as represented by (6) and (8) increases with the bandwidth, whereas with the finite difference methods it is independent of the overall size of the model. The relatively low computational efficiency of these finite element methods stems from the sparseness of the global matrix. That is, although most of the matrix elements within the populated band are zero, generally this fact does not lead to any computational simplification in the Gaussian elimination. The number of nodes required for a finite element solution will in general be less than the number of grid points required for a finite difference solution of the same problem. However, with the use of a patched grid finite difference scheme (that is one in which the uniform grid is reduced in scale where necessary by discrete jumps, i.e. by a factor of 2 or 3, say) the saving is unlikely to be greater than a factor of 2 for most problems. Allowing for this factor of 2 it follows from (6), (8), (1 1) and (12) that to be economically viable, the finite element methods must satisfy the following inequalities: (i) explicit methods, iterative schemes 18mi< 2X40,
i.e.
mi<4,
(13)
where i, the number of iterations/time step, is typically (ii) implicit methods, direct scheme 21 urn’
+
18~ < 2 X 50,
i.e.
m<2.
- 4 to 7 (see [ 1,2]),
(14)
These restrictions are, if anything, generous, since only the work involved in the central task of the finite element solution has been estimated. The cost of setting up the vector Fi in equation (5) or the matrices M+ and M- in equation (7) together with general ‘housekeeping’ costs (boundary conditions, for example) have been ignored for the finite element method. This being so and bearing in mind that m, the semiband-width, is the maximum number of nodes across the width of the
356
T. J. Weare, Finite element or finite difference
methods
two-dimensional region being modelled, it is clear that the above economic constraints are so severe as to prohibit the use of present finite element methods for most practical problems in tidal hydraulics. For example [ 91, the finite element method of [ 1 ] when run on a CDC 7600 machine took 174 seconds of CPU time to perform a calculation of 480 time steps for a model with n = 64 nodes and a semiband-width YIZ= 22. Allowing for the fact that the CDC 7600 has at least 20 times the speed of the ICL 1904S, the explicit finite difference scheme [ 5 I would require only 6 seconds for the same number of time steps but with twice as many grid points.
5. Other considerations It is genarlly argued [ l-33 that the finite element method is better suited than the finite difference method for the problem of modelling complex geometrical configurations and boundaries. That fewer nodes will be needed in a finite element representation of a given coastline is certainly true. However, compared with a patched grid finite difference model the saving is unlikely to amount~to more than a factor of 2. Furthermore, the difficulty with boundaries is not just one of achieving adequate spatial resolution but also a problem of modeling their time dependence (beaches and intertidal flats, for example). This problem is readily tackled by finite difference methods but is likely to be less amenable to the finite element approach. For example, treating a moving coastline by the obvious method of adjusting the shape of the boundary elements would involve the costly procedure of inverting a new global matrix M at each time step. Although generally of secondary importance, considerations of computer storage requirements can sometimes be inhibiting. The minimum storage requirement for the finite element methods described in section 2 is dictated largely by the size of the banded matrices involved, i.e. - 3n X 3m = 9mn real numbers, that is 9~2 real numbers per node. Note that the effectiveness of the storage decreases as the band-width increases. The storage requirements of finite difference methods vary from method to method, but a typical figure is - 15n real numbers, that is 15 real numbers per grid point.
6. Conclusions The considerations outlined above show that the finite element methods currently being applied to two-dimensional models in tidal hydraulics suffer a severe economic disadvantage when compared with finite difference methods. The unit cost of these finite element methods increases with the band-width of the computational mesh, whereas for the finite difference methods it is independent of the size of the model. This weakness is not intrinsic to the finite element method. It arises from the use of band algorithms to solve the sparse matrix equations which result in this method of discretisation. Future development in this field should clearly concentrate on applying more efficient matrix handling techniques. Perhaps the well known Gauss-Seidel method (see 141) could be used. Alternatively, by using elements which conform to a distorted uniform grid of nodes with a constant width, it may be possible to exploit the sparseness by adopting an alternating direction solution in much the same way as in the implicit finite difference method. Another attractive possibility is the lumped mass system [ 101. It is essential that the disadvantage of a band-width dependent efficiency should be overcome
T. J. Weare, Finite element or finite difference
before the full potential cations.
of the finite element
method
methods
can be realised in hydraulic
351
engineering
appli-
Acknowledgement The author is indebted to Drs. G.V. Miles, P. Hood and J.M. Davis for helpful discussions. The work described in this paper has been carried out as part of the research program of the Hydraulics Station, and is published with the permission of the Director of Hydraulics Research.
References [ 1] C. Taylor and J. Davis, Tidal and long wave propagation - a finite element approach, Civil Eng. Dept., Univ. Wales, Swansea C/RI 89112. [ 21 J. Connor and J. Wang, Finite element modelling of hydrodynamic circulation, International Conference on Numerical Methods in Fluid Dynamics, Univ. Southampton, Sept. 1973. [3] G. Grotkop, Finite element analysis of long-period water waves, Comp. Meth. Appl. Mech. Eng. 2 (1973) 147-157. [4] O.C. Zienkiewicz, The finite element method in engineering science (McGraw-Hill, London, 1971). [5] T. J. Weare, Report on the results from the numerical models of the Outer Thames estuary, Hydraulics Research Station report EX 656, Jul. 1974. [6] D. Prandle, A numerical model of the Southern North Sea and River Thames, Inst. Oceanographic Sciences, Bidston Observatory, Report No. 4, 1974. [7] A.R. Gourlay and J.L. Morriss, Finite-difference methods for nonlinear hyperbolic systems, Math. Comp. 22 (1968) 28-39, 549-556. [S] J. J. Leendertse, A waterquality simulation model for well-mixed estuaries and coastal seas: Vol. 1, Principles of computation, The Rand Corporation, RM-6230-RC, Feb. 1970. [ 91 J.M. Davis, Water Research Centre, Medmenham, England. Private communication. [lo] S.L. Smith and CA. Brebbia, Finite element solution of Navier-Stokes equations for transient two-dimensional incompressible flow, J. Comp. Phys. 17 (1975) 2355245.
FINITE ELEMENT OR FINITE DIFFERENCE METHODS FOR THE TWO-DIMENSIONAL SHALLOW WATER EQUATIONS? T.J. WEARE Hydraulics Research Station, Oxon OX1 0 8BA, Wallingford, UK
Received 21 April 1975 Revised manuscript received 30 June 1975
The computational costs of some finite element and finite difference methods currently used to solve the shallow water equations are compared on theoretical grounds. It is shown that because band algorithms are employed, the finite element methods considered are not economically attractive for practical calculations. This weakness is not intrinsic to the finite element approach and ways in which it might be avoided are briefly discussed.
1. Introduction In recent years considerable effort [ l-31 has been devoted to the problem of applying the method of finite elements [4] to the solution of the two-dimensional shallow water equations. Hitherto most workers have been concerned with those matters of a fundamental nature which arise when adapting a well established method to a completely new class of problems. Less attention has been paid to the question of the relative economics of the finite element method compared with the widely used method of finite differences. The following theoretical considerations indicate that the finite element methods currently proposed in the literature [ l-31 are uneconomic for practical calculations because band algorithms are employed in the solution. The usual starting point for such models is the two-dimensional shallow water equations which are applicable to nearly horizontal flow under well mixed conditions. These may be written in terms of the surface elevation f and the x, y components of the depth-averaged velocity U, u as follows:
&.$pd)-$(ud)-f,, (1)
where i - at/at etc., f is a friction factor, ter. The problem is to solve these equations
d is the total water depth and a is the Coriolis paramewhen the initial conditions
of the variables {, u and u
352
T.J. Weare,Finite element or finite difference methods
and appropriate boundary conditions are given. To facilitate consideration of the relative economics of finite difference and finite element methods, a brief description of the salient features of each is given. In both methods the continuous variables {, U, u etc. are replaced by a set of values on a discrete space-time mesh.
2. The finite element method In the finite element method each variable is expanded in terms of a discrete set of basis functions Nj(x, y) and nodal values ci(t), z+(t) and vi(t). The system of coupled equations which results is assembled into a single matrix equation of the form Mir-F(U)=o,
(2)
where U is the vector of nodal values n = total number of nodes . The vector F is a nonlinear function of the nodal values. The 3n X 3n global matrix M is a function of the element geometry only. in general M is symmetric and banded with a semiband-width 3m, where m is the largest difference in node numbers that occurs in any element. The central task in the finite element method is thus that of integrating in time a set of coupled equations (2) for the nodal values from given initial conditions and subject to specified boundary constraints. The question of accommodating the boundary constraints is not relevant to the present discussion and will not be considered here, Many numerical integration schemes have been developed to handle this task. They can be divided into two basically different approaches - iterative and direct methods. Since the equations to be integrated are nonlinear, care must be exercised in order to ensure that the integration is stable. Iterative integration
schemes
The models described in [ 1] and [ 21 handle the nonlinearity of equation (2) by an iterative procedure. For example, in the trapezoidal scheme the first approximation for the unknown values at time t + At, denoted UT, is obtained from the known values of U at time t through the approximation
and successive approximations Mz
1
(u;l
- u) - ;(F(U;)
are obtained by + F(U)) = o .
(4)
Convergence of this iterative procedure is usually accelerated by a relaxation technique. Both approximations (3) and (4) can be written in the form
T.J. Weare,Finite element or finite difference methods
353
M(Ui+I - u> = Fi 3
(5)
where Fi = i At(F(Ui’> + F(U)) so that the vectors U and Fi are known. The standard method of solving the matrix equation (5) is by Gaussian elimination. Since the matrix to be inverted in order to solve the coupled equations in (5) is constant, the inversion has to be carried out once only. However, at each iteration the vector Fi has to be used in a forward substitution, and a back substitution must be performed (that is the matrix multiplication M-‘Fi has to be repeated). If i iterations are necessary per time step, the number of operations involved is i X 2 X (3~) X (3m)/n
= 18 mi /node/time
step .
(6)
The essential point to note is that the cost of generating to the band-width, i.e. to m.
each new nodal value is proportional
Direct integration schemes In the second approach [ 31 the nonlinear term F(U) is written and equation (2) is represented by the Crank-Nicolsin form M it
(U+-
U) - $U);(U++
in a quasilinear
form E(u)U,
U) = o ,
i.e. M+U+=M-U,
(7)
where M’ =M T ;At&J)
.
Since the matrix M+ is time dependent, the matrix inversion involved in solving equation (7) has to be repeated at each time step. If we assume that M+ is a 3n X 3n symmetric banded matrix with semiband-width 3m, then the number of operations involved in the solution of (7) is 3n($(3m)2 + 2 X 3m)/n = ym2 + 18m /node/time Thus this form of implicit finite element depends quadratically on the band-width.
method
step . leads to a computational
(8)
cost per node that
3. The finite difference method There are many finite difference methods in current use. The usual technique is to approximate the partial derivatives in equations (1) directly by differences of values on a uniform grid. A continuous variable such as {(x, y, t) is replaced by a discrete set of variables {(iAx, kAy, ZAt) = cjlk on a uniform space-time mesh with spacings of Ax, Au and At in the x, y and t directions respectively. Partial derivatives can then be replaced by second order approximations such as
3.54
T. J. Weare, Finite element or fiizite difference methods
i-;_s,k)+ @Ax2). As with the finite element method the finite difference approximation or implicit numerical integration schemes.
may lead to either explicit
Explicit finite differences rnet~~ds
Matrix techniques are not generally required in explicit finite difference methods, and a stable solution is possible without resorting to an iterative procedure [ 5,6]. The essential feature of this approach is that the right hand side of equations (1) are represented in terms of known time levels. For example the first equation might be represented by the finite difference equation
(10) Thus if the variables at time levels 1 and 1 + fr are known, equation (10) can be used to predict the elevation {k’ at time level E+l directly. In this way the computational effort required to advance a solution is simply proportional to the number of grid points involved. In order to draw a comparison with the finite element methods discussed in the previous section, it is necessary to define a unit of work. In any FORTRAN program set up to handle the triangularisation and back substitution involved in the Gaussian elimination the basic operation is of the form DO 10 I = ...
10
B (I) = B (I) - A(1, J) * B(J) .
For an ICL 1904s computer this is equivalent to approximately 7 multiplications. A basic unit of work of 5 multiplications is assumed for simplicity. Obviously this unit is machine dependent, but since the same unit is used throughout, the comparisons should be generally valid. An explicit finite difference scheme used by the author [ 51 requires approximately 40 units of work/grid point/time The computational
(11)
step .
cost per grid point is independent
of the overall size of the grid.
Implicit finite difference methods
To avoid the sometimes severely restrictive stability constraint of explicit schemes, it is often desirable to use an implicit integration scheme. This can be achieved by representing the main spatial derivative terms in equations (1) by a Crank-Nicolson form, that is by averaging between the known and unknown time levels. However, using this representation simultaneously for both the x and y derivatives leads to a two-dimensional coupling of the equations for the unknowns. That is, the finite difference equations will be similar to the matrix equation (7) arising in the implicit finite element method, and consequently the process of matrix inversion and multiplication
T. J. Weare, Finite element or finite difference
methods
355
will be required at each time step. Thus the computational cost of this scheme is the same as given in (8) where m now corresponds to the number of grid points across the width of the finite difference mesh. The usual way of circumventing this two-dimensional coupling is to split the implicit representation of the spatial derivative terms into two steps by the alternating direction method [ 5,7,8]. In the first step the appropriate x derivatives (say) are treated implicitly. The integration is completed in the second step when the y derivatives are treated implicitly. The coupling is thereby reduced to one dimension in each step and the matrices appearing in equation (7) become tridiagonal. Since the band-width is then fixed, the computational effort required with the alternating direction implicit finite difference method is simply proportional to the number of grid points, as with the explicit finite difference scheme. An implicit scheme used by the author [ 51 requires approximately 50 units of work/grid
4. Relative economics
point/time
step .
of finite elements
(12)
and finite differences
The most important feature of the economics of both types of finite element methods described in (2) is that the unit cost of computation as represented by (6) and (8) increases with the bandwidth, whereas with the finite difference methods it is independent of the overall size of the model. The relatively low computational efficiency of these finite element methods stems from the sparseness of the global matrix. That is, although most of the matrix elements within the populated band are zero, generally this fact does not lead to any computational simplification in the Gaussian elimination. The number of nodes required for a finite element solution will in general be less than the number of grid points required for a finite difference solution of the same problem. However, with the use of a patched grid finite difference scheme (that is one in which the uniform grid is reduced in scale where necessary by discrete jumps, i.e. by a factor of 2 or 3, say) the saving is unlikely to be greater than a factor of 2 for most problems. Allowing for this factor of 2 it follows from (6), (8), (1 1) and (12) that to be economically viable, the finite element methods must satisfy the following inequalities: (i) explicit methods, iterative schemes 18mi< 2X40,
i.e.
mi<4,
(13)
where i, the number of iterations/time step, is typically (ii) implicit methods, direct scheme 21 urn’
+
18~ < 2 X 50,
i.e.
m<2.
- 4 to 7 (see [ 1,2]),
(14)
These restrictions are, if anything, generous, since only the work involved in the central task of the finite element solution has been estimated. The cost of setting up the vector Fi in equation (5) or the matrices M+ and M- in equation (7) together with general ‘housekeeping’ costs (boundary conditions, for example) have been ignored for the finite element method. This being so and bearing in mind that m, the semiband-width, is the maximum number of nodes across the width of the
356
T. J. Weare, Finite element or finite difference
methods
two-dimensional region being modelled, it is clear that the above economic constraints are so severe as to prohibit the use of present finite element methods for most practical problems in tidal hydraulics. For example [ 91, the finite element method of [ 1 ] when run on a CDC 7600 machine took 174 seconds of CPU time to perform a calculation of 480 time steps for a model with n = 64 nodes and a semiband-width YIZ= 22. Allowing for the fact that the CDC 7600 has at least 20 times the speed of the ICL 1904S, the explicit finite difference scheme [ 5 I would require only 6 seconds for the same number of time steps but with twice as many grid points.
5. Other considerations It is genarlly argued [ l-33 that the finite element method is better suited than the finite difference method for the problem of modelling complex geometrical configurations and boundaries. That fewer nodes will be needed in a finite element representation of a given coastline is certainly true. However, compared with a patched grid finite difference model the saving is unlikely to amount~to more than a factor of 2. Furthermore, the difficulty with boundaries is not just one of achieving adequate spatial resolution but also a problem of modeling their time dependence (beaches and intertidal flats, for example). This problem is readily tackled by finite difference methods but is likely to be less amenable to the finite element approach. For example, treating a moving coastline by the obvious method of adjusting the shape of the boundary elements would involve the costly procedure of inverting a new global matrix M at each time step. Although generally of secondary importance, considerations of computer storage requirements can sometimes be inhibiting. The minimum storage requirement for the finite element methods described in section 2 is dictated largely by the size of the banded matrices involved, i.e. - 3n X 3m = 9mn real numbers, that is 9~2 real numbers per node. Note that the effectiveness of the storage decreases as the band-width increases. The storage requirements of finite difference methods vary from method to method, but a typical figure is - 15n real numbers, that is 15 real numbers per grid point.
6. Conclusions The considerations outlined above show that the finite element methods currently being applied to two-dimensional models in tidal hydraulics suffer a severe economic disadvantage when compared with finite difference methods. The unit cost of these finite element methods increases with the band-width of the computational mesh, whereas for the finite difference methods it is independent of the size of the model. This weakness is not intrinsic to the finite element method. It arises from the use of band algorithms to solve the sparse matrix equations which result in this method of discretisation. Future development in this field should clearly concentrate on applying more efficient matrix handling techniques. Perhaps the well known Gauss-Seidel method (see 141) could be used. Alternatively, by using elements which conform to a distorted uniform grid of nodes with a constant width, it may be possible to exploit the sparseness by adopting an alternating direction solution in much the same way as in the implicit finite difference method. Another attractive possibility is the lumped mass system [ 101. It is essential that the disadvantage of a band-width dependent efficiency should be overcome
T. J. Weare, Finite element or finite difference
before the full potential cations.
of the finite element
method
methods
can be realised in hydraulic
351
engineering
appli-
Acknowledgement The author is indebted to Drs. G.V. Miles, P. Hood and J.M. Davis for helpful discussions. The work described in this paper has been carried out as part of the research program of the Hydraulics Station, and is published with the permission of the Director of Hydraulics Research.
References [ 1] C. Taylor and J. Davis, Tidal and long wave propagation - a finite element approach, Civil Eng. Dept., Univ. Wales, Swansea C/RI 89112. [ 21 J. Connor and J. Wang, Finite element modelling of hydrodynamic circulation, International Conference on Numerical Methods in Fluid Dynamics, Univ. Southampton, Sept. 1973. [3] G. Grotkop, Finite element analysis of long-period water waves, Comp. Meth. Appl. Mech. Eng. 2 (1973) 147-157. [4] O.C. Zienkiewicz, The finite element method in engineering science (McGraw-Hill, London, 1971). [5] T. J. Weare, Report on the results from the numerical models of the Outer Thames estuary, Hydraulics Research Station report EX 656, Jul. 1974. [6] D. Prandle, A numerical model of the Southern North Sea and River Thames, Inst. Oceanographic Sciences, Bidston Observatory, Report No. 4, 1974. [7] A.R. Gourlay and J.L. Morriss, Finite-difference methods for nonlinear hyperbolic systems, Math. Comp. 22 (1968) 28-39, 549-556. [S] J. J. Leendertse, A waterquality simulation model for well-mixed estuaries and coastal seas: Vol. 1, Principles of computation, The Rand Corporation, RM-6230-RC, Feb. 1970. [ 91 J.M. Davis, Water Research Centre, Medmenham, England. Private communication. [lo] S.L. Smith and CA. Brebbia, Finite element solution of Navier-Stokes equations for transient two-dimensional incompressible flow, J. Comp. Phys. 17 (1975) 2355245.