open TELEMAC-MASCARET The mathematically superior suite of solvers

Hello,

we have now a number of morphodynamic river models applying lower resolution meshes, computing in parallel. Unfortunately the errors with which Telemac users learned to live with in the past, have started to affect us again. The results obtained in serial and with different number of mesh partitions differ. The errors accumulate with the time of simulation. Reproducibility of the results (especially after restarts) is haphazard.

A closer look reveals that the cause is one of the still unresolved problems in Telemac, the parallel "tidal flats" or "wetting & drying" algorithms. These algorithms (three versions) require that a treated element on the shoreline has full information about the state ("wet/dry", waterlevels) in all neighbouring elements. This is not the case in parallel runs.

If the element in question lies not only on the shoreline, but also on the boundary between partitioned meshes, it does not become the full information on the elements behind the partition boundary ("interface"), because only the nodal values are exchanged via communication between processors. In the result the "wet/dry" state in elements where the shoreline and the partition interface can differ compared to the case when the interface is not there.

Additionally, the wet/dry states are steered by threshold values (minimum depth). How these thresholds are reached can differ depending on the number of partitions, because although the iterative slovers for the waterlevel have the same accuracy, this accuracy is reached in a different way and the values on the interface nodes can differ slightly anyway.

The result are isolated areas with differences between waterlevels and currents nearby the shoreline depending on the number of processors one uses (and which processors as well). If you compute stationary flows, these errors remain the same and you might mentally ignore them if you are rather interested in an overall, global sight of the waterlevels and not in the current details.

HOWEVER, if you compute the bed load transport, the errors accumulate during the computation. Then, if larger erosion or sedimentation appears, the wetting and drying algorithms are asked to adjust the shoreline again -- and the artificial errors in parallel grow, grow, grow. And we have another factor influencing the error: the simulation length in time.

This inconsistency has been reported -- so far as I remember -- more than ten years ago, there were some cosmetic improvements, but rather concerning treatment of some specific exceptional cases in wetting & drying algorithms than the parallel problems. The users have learned to live with it, refining the meshes in the problematic areas, making the meshes more regular there, computing with increased accuracy, manipulating the minimum depth threshold, cutting forelands off and switching all these tidal flats algorithms off.

I assume the problem can affect also other (newer?) algorithms depending on the correct or consistent treatment of the shoreline, like the "positive depths"...

The problem has been thoroughly discussed in the past. Modifying of the existing algorithms for parallelism has been excluded as too bothersome and awkward. The solution applying non-linear Newton iterations with the solver for the free surface has been excluded as well, because in case of Telemac-2D it would be inefficient (V. Casulli, "A high-resolution wetting and drying algorithm for free-surface hydrodynamics", Int. J. Numer. Meth. Fluids 2009; 60:391–408). The resigned users have stopped complaining, anyway.

QUESTION: Are there presently any ideas, how to improve the wetting and drying algorithms in Telemac? Especially coupled with Sisyphe (the errors accumulate!)

Best regards,
jaj

Hello Jacek,

I think you mention the treatment of tidal flats with masking of elements, which for me should be avoided for the reasons you give, and also because mass conservation is not exactly ensured. We do not remove it because every time you suppress an obsolete option there is always someone who complains that he/she was using it successfully, but if you insist we'll remove it. The other option (option 1) should be OK and if not please send the corresponding case so that we have a look at it.

Actually we have a PhD student working in the university of Perpignan on reproducibility. It is sufficiently advanced to know that there are no more parallel mistakes in Tomawac (there is now an option in version 7.0 which guarantees the same digits in serial and parallel). In this case only the finite element assembly could introduce truncation errors and this can be solved with two solutions: using integers in some sums, or what is called "compensated algorithms".

Officially the only left reasons for having differences between serial and parallel is the truncation errors in finite element assembly and in the computation of dot products. The PhD main task is to deal with this causes of errors and track finite element assembly and dot products in the different modules. Once this is secured we shall decide that all remaining differences are a real bug that must be looked at.

Note also that a recent important improvement was to ensure all digits equality in the streamline module. This was very fruitful as the last remaining reasons for differences did not affect at all the results but only the computer time.

Of course any truncation error may lead in the end to large differences, this is also due to our equations and we must live with this. There will always be a high sensitivity of cases with von Karman eddies and in that sense there is no reproducibility in nature.

With best regards,

Jean-Michel Hervouet

Hello JMH,

the relatively coarse river model I have looked at more thoroughly was built with Telemac2D coupled with Sisyphe v6p3r2. Parameters for the steering of numerical algorithms (and, this is our "favorite parameter set" I presume), I also leave the commented-out options:

/

---

/ Numerische Parameter
/

---

SOLVER ACCURACY = 1.E-06
PRECONDITIONING = 2
TYPE OF ADVECTION = 14;5
/5="SCHEMA PSI CONSERVATIF"
/14="EDGE-BASED N-SCHEME"
/ NUMBER OF SUB-ITERATIONS FOR NON-LINEARITIES = 3
MAXIMUM NUMBER OF ITERATIONS FOR SOLVER = 300
SUPG OPTION = 2;0 / 2;2
SOLVER = 1
SOLVER OPTION = 3
TREATMENT OF THE LINEAR SYSTEM = 2
DISCRETIZATIONS IN SPACE = 11;11
MASS-LUMPING ON H = 1
TURBULENCE MODEL = 2
/ TURBULENCE MODEL FOR SOLID BOUNDARIES = 2
/ ROUGHNESS COEFFICIENT OF BOUNDARIES = 0.001
/ INFORMATION ABOUT K-EPSILON MODEL
/ VELOCITY DIFFUSIVITY = 1.E-06
TIDAL FLATS = YES
OPTION FOR THE TREATMENT OF TIDAL FLATS = 1
COMPATIBLE COMPUTATION OF FLUXES = YES
CONTINUITY CORRECTION = YES
TREATMENT OF NEGATIVE DEPTHS = 2 / geht nur mit tidal options 1
MATRIX STORAGE = 3
MATRIX-VECTOR PRODUCT = 1
IMPLICITATION FOR DEPTH = 1.0 / 0.8
/ IMPLICITATION FOR VELOCITY = 1.0 / 0.6
OPTION FOR THE DIFFUSION OF VELOCITIES = 2
/ MINIMUM VALUE OF DEPTH = 0.05
H CLIPPING = NO
FREE SURFACE GRADIENT COMPATIBILITY = 0.8
/

---

In this version "dictionary" we have:

AIDE1 = 'Used if TIDAL FLATS is true
1 : EQUATIONS SOLVED EVERYWHERE WITH CORRECTION ON TIDAL FLATS
2 : DRY ELEMENTS FROZEN
3 : LIKE 1 BUT WITH POROSITY (DEFINA METHOD)'

So far as I remember everyone here uses option "1", with "2" we have had bad experiences in the past (flow blocking), at least for coarse meshes and there were some issues with "3", porosity a la Defina.

I also suspect that some of that fancy advection schemes "14;5" can be affected, I have limited experience here, or we have overlapping inconsistencies.

I will ask here if we can send you the particular model, so far as I know it has been used also for benchmarking by external experts, therefore I presume it is not confidential in any way.

Best regards,
jaj

Hello,

This is the same model which has been applied for a presentation prepared for the last Scientific Committee meeting, although I am not sure it has been actually presented there, there is a number of plots available exactly on the subject and an analyse of a colleague of mine. Maybe you have already this presentation?

Best regards,
jaj

Hello,

I include as the illustration of the symptoms a sketchy plot done after a relatively small number of time steps. This is a realively coarsely resolved slope attacked by the current with a shoreline along it. The water level difference between two parallel runs with different number of processors (variable v20, in meters) is relatively small, but large enough to influence the current deflected by this shore, see differences between red and black vectors. By bedload transport this is enough to cause differences accumulating in time, especially by long-term computations.

Well I am very realistic about the reproducibility of the results between serial and parallel runs with variable number of partitions and I do accept some methodic errors. However I would prefer to have methods in which I can clearly influence the accuracy with some steering variables in order to find some balance between the cost of computations and the degree of correctness which is acceptable to me. I think most of the users would agree.

Best regards,
jaj

Hello,

here is a quick update of further investigations.

(1) With Sisyphe decoupled from Telemac some of the locally occuring differences remain, some disappear, some new appear. Seem to be consistent with occuring erosion/deposition somehow.

(2) With teh TYPE OF ADVECTION =1;5 (changed for velocities to the characteristics)... some of the locally occuring differences remain, some disappear, some new appear. But it seems to be inconsistent and not understandable why here and not there.

(3) After observing a lot of plots with differences between waterlevels in runs with different number of partitions: The free surface is very wavy, and these small instabilities also differ with the number of partitions which is clearly seen. With collocated linear elements:

DISCRETIZATIONS IN SPACE = 11;11

we have always some small instabilities (well, they can be instable due to the theory...), although the users precautionally set implicit depth computations and an artificial filter for instabilities.

IMPLICITATION FOR DEPTH = 1.0
FREE SURFACE GRADIENT COMPATIBILITY = 0.8

The problem is, the amplitude of these waves can be already in the order of minimum depth threshold for wetting or drying... Or overlap on waterlevels causing disturbances seemingly depending on the number of partitions.

I will continue with the investigations, but now they are inconclusive.

Best regars,
jaj

Hello,

the gruesome news from my further computations is that changing the strength of the filter for the free surface:

FREE SURFACE GRADIENT COMPATIBILITY

from 1.0 to 0.8 and then 0.5 and 0.1 brings astonishing(?) results: from a great number of larger erroneous places (1.0) over quite a few in varying places (0.9->0.5) to two-three isolated single nodes on dried-out surfaces (0.1).

This is ugly, because it is connected with the largest method weakness. It gives us a hint that the oscillations due to "natural" instabilities in (collocated linear) finite elements have something to do with it. Or things which trigger these instabilities. Ugly!

Best regards,
Jacek

Hello,

a small update here. It seems that

FREE SURFACE GRADIENT COMPATIBILITY = 0.0

is the hammer for most troubles. We have then inconsistencies in parallel runs only in some awkwardly resolved small channeles. min-bays, where the wetting-and-drying algorithms have their troubles anyway. And away from mesh partition boundaries.

There are two kinds of troubles it seems:

(1) Ones definitely connected with the crossing of the shoreline and an interface between partitions, preferably on very steep shores resolved coarsely. But strangely, they would disappear when the instabilities are damped.

(2) Ones appearing also away from interfaces. They seem to be caused by the instabilities of the free surface and can be influenced by controlling the way the gradients of the free surface are computed. They would disappear by the "piecewise constant" way in contrast to the "linear" way, i.e. setting the mentioned above parameter from its default value of 1.0 drastically down to 0.0.

Best regards,
jaj

Hello Jacek

Did you try with finite volume kernel? the use the HLLC scheme is recommanded.
The wetting and drying algorithm is stricktly different from the finite elements one and thus, if this latter is the main cause of troubles, we can focus on it.

with my best regards
Riadh

Dear Riadh,

has the "finite volume kernel" the same features as the "finite element kernel"? Does it solve the same equations with the same initial and boundary conditions? With the same theory?

I haven't worked with Telemac for quite a while but I connect the "FV" branch with the hyperbolic theory for advection-dominated flows (all that Godunov-, Riemann-, Roe-approaches, etc.). You seem to confirm this naming the HLLC solver, so a pure example of this kind of hyperbolic approach for a specific flow classes.

Unfortunately we deal here with boringly quasi-stationary lowland river flows controlled mainly by friction and only partially in shear zones by turbulent dissipation... All that catastrophic dam breaks are definitaly not our job, save for some locally appearing effects over dykes (obstacles, etc.) on forelands while computing floods, but floods slowly increasing and decreasing... "We only need advection to get our river meanders right..."

So respect, respect for dealing with all that shock waves, but we -- for at least well explained theoretical reasons -- would stick to solvers for "general environmental flows".

Anyway thank you for your message.

Best regards,
jaj