New general feature for v6p0
links from [[:new_version_features|New version features]]
\\
===== v6p0 General =====
===== 1. Telemac-2D and 3D: positive depths on dry zones =====
==== 1.1 Principle ====
The problem of negative depths in Telemac-2D and 3D has always been the price to
pay to have fast and implicit schemes, whereas explicit techniques such as the fi
nite
volume option with kinetic schemes were able to ensure a positive depth, but at a
considerably higher computer time, due to much smaller time steps. To cope with
negative depths, a speci
c smoothing algorithm had been designed, with sometimes
a well-known drawback effect: water slowly creeping above dykes when they were
discretized with too few points. We now present another solution which consists of
limiting the fluxes between points. It is actually a post-treatment which ensures both
mass-conservation and positivity of depth. The continuity equation in the sense of
fi
nite volumes (e.g. as transmitted to Delwaq) is still ensured, the continuity equation
in the sense of fi
nite elements is spoiled because the original velocities are not changed
accordingly to the new depths.
The main idea is summed up here and consists of 3 steps :
* The fluxes between points are computed. We use here the ideas of LeoPostma, already implemented in the interface to Delwaq.
* Starting from depths at time $n$, water corresponding to the uxes are transferred between points, provided that the depth remains positive, otherwise the fluxes are locally limited (fluxes which are not used are kept for a further iteration). This is done in a loop over triangle edges, which can be repeated until there is no more possible water to transfer.
* The remaining fluxes are left over.
We shall now get into the details of the technique. We start from the 2D continuity
equation:
\begin{equation*}
\frac{h^{n+1}-h^n}{\Delta t} + \mathrm{div} (h^{prop}[\theta_u \vec{u}^{n+1}+(1-\theta_u)\vec{u}^n]) = Sce
\end{equation*}
where $Sce$ stands for the discharge sources at points (culverts and so on). This
equation is discretized in the following matrix form:
\begin{equation*}
\frac{m}{\Delta t} (H^{n+1} - H^n) = BM1 \ U + BM2 \ V + RHS
\end{equation*}
(1.1)
where $M$ is the mass-matrix and $\Delta t$ the time-step. $RHS$ is a right-hand side accounting for the boundary fluxes (stemming from an integration by part of the divergence
term) and the source terms stemming from $Sce$. If we use mass-lumping (it will be
mandatory here), this will give for every degree of freedom $i$:
\begin{equation*}
\frac{S_i}{\Delta t} (h_i^{n+1}-h_i^n) + (BM1 \ U + BM2 \ V - RHS)_i = 0
\end{equation*}
(1.2)
where $S_i = \int_\Omega \Phi_i \mathrm{d} \Omega$ is as well the area of the fi
nite volume around point $i$, and
what we call the volume of basis $i$. Options like the wave equation and the speci
fic
treatment of the free surface gradient (keyword "free surface gradient compatibility")
also fi
t within this framework, at the cost of changing $U$ and $V$ into modifi
ed velocities
which may be partially treated as piece-wise constant. These modifi
ed velocities are
denoted $UDEL$ and $VDEL$ in Telemac (because they are used by the interface to
Delwaq).
The quantity:
\begin{equation*}
(BM1 \ U + BM2 \ V - RHS)_i
\end{equation*}
can be interpreted as the flux leaving point $i$. It includes source terms and fluxes at
the boundaries, which are in $RHS$. The terms:
\begin{equation*}
BM1 \ U + BM2 \ V
\end{equation*}
at element level are:
\begin{equation*}
\Phi_i^{el} = - \int_\Omega h \vec{u} . \overrightarrow{\mathrm{grad}}(\Phi_i) \mathrm{d} \Omega
\end{equation*}
(1.3)
and we have explained in references [2] and [3] how to transform them into fluxes
between points, so that the complete continuity equation becomes:
\begin{equation*}
\frac{S_i}{\Delta t}(h_i^{n+1}-h_i^n) + \sum_j \Phi_{ij} + b_i = Sce_i
\end{equation*}
(1.4)
where $b_i$ are the fluxes at the boundaries (denoted $FLBOR$ in Fortran sources). Note
also that term $Sce_i/S_i$ is $SMH$ in Telemac-2D Fortran sources. We have now:
\begin{equation*}
h_i^{n+1} = h_i^n - \frac{\Delta t}{S_i} \left(-Sce_i + \sum_j \Phi_{ij} + b_i \right)
\end{equation*}
(1.5)
==== 1.2 Limiting internal fluxes ====
Let us
first deal with fluxes between points, regardless of other boundary and source
terms. Starting from $h^n$ we want to construct a new depth at time $n + 1$, and the
depth "in construction" is denoted here $\tilde{h}$. In a loop over all segments, we get every
time a specifi
c $i$ and $j$ (apices of the segment), and we would like to apply the formula:
\begin{equation*}
\tilde{h}_i \ \mathrm{replaced} \ \mathrm{by} \ \tilde{h}_i - \frac{\Delta t}{S_i} \Phi_{ij}
\end{equation*}
(1.6)
\begin{equation*}
\tilde{h}_j \ \mathrm{replaced} \ \mathrm{by} \ \tilde{h}_j - \frac{\Delta t}{S_j} \Phi_{ij}
\end{equation*}
(1.7)
but there is a risk of negative $\tilde{h}_i$. If there is a risk, i.e. if $\Phi_{ij} > \frac{S_i \tilde{h}_i}{\Delta t}$, then the flux is limited at least by a factor:
\begin{equation*}
\theta_1 = \frac{S_i \tilde{h}_i}{\Phi_{ij} \Delta t}
\end{equation*}
As a matter of fact there is also a limiting factor linked to the positivity of ehj when ij
is negative, namely we need to have: (), which gives a second coefficient:
\begin{equation*}
\theta_2 = \frac{S_j \tilde{h}_j}{\Phi_{ij} \Delta t}
\end{equation*}
Only one of both limitations at most will have to be enforced, and at the end there
will be a limiting factor $\theta = \mathrm{min}(\theta_1,\theta_2)$ which will be in the range [0,1]. We then do:
\begin{equation*}
\tilde{h}_i \ \mathrm{replaced} \ \mathrm{by} \ \tilde{h}_i - \theta \frac{\Delta t}{S_i} \Phi_{ij}
\end{equation*}
(1.8)
\begin{equation*}
\tilde{h}_j \ \mathrm{replaced} \ \mathrm{by} \ \tilde{h}_j - \theta \frac{\Delta t}{S_j} \Phi_{ij}
\end{equation*}
(1.9)
which ensures the conservation of water, and:
\begin{equation*}
\Phi_{ij} \ \mathrm{replaced} \ \mathrm{by} \ (1 - \theta) \Phi_{ij}
\end{equation*}
which stores in $\Phi_{ij}$ the flux that has not yet been taken into account (it is likely to be
used in the next loop over all segments). After a number of iterations, the situation
remains unchanged, i.e. a criterion like $\sum \mathrm{abs}(\Phi_{ij})$ is no longer decreasing. The
remaining $\Phi_{ij}$ are then left over. They may be transmitted to the interface to Delwaq
so that they are taken into account to form a perfect continuity equation with depths
and fluxes in accordance.
==== 1.3 Boundary and source terms ====
Boundary and source terms are not likely to be interpreted in terms of fluxes between
points. Moreover a sink term may lead to negative depths, which brings in fact a
specific CFL number for the time step. We have applied the following algorithm:
* Step 1: taking into account the source and boundary terms bringing water (the depths are increased)
* Step 2: applying the limitation of internal fluxes described above, this leads to positive depths.
* Step 3: taking into account the source and boundary terms removing water (the depths are decreased)
Step 3 may raise problems, thus a limiting factor of the source terms or boundary
terms may be applied also at this level.
==== 1.4 Extension to 3D ====
**Horizontal fluxes**
The extension to 3D raised a priori little difficulty. When the wave equation option is
used (i.e. the only option left in version 6.0) the shallow water continuity equation is
solved first and we can work on it as it has been described above. The only difference
is that the compatible 2D depth averaged velocity f
ield is not known, as we work on
3D velocity fields and then do the integration on depth at the discrete level. The fluxes
are thus first computed in 3D, then assembled on the vertical. Instead of computing
fluxes with equation 1.3, we compute in 3D:
\begin{equation*}
\Phi_i^{el} = - \int_\Omega . \vec{u} . \overrightarrow{\mathrm{grad}}(\Phi_i) \mathrm{d} \Omega
\end{equation*}
without assembling the element by element fluxes, and add them on the vertical,
to get non assembled fluxes on triangles. For every layer of 3D elements, the contributions of the 6 points of a prism are added on the 3 points of a triangle in the following way:
* 3D points 1 and 4 on 2D point 1
* 3D points 2 and 5 on 2D point 2
* 3D points 3 and 6 on 2D point 3
then we can compute the uxes between 2D points as already explained.
**Vertical velocities**
* A posteriori, a new difficulty appeared. In 3D, the horizontal fluxes (computed
in the transformes mesh) are used to compute the vertical velocities in the
transformed mesh. The horizontal fluxes must then be strictly compatible to the
depth. As we limit only the point to point fluxes and not the original advecting
field, there is an error here that could translate into wrong vertical velocities.
The solution consisted of replacing the element by element horizontal fluxes
by edge by edge horizontal fluxes to compute the internal fluxes (array called
$FLUINT$ in Telemac-3D). The edge by edge fluxes can be limited, provided that
the 2D limiting coefficient has been kept. This 2D limiting coefficient (called
$FLULIM$) is given per 2D segment and is applied for all planes on the vertical.
This is done in subroutine $FLU3DLIM$.
==== 1.5 Domain decomposition in parallelism ====
The algorithm raises an important problem in parallelism with domain decomposition.
During the loop on all segments which changes the depths under construction, these
changes, if done on an interface point, should be immediately transmitted to the
relevant neighbouring sub-domains. This has been considered a too heavy way to
proceed. So far we resorted to the following procedure: the depths of interface points
are merely shared between processors (structure $MESH%FAC%R$, the inverse of the
number of sub-domains a point belongs to, used for parallel dot products, is available
for this). This means that if a point belongs to 2 sub-domains, the loop on segments
will start locally with half the real depth. All sub-domain will ensure the positivity
of their part of depth. The modified depths will then be summed after the loop (this
also ensures the digit-to-digit equality of depths). Tests show that in this way we have
also a convergence of the algorithm, but slower. It is important here that the fluxes
used are not assembled (other sub-domains are ignored when they are built). The
parallel assembly is in fact done when the depths are summed on interface points.
It could be that the algorithm is hindered by the fact that fluxes are not assembled.
For example, on either side of an interface segment fluxes could be opposite and sum
to 0, which is easier to avoid negative depths. In fact any combination of fluxes
that have the same sum when assembled should work, but maybe with a different
efficiency. We have eventually chosen to take the average on either side. This seems
consistant with the fact that we share the depth, and it will become mandatory with
tracers, because the sign of the flux will give the way to do the upwinding. A specific
subroutine "mult_interface_segments" has been designed for this. Tests show that
there is no overcost, which, given that it costs an extra parallel communication, is a
hint that it speeds up the rest of the process.
==== 1.6 Tests and applications ====
=== 1.6.1 The Malpasset dam break ===
**With Telemac-2D:**
Figure 6.1 compares the previous smoothing algorithm and the one presented here,
in scalar mode and in parallel mode with 16 processors. The computation is done on
a Unix HP C3700 with a frequency of 750 MHz and on a Dell Linux machine with 2
processors. For the sake of simplicity the 16 processors have only been simulated as
16 different runs on a single HP workstation.
Table 1.1 gives the computer times of different options and machines with 1 or 2
processors. Kinetic schemes have also been tested, as an alternative to get positive
depths, but the times are much higher, due to the fact that it is an explicit scheme
with a CFL limitation which leads to time steps sometimes 100 times smaller than
the 4 s used on fi
nite element side. On Figure 6.1 a threshold of 1 cm has been chosen
| ^ smoothing ^ flux correction ^ kinetic scheme ^
^ HP C3700 750 MHz | 133 s | 168 s | > 6 hours |
^ Dell (Calibre 5) 1 processor | 60 s | 68 s | 2 h 32 ' 18" |
^ Dell (Calibre 5) 2 processors | 44 s | 50 s | no parallelism |
Table 1.1: computer times with the Malpasset dam break test-case
for displaying depths. This allows a fair comparison of inundation extents. What is
then not seen in the fi
gure at the top (smoothing) is the fact that depths of -0.11 m
are observed in the results. The -1 mm line extends far from the normally flooded
area, due to the number of smoothings applied. Another important point is that
due to truncation errors, the smoothing algorithm starts working also for the water
at rest in the sea. At the end of the computation, a boundary point o¤shore which
should have a depth of 20 m has only 19.9994 m. This is negligible but may lead to
noticeable di¤erences in long term computations. The new algorithm suppresses all
these drawbacks. The mass conservation is excellent in both cases. The water volume
lost with smoothing is -1858.743 m3, and it is 0.298 10,-7 m3 with flux correction. If
we use an exact solver for linear systems, the error becomes respectively 0.596 107
m3 (computer time 217 s) and -0.11 10-6 m3 (computer time 276 s).
**With Telemac-3D**
The behaviour is quite comparable to what is obtained in 2D, the computer times
obtained with the 2-plane case are given in Table 1.2.
| ^ smoothing ^ flux correction ^
^ HP C3700 750 MHz | 769 s | 768 s |
^ Dell (Calibre 5) 1 processor | 341 s | 361 s |
^ Dell (Calibre 5) 2 processors | 238 s | 252 s |
Table 1.2: computer times with the Malpasset dam break test-case
The difference of computer time is not significant.
Again with the smoothing option, the total mass lost at the end of the compu-
tation is -143.8726 m3, whereas with positive depths it is -0.59 10-7 m3, the solver
relative accuracy being 10-6. If we use an exact solver for linear systems, it becomes
respectively -0.29 10-7 m3 (computer time 822 s) and -0.298 10-7 m3 (computer time 850 s). We have thus the same behaviour than in 2D and the explanation is the
following: the flux correction algorithm recomputes the depth with the continuity
equation. In the original continuity equation, which is the wave equation, we have
to solve a linear system with a matrix that contains a Laplacian, due to the implicit
treatment of the free surface gradient stemming from the momentum equation. Then
this free surface gradient is incorporated in the advection field and is thus explicit,
the matrix become a diagonal (mass-lumping is mandatory). The resulting continuity
equation is thus much simpler to solve in an exact way. The conclusion is that when
this flux correction is applied, the continuity equation is exactly solved, whatever the
accuracy asked for the initial linear system. Consequently, even if there is no tidal
flat, the algorithm will have an effect on depth to get an exact continuity equation.
=== 1.6.2 The Wesel river ===
**With Telemac-2D**
This river hydraulics case, a steady state flow, was provided by the BAWKarlsruhe. It
is numerically speaking very tough, with very high Courant numbers (time steps of 2
minutes) and a locally highly refined mesh (17340 elements) including many groynes.
There are 360 steps in the computation. Table 1.3 summarizes the minimum depth
and computer times for 3 techniques, including masking of dry elements. These masked
elements may contain hidden negative depths which reappear in the post-treatment,
as is the case here. Note that the smoothing technique manages to limit the negative
depths at -1 cm. Figure 6.3 shows the water depth in a small part of the domain, the
smoothing technique and the flux correction give very similar results. The difficulty of
the problem shows in the computer time of the flux correction technique. 12 iterations
are necessary here. The minimum depth to be corrected is -0.22 m. The initial sum
of absolute values of all fluxes is 337119.13 m3/s. At the end of the process 6.67 m3/s
only are discarded as generating negative depths.
| ^minimum depth ^ computer time (HP C3700) ^
^ smoothing | -0.01 m | 42 s |
^ flux correction | 0 m | 65 s |
^ masking | -0.085 m | 43 s |
Table 1.3: minimum depth and computer time on the Wesel test-case
**With Telemac-3D**
A remaining problem is that this test-case does not work properly in 3D, unless the
time step is reduced to 50 s instead of 100 s, and with a diffusion of 2 m2/s, whereas the
smoothing depth option works without diffusion. Large horizontal velocities appear
on groynes, i.e. in highly r
efined zones, and we notice a sensitivity to the constant
chosen to change finite element fluxes into finite volume fluxes. This is exemplified
in Figure 6.8 which shows the velocity field near a groyne in red, with the underlying
free surface elevation in colored surfaces. The conclusion is different from the case
with tracers in 2D, where Leo Postma's constant was the best.
==== 1.7 Discussion ====
It is worth noticing that this algorithm may be applied also when there is water
everywhere. In this case we may have also a temporary limitation of fluxes during
the iterations (especially with large Courant numbers) but eventually all the fluxes
should be accepted. A test $\sum \mathrm{abs}(\Phi_{ij}) = 0$ has thus been added to the cases that
stop the iterations.
Another point is that it is very amazing to see that there is no specific CFL
limitation of the process, which is on the contrary a strong limitation of explicit
schemes. The price to pay is in the iterations, but it happens that it is much less
restrictive than real explicit schemes when they want to ensure positive depths. This
fact will occur again in the next chapter.
It seems easy at first sight to extend this technique to 3D, as the continuity equation
is broadly the same. As a matter of fact the smoothing of negative depths was also
used in 3D without extra problems.
Last but not least: the extension to 3D is straightforward as it works on the same
continuity equation.
===== 2. Sisyphe: the positive depths algorithm applied to non erodable beds =====
The Exner equation in Sisyphe, for bed-load transport, is formally similar to the
Saint-Venant continuity equation as it reads:
\begin{equation*}
\frac{\partial Z_f}{\partial t} + \mathrm{div}(\vec{Q}_s) = 0
\end{equation*}
where $Z_f$ is the bottom and $\vec{Q}_s$ the solid discharge. When there are non erodable
beds, a new constraint is that $Z_f$ must not be lower than $Z_r$, the elevation of non
erodable bed. So far non erodable beds were dealt with by an a priori treatment of
$\vec{Q}_s$ ensuring the required property. A limitation factor $g$ for $\mathrm{div}(\vec{Q}_s)$ was first computed, then a second limiting factor $f$ was deduced such that $\mathrm{div}(f \vec{Q}_s) < g \mathrm{div}(\vec{Q}_s)$, and the form $\mathrm{div}(f \vec{Q}_s)$ was eventually used in Exner equation, thus ensuring mass
conservation (see reference [6]). This was however a rather tedious procedure. We
now consider that $Z_f$ - $Z_r$ is an available layer of sediment, that must not become
negative, the new equivalent equation is:
\begin{equation*}
\frac{\partial (Z_f-Z_r)}{\partial t} + \mathrm{div}(\vec{Q}_s) = 0
\end{equation*}
with the same constraint than the water depth for $Z_f$ - $Z_r$. In the variational
formulation the term $\mathrm{div}(\vec{Q}_s)$ is integrated by parts and thus split into boundary
fluxes and internal fluxes. These internal fluxes can be tranformed into point to
point fluxes like in previous chapter, then the problem is exactly like moving water
between points, the sediment replaces water. For graded sediment, each class has its
own layer (which in Sisyphe is a product AVAIL * ELAY) and it leads to a series of
problems where every layer must remain positive. With a saturated bed load equation,
discarding some uxes is even more natural than with the Saint-Venant continuity
equation. As a matter of fact, the solid discharge is only a potential discharge that
assumes that sediment is available for movement. When there is no sediment, it is
normal to ignore that flux. Iterations are necessary, just in case sediment is brought
within the time step to a point that was previously bare. Figure 6.11 compares the
old and the new technique on a Sisyphe schematical test-case. Sediment in a square
hole is washed away by the flow. No limitation is prescribed on the bed slope here so
the heap of sediment is not very physical, but shows the similar behaviour of both
techniques which, given their very different approach, is a good cross-validation. In
this example the new technique is much faster.
===== 3. Telemac-2D and 3D: a new advection scheme designed for tidal flats =====
==== 3.1 Principle ====
The upwind explicit finite volume scheme was the only one in Telemac-2D to ensures
mass conservation of tracers in the sense of depth-averaged concentrations or temperatures.
This scheme could not be used so far with tidal flats because there is a
division by the depth in the derivation, and the CFL number tends to infinity on dry
zones. We had then to resort to a masking of dry elements, along with a clipping of
depth, which was not mass-conservative. The flux correcting technique presented in
the first chapter leads in fact straightforwardly to a new advection scheme which is
not sensitive to dry zones.
In the reference [3] we have presented the upwind explicit finite volume advection
scheme. With the same notations as in Chapter 1, the new concentrations at a point
$i$ were given by the formula:
\begin{equation*}
C_i^{n+1} = (1+ \frac{\Delta t}{h_i^{n+1}S_i} \sum_{negative \ \Phi_{ij}} \Phi_{ij})
C_i^n - \frac{\Delta t}{h_i^{n+1}S_i} \sum_{negative \ \Phi_{ij}} C_j^n \Phi_{ij}
\end{equation*}
(3.1)
The monotonicity condition was that:
\begin{equation*}
\Delta t < \frac{h_i^{n+1}S_i}{\displaystyle{\sum_{negative \ \Phi_{ij}} |\Phi{ij}|}}
\end{equation*}
(3.2)
which may turn into 0 on dry zones. The reason is that the inputs and outputs on
one point may result in a negative or zero depth, while the quantity of tracer may not
be zero. The fundamental reason is that all fluxes to and from a point are considered
at the same time. Imagine now that we do it edge by edge, thus considering only two
points at a time. Let us number these points 1 and 2 and let us assume that the
flux $\Phi_{12}$ is positive, i.e. the water goes from point 1 to point 2. The conservative
tracer equations of both points read simply (we omit boundary and source terms for
simplicity):
\begin{equation*}
\frac{S_1}{\Delta t}(H_1^{n+1}C_1^{n+1}-h_1^n C_1^n) + \Phi_{12} C_1^n = 0
\end{equation*}
(3.3)
\begin{equation*}
\frac{S_2}{\Delta t}(H_2^{n+1}C_2^{n+1}-h_2^n C_2^n) + \Phi_{12} C_1^n = 0
\end{equation*}
(3.4)
$C_1^n$ appears in both equations because the flux goes from 1 to 2 (upwind scheme).
If we also treat the continuity equation only for this edge (what we do in the flux
limiting algorithm), we have also:
\begin{equation*}
H_1^{n+1} = h_1^n - \frac{\Delta t}{S_1} \Phi_{12}
\end{equation*}
\begin{equation*}
H_2^{n+1} = h_2^n - \frac{\Delta t}{S_2} \Phi_{12}
\end{equation*}
It then turns out that equations 3.3 and 3.4 become simply:
\begin{equation*}
C_1^{n+1} = C_1^n
\end{equation*}
(3.5)
\begin{equation*}
C_2^{n+1} = \frac{h_2^n}{h_2^{n+1}} C_2^n + (1 - \frac{h_2^n}{h_2^{n+1}})C_1^n
\end{equation*}
(3.6)
In this context there is no risk of division by 0 because we started from a positive
depth $h_2^n$ which was increased by the positive quantity $\frac{\Delta t}{S_2} \Phi_{12}$.
$\frac{h_2^n}{h_2^{n+1}}$is then in the range [0,1]. The positivity and monotonicity of tracers is ensured even on dry zones (in case of zero depth the concentrations remain unchanged). This edge by edge
treatment (only a few lines of Fortran) may be inserted within the previous tidal flats
algorithm.
==== 3.2 Domain decomposition in parallelism ====
We have mentioned in Chapter 1 that the internal uxes were not assembled in parallel,
but could be averaged on interfaces, so that the upwinding information is the
same on either sides. This trick is used here for the tracers. However it is not the
only thing to do. Let us assume that a point is shared between 3 processors, $a$, $b$ and
$c$. It thus exists in 3 locations in memory, with concentrations $C_a$, $C_b$, $C_c$ and depths
$h_a$, $h_b$, $h_c$. At the end of the parallel communication the 3 depths are added so that
every processor gets $h_a + h_b + h_c$. The relevant conservative merging of concentrations gives the common value of $\frac{C_a h_a + C_b h_b + C_c h_c}{h_a + h_b + h_c}$ . This is obtained by creating an
array containing the product $Ch$, and running the parallel communication that adds
contribution of interface points (PARCOM with option 2). The same is done with
$h$ and we then do the division. This is a case where we have a division by a depth,
which is done only if the depth is strictly positive.
==== 3.3 Extension to Sisyphe ====
In Sisyphe, the velocity
eld for advection may be corrected by a factor $F$ in the
range [0,1] to account for the fact that the suspended sediment is concentrated near
the bottom, where the velocity is smaller than the depth-averaged velocity. A special
care is then necessary: the fluxes for computing the positive depths must be done
with the original velocity
eld, which is consistent with the continuity equation. Then
the fluxes to consider for the tracer advection must have the coefficient $F$. We have
then to come back to the original tracer equations 3.3 and 3.4, which now read:
\begin{equation*}
\frac{S_1}{\Delta t}(h_1^{n+1} C_1^{n+1} - h_1^n C_1^n) + F \Phi_{12} C_1^n = 0
\end{equation*}
(3.7)
\begin{equation*}
\frac{S_2}{\Delta t}(h_2^{n+1} C_2^{n+1} - h_2^n C_2^n) + F \Phi_{12} C_1^n = 0
\end{equation*}
(3.8)
if we assume that $\Phi_{12}$ is positive and apply upwinding also on the $F$ coefficient. It
gives now:
\begin{equation*}
C_1^{n+1} = \left[ \frac{h_1^n}{h_1^{n+1}} + F_1 (1 - \frac{h_1^n}{h_1^{n+1}}) \right] C_1^n
\end{equation*}
(3.9)
\begin{equation*}
C_2^{n+1} = \frac{h_2^n}{h_2^{n+1}} C_2^n + F_1 (1 - \frac{h_2^n}{h_2^{n+1}}) C_1^n
\end{equation*}
(3.10)
This solution does not ensure monotonicity of point 1, which is a known drawback of
using a corrected velocity fi
eld. It does not seem practicable here as the interpolation
coefficients may tend to infi
nity. A solution could consist of dealing with erosion and
deposition at the same time, in order to get a more physical behaviour. Advection
would be done on the depth-integrated quantity of sediment $hC$. Then, depending
on the final depth, it would be decided if this quantity would remain in the water or
would deposit.
==== 3.4 Extension to 3D ====
Once fluxes between points are known (this is done in the next chapter for a finite
volume scheme) there is absolutely no difference between 2D and 3D. In 2D quantities
like $S_1 h_1^{n+1}$ and $S_2 h_2^{n+1}$ are volumes of water carried by points 1 and 2. In 3D the
volumes carried by points will be simply the integral of test functions (which have
been purposedly called VOLUN and VOLU in Telemac-3D). Note that the sum of all
these volumes is the integral of 1 over the whole domain, hence the total amount of
water. At the beginning of a computation these integrals are VOLUN, when all the
fluxes between points have been transferred they are equal to VOLU. The algorithm
in 3D is otherwise exactly the same as in 2D. This new advection solver may be used
on the velocities, either in 2D or 3D. In this case the continuity equation is not yet
done and the new depth will only be compatible with the advection velocity at the
beginning of the time step, hence it will be only a predictor of the final depth.
===== 3.5 Tests ====
==== 3.5.1 Tracer in the bridge piers test-case ====
Our first test has no tidal flats at all. It is just a comparison of advection schemes to
assess their numerical diffusion. We use the bridge piers test-case and enter a tracer
with a value of 1 at 3 points in the entrance. The time step is 0.8 s and there are 100
steps in the computation. There is no diffusion. Six different results are plotted on
Figure 6.2. On the right are displayed the results given by the upwind explicit finite
volume scheme, the method of characteristics, and the Positive Streamwise Invariant
distributive scheme. This latter scheme is reputed mass-conservative, but not here
in the context of shallow water equations, because it is applied to the concentration,
which is not the conservative variable. On the left 3 variants of the present scheme,
depending on the choice of the constant to get the element fluxes (see discussion
in reference [3] page 33 and in reference [2]). It happens that the choice of this
constant. Only with the original solution introduced by Leo Postma do we get a
correct advection scheme, comparable to other schemes. The conservation of mass
(of water and tracer) is ensured at the accuracy of the machine (provided that direct
solvers are chosen). It can be noted that the method of characteristics and the PSI
scheme are less diffusive.
==== 3.5.2 Thermal plume in tidal conditions ====
We study here a thermal plume in 2 dimensions, i.e. with depth-averaged temperature,
on the small domain. The boundary conditions, subjected to tidal conditions, are given
by a larger model. For checking the conservation of heat, the diffusion is removed. As
a matter of fact this step is actually applied on the temperature and is conservative for
this variable, which is not the integral of temperature on the vertical (which is the real
conservative variable). The hot water is released by 16 source points of 10.5 m3/s each
with a velocity of 1.51 m/s towards West, and an increment of temperature of 10.6ºC.
There is no exchange with atmosphere. The computation consists of 20000 steps of 50
s each, i.e. 11 days 13 hours 46 mn and 40 s. The total heat (temperature multiplied
by volume) entered in the domain is 0.17808 1010. When the flow is entering the
domain an increment of temperature of 0 is prescribed. At low tide there are dry zones
everywhere in the domain, including on the boundaries. Figure 6.4 shows the thermal
plume at the end of the computation, obtained with the method of characteristics,
the PSI scheme, and the present edge by edge flux correcting scheme. Characteristics
and flux correction are in broad agreement, though the latter appears to be more
diffusive. The PSI schemes overestimates temperatures. The figures of Table 3.1 give
more explanations on the comparison. CPU time is given on a HP C3700 workstation.
Both the PSI scheme and characteristics create an excess of about 20% of heat, but
this excess gets out of the domain in the case of characteristics. It is worth noticing
for further studies that the quantity of heat exiting the domain is here not supposed
to re-enter. The very small loss of heat of the new scheme is however probably a
remaining mistake in theory. It occurs only when the types of boundary conditions
are changed from input to output and vice-versa. It is also worth noting that, to avoid
loss or gain of heat, the Dirichlet boundary conditions of temperature are discarded.
Only the correct fluxes are considered. The values observed at boundaries may thus
appear to be slightly different from prescribed values (this point is also very important
in the Berre lake study).
| ^ CPU time ^ min. depth ^ % heat lost ^
final heat ^ exited heat ^
^ characteristics | 5116 s | -0.0114 m | -19.3 | 0.467 109 | 0.166 1010 |
^ PSI scheme | 5623 s | -0.0114 m | -17.7 | 0.836 109 | 0.126 1010 |
^ flux correction | 7202 s | 0 m | -7 10-11 | 0.646 109 | 0.113 1010 |
Table 3.1: comparing advection schemes on a thermal plume test case
On a 2-core Dell Linux workstation, the CPU time in parallel, for the 3 advection
chemes, is summarized in Table 3.2. As the method of characteristics is already used
for the advection of depth in the momentum equation, it is obviously the most efficient
technique here in terms of computer time, as the advection of tracer is only an extra
interpolation. It is also worth noticing that the flux correction procedure is done
twice, for simplicity of implementation: once for getting positive depths, once again
for the tracer. It is clear by comparing Table 3.1 and Table 3.2 that the efficiency
of algorithms highly depend on the machine (which is a combination of architecture
+ compiler). On Dell in parallel the PSI scheme and the new scheme have the same
cost.
| ^ CPU time with 2 processors on Dell Linux ^
^ characteristics | 1286 s |
^ PSI scheme | 1495 s |
^ flux correction | 1494 s |
Table 3.2: comparing advection schemes in parallel
==== 3.5.3 Dam breaks ====
The 1D dam break exact solution will enable us to show the interest of the new
scheme when applied to the advection of velocities. Figure 6.7 compares the free
surface elevation obtained with the method of characteristics (in red) with the new
scheme (in green). The exact solution is in blue. It appears clearly that the method of
characteristics is hindered by tidal flats. In this ideal case, the wave front progresses
only with advection. There is no propagation velocity on the front because the depth
is zero.
We then tested the Malpasset dam break in 3D, with 2 planes on the vertical. As
in the previous case, the new finite volume advection solver is applied to velocities.
Figure 6.6 compares the depths after 40 mn, with the method of characteristics or
with the new scheme (which is scheme option 9 for key-words SCHEME FOR ADVECTION OF VELOCITIES or SCHEME FOR ADVECTION OF TRACERS). The
surprise is the difference which does not appear so large in 2D, but this difference is
in favour of the new scheme, since it is known that the wave reached the sea before 40
mn. As in the 1D example this may be explained by the fact that the method of backward characteristics is not very good on tidal flats, since no characteristic on a tidal flat (thus with no velocity) will go backward to get information on the coming wave.
The difference with 2D results was investigated, and an explanation could be found
and is explained here: in 2D the friction terms in the non conservative momentum
equation are multiplied by a factor $1/h$ where $h$ is the depth. This depth was so far
taken as the nodal value. In 3D there is no obvious $1/h$ factor but it implicitly hidden
in a factor:
\begin{equation*}
\frac{\int_{\Omega 2D} \Phi_i \ \mathrm{d}(\Omega 2D)}{\int_{\Omega 3D} \Phi_i \ \mathrm{d}(\Omega 3D)}
\end{equation*}
which stems from the variational formulation of the boundary terms. The difference is that when $h$ is zero in 2D friction becomes infinite, velocity is cancelled and
it reduces the wave celerity. When h is zero in 3D, the denominator may not be zero
if one neighbour of the given point has a depth, because the volume associated with
the corresponding test function will not be zero. In this case the friction will not be
infinite and this will ease the wave propagation. This is exempli
ed in Figure 6.9,
where the effect of the new advection scheme is also clear.
==== 3.6 Discussion ====
The new scheme is perfectly mass-conservative in 2D and 3D, it ensures monotonicity
and is stable on dry zones. Further tests are necessary to assess its numerical di¤usion.
A possible drawback is a sensitivity to the mesh, probably by construction because
we use the edges as a way of transit of tracer. An obvious odd property is that the
result certainly depends on the numbering of edges. A necessary development before
other quantitative tests would be adding the diffusion step in the depth-averaged
conservative concept. It would simply consist of evaluating diffusion terms as an
extra advection.
===== 4. Telemac-3D: a finite volume advection solver =====
The main idea is to design an advection solver close to what is done with the distributive
scheme, i.e. with a splitting of time to ensure the CFL condition, with the
only modification that the fluxes taken into account are not the finite element fluxes,
but rather the finite volume point to point fluxes which are currently built for the
Delwaq interface and for the treatment of negative depths in 3D. In 3D and with distributive
schemes (full explanations on distributive schemes are given in reference [1]
from page 183 on) we just recall here that we end up in a discretization of advection
equation in the form:
\begin{equation*}
S_i (\frac{C_i^{n+1}-C_i^n}{\Delta t}) = -\beta_i \Phi_T
\end{equation*}
(4.1)
Where $S_i$ is the integral of the test function of point $i$, $\beta_i$ is a distribution coefficient, $\Delta t$ the time step, $C_i^{n+1}$ the value of function $C$ at point $i$ after advection, $C_i^n$ the
value of function $C$ at point $i$ before advection. $\Phi_T$ actually represents $S_i(\vec{u} . \mathrm{\overrightarrow{grad}}(C))$
where \vec{u} is the advecting fi
eld. The only thing to know here is that $\beta_i \Phi_T$ is eventually
put in the form:
\begin{equation*}
\beta_i \Phi_T = \sum_j \lambda{ij} (C_i - C_j)
\end{equation*}
(4.2)
where the coefficients $\lambda_{ii}$ are zero and all coefficients $\lambda_{ij}$ are positive or zero. The
sum is done on all the neighbouring points of point $i$, i.e. all the points that belong
to an element containing $i$.
The final value of function $C$ at point $i$ is thus:
\begin{equation*}
C_i^{n+1} = C_i^n \left( 1 - \sum_j \lambda{ij} \frac{\Delta t}{S_i} \right) + \sum_j \frac{\Delta t}{S_i} \lambda_{ij} C_j^n
\end{equation*}
(4.3)
Classically, it is said then that the monotonicity criterion consists of ensuring that all the coefficients of $f_i^n$ and $f_j^n$ be positive, which yields, given the fact that all $\lambda_{ij}$ are positive or zero:
\begin{equation*}
\Delta t \le \frac{S_i}{\displaystyle{\sum_j \lambda_{ij}}}
\end{equation*}
(4.4)
In finite volumes with fluxes between points, the corresponding equation would
be:
\begin{equation*}
C_i^{n+1} = \left(1 - \frac{\Delta t}{S_i} \sum_j \mathrm{max}(\Phi_{ij},0) \right) C_i^n + \sum_j \frac{\Delta t}{S_i} \mathrm{max}(\Phi_{ij},0)C_j^n
\end{equation*}
assuming that we start from element by element fluxes in the form:
\begin{equation*}
\Phi_i^{el} = \int_\Omega \vec{u} . \overrightarrow{\mathrm{grad}}(\Phi_i) \mathrm{d} \Omega
\end{equation*}
this explains the difference of signs with what was done in 2D where we start from
fluxes in the form:
\begin{equation*}
\Phi_i^{el} = \int_{\Omega 2d} h \vec{u} . \overrightarrow{\mathrm{grad}}(\Phi_i) \mathrm{d} (\Omega 2d)
\end{equation*}
which are fluxes leaving points. This leads us to the following CFL condition:
\begin{equation*}
\Delta t \le \frac{S_i}{\displaystyle{\sum_j \mathrm(max)(\Phi_{ij},0)}}
\end{equation*}
(4.5)
Implementation should be the same provided that we change everywhere $\lambda_{ij}$ by
$\mathrm{max}(\Phi_{ij},0)$.
In fact, the distributive scheme implementation in subroutine murd3d.f in Telemac
only uses $\displaystyle{\sum_j \lambda_{ij}(C_j^n - C_i^n)}$ which is array TRA02 and $\displaystyle{-\sum_j \lambda_{ij}}$ which is put in array DB. A finite volume upwind scheme can then be easily implemented if we take:
\begin{equation*}
DB = -\sum_j \mathrm{max}(\Phi_{ij},0)
\end{equation*}
and
\begin{equation*}
TRA02 = \sum_j \mathrm{max}(\Phi_{ij},0)(C_j^n - C_i^n)
\end{equation*}
Starting from the coefficients $\Phi_i^{el}$ (six per element), the point to point coefficients $\Phi_{ij}$ are obtained as is done in the interface to Delwaq, i.e.:
1) assembling element by element fluxes on points: horizontal fluxes will be a sum
of two contributions (from the upper and lower layer, except on the bottom and the
free surface).
2) changing the horizontal fluxes of points 1, 2, and 3 in the prism into point to
point fluxes.
Points 1 and 2 are done in a subroutine called fluxvfleo_3d, which calls the 2D
subroutine fluxvfleo.f
3) cancelling all crossed fluxes (between points 1 and 5, 1 and 6, 2 and 4, 2 and 6,
3 and 4, 3 and 5)
4) finding the vertical fluxes that solve the continuity equation. Comparing what
is done in subroutine tridw2.f and in subroutine tel4del.f, we find that the vertical
fluxes are such that:
\begin{equation*}
\Phi_{ij} = -\Delta z \ W^* \int_{\Omega 2d} \Phi_k \mathrm{d}(\Omega 2d)
\end{equation*}
where $\Phi_{ij}$ is the vertical flux between to points $i$ and $j$ located in the same prism,
on the vertical of the 2D point $k$. On a total of 15 possible fluxes in the prism, we
thus retain 9 non-zero terms, which is exactly what is done also with the N-scheme,
with the same vertical fluxes. Both schemes are actually very close, though different,
and it was very easy to add the new finite volume scheme as a variant of distributive
scheme, in subroutine murd3d.f. Figure 6.5 compares the 3 available options on the
test-case of the lock-exchange.
===== 5. FE and FV fluxes: Leo Postma scheme is the N-scheme =====
The way to get point to point uxes from
nite element uxes has already been
discussed in reference [3] (page 29). The problem is not well posed and a constant
must be introduced. In our tests the constant chosen by Leo Postma seemed to give
the best results on the rotating cone test. More recently it became obvious that the
distributive schemes were another way to get the point to point uxes and the N-
scheme was tested. Though the formulas are at
rst sight very di¤erent, it happens
that the Leo Postma scheme is exactly the N-scheme. We show it hereafter. We keep
the notations of reference [3], recalled in the sketch below (a negative 1 means that
point 1 loses water).
The original Leo Postma formula can be implemented for every element in the
following Fortran-like form (the assembly on segments has been removed to allow a
better understanding):
\\
IF(ABS($\Phi_1$).GE.ABS($\Phi_2$).AND.ABS($\Phi_1$).GE.ABS($\Phi_3$)) THEN \\
$\Phi_{12}$ = $\-\Phi_2$ \\
$\Phi_{23}$ = 0 \\
$\Phi_{31}$ = $\Phi_3$ \\
ELSEIF(ABS($\Phi_2$).GE.ABS($\Phi_1$).AND.ABS($\Phi_2$).GE.ABS($\Phi_3$)) THEN \\
$\Phi_{12}$ = $\Phi_1$ \\
$\Phi_{23}$ = $-\Phi_3$ \\
$\Phi_{31}$ = 0 \\
ELSE \\
$\Phi_{12}$ = 0 \\
$\Phi_{23}$ = $\Phi_2$ \\
$\Phi_{31}$ = $-\Phi_1$ \\
ENDIF \\
The N-scheme implementation is now the following:
\begin{equation*}
\Phi_{12} = \mathrm{ MAX(MIN(-\Phi_1,\Phi_2),0.D0) - MAX(MIN(-\Phi_2,\Phi_1),0.D0) }
\end{equation*}
\begin{equation*}
\Phi_{23} = \mathrm{ MAX(MIN(-\Phi_2,\Phi_3),0.D0) - MAX(MIN(-\Phi_3,\Phi_2),0.D0) }
\end{equation*}
\begin{equation*}
\Phi_{31} = \mathrm{ MAX(MIN(-\Phi_3,\Phi_1),0.D0) - MAX(MIN(-\Phi_1,\Phi_3),0.D0) }
\end{equation*}
Figure 5.1: fluxes from and between points
To prove the equivalence we first remark that we can always consider that point
1 has a positive $\Phi_1$ (we could otherwise change the local numbering of points without
changing the result). Then we are left with 3 cases (because the sum of all fluxes is
0):
case 1: $\Phi_1$ > 0 and $\Phi_2$ > 0 and $\Phi_3$ < 0 \\
case 2: $\Phi_1$ > 0 and $\Phi_2$ < 0 and $\Phi_3$ < 0 \\
case 3: $\Phi_1$ > 0 and $\Phi_2$ < 0 and $\Phi_3$ > 0 \\
It is easy (but tedious, headache guaranteed!) to check that both methods will
give:
case 1: $\Phi_{12}$ = 0 and $\Phi_{23}$ = $\Phi_2$ and $\Phi_{31}$ = -$\Phi_1$ \\
case 2: $\Phi_{12}$ = -$\Phi_2$ and $\Phi_{23}$ = 0 and $\Phi_{31}$ = $\Phi_3$ \\
case 3: $\Phi_{12}$ = $\Phi_1$ and $\Phi_{23}$ - $\Phi_3$ and $\Phi_{31}$ = 0 \\
Note that the N-scheme was designed to get positive fluxes, which is the case here.
The N-scheme uses 6 MAX functions and 6 MIN functions.
The Leo Postma method can be reduced to 3 ABS functions and 2 or 4 comparisons
.GE., it is cheaper, but the proof given here could be translated into: \\
\\
IF($\Phi_1$.GE.0) THEN \\
IF($\Phi_2$.GE.0) THEN \\
This is case 1 \\
ELSE \\
IF($\Phi_3$.GE.0) THEN \\
This is case 3 \\
ELSE \\
This is case 2 \\
ENDIF \\
ENDIF \\
ELSE \\
IF($\Phi_2$.GE.0) THEN \\
IF($\Phi_3$.GE.0) THEN \\
This is case 4 \\
ELSE \\
This is case 5 \\
ENDIF \\
ELSE \\
This is case 6 \\
ENDIF \\
ENDIF \\
\\
This way of writing the tests limits to 3 comparisons .GE. in the worst case
(probably not worth the effort).
===== 6. Telemac-3D: smashed elements now possible =====
Up to version 5.9, only smashed elements on tidal flats or dry zones were treated. The
fact that such elements had no volume precluded the use of finite-volume like numerical schemes, hence only the method of characteristics could be applied for advection.
It appeared that most of the treatments done in the case of tidal flats could be applied
also in the case of smashed elements with water above. There is indeed a risk of such
a situation when the generalised sigma transformation is used, and a constant eleva-
tion of planes is prescribed. When the bottom goes above the prescribed elevation of
a plane, a number of elements is smashed and this caused a crash of computations.
To avoid this the parameter DISMIN in subroutine CALCOT guaranteed a minimum
height in the elements. From version 6.0 on, it is now possible to have DISMIN=0.
To be more precise there is now a DISMIN_BOT and a DISMIN_SUR, respectively
for bottom and free surface, and DISMIN_BOT may be 0. The key modification for
achieving is an array of integers IPBOT, of size NPOIN2 (number of points in the 2D
mesh), giving the rank of the last layer with no height. if NPLAN is the number of
planes on the vertical, we have for a 2D point I: \\
IPBOT(I)=0 if all layers are normal (height greater than 1 mm) \\
IPBOT(I)=3 if plane 3 and 4 are closer than 1 mm and distance between plane 4
and plane 5 is greater than 1 mm. \\
IPBOT(I)=NPLAN-1 if all the planes are smashed (case of tidal flats). NPLAN
is the number of planes. \\
This new array allowed a number of speci
c treatments listed below:
* Friction is applied at the level IPBOT(I)+1 and all points below.
* In diffusion all points below the real bottom are treated as Dirichlet points, with the previous value as prescribed Dirichlet value. After solving the linear system the points below the real bottom are given the value of the real bottom.
* The Poisson equation for the non-hydrostatic pressure is treated in the same way as diffusion.
A general principle is that points with the same elevation on a vertical must even-
tually have the same physical value, so that no artificial infinite gradient is created.
1 mm is an arbitrary but reasonable choice and once it is done all the tests are on
IPBOT. Figure 6.10 is an example of flow over a bump with an intermediate plane
(number 4) imposed at the elevation -0.2 (it appears at elevation -1 on the picture
due to a distortion of 5 on the vertical). On the bump planes 1, 2, 3 and 4 are
superimposed. Because the velocities of superimposed points are the same, only a
single vector appears on the
gure. For this case the new advection scheme for tidal
flats is used and can deal with elements without volume. Mass conservation of salt is
verified and monotonicity (here salinity between 30 and 40 g/l) is strictly ensured.
====== Bibliography ======
[1] HERVOUET J.-M., Hydrodynamics of free surface flows, modelling with the finite element method. Wiley & sons. 2007 \\
[2] POSTMA L., HERVOUET J.-M., Compatibility between finite volumes and finite elements in solutions of shallow water and Navier-Stokes equations. International Journal for Numerical Methods in Fluids. 53(9), pp.1495-1507. 2007. \\
[3] HERVOUET J.-M., PHAM C.-T.: Telemac version 5.7, release notes. Telemac-2D and Telemac-3D. 2007 \\
[4] HERVOUET J.-M., RAZAFINDRAKOTO E., VILLARET C.: Telemac version 5.8, release notes. Telemac-2D, Telemac-3D and Sisyphe. 2008 \\
[5] HERVOUET J.-M.: Telemac version 5.9, release notes. Bief, Telemac-2D, Telemac-3D and Sisyphe. 2009 \\
[6] HERVOUET J.-M., MACHET C., VILLARET C.: dealing with rigid beds in saturated bed load transport equations. Proceedings of Coastal Engineering 2003, Cadiz, 23-25 June 2003, Spain. 2003 \\