open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear all,

I know my question is a bit too philosophic, but maybe useful for othe forum members. I must admit I have always felt confused in the depths of subroutine propag.f of Telemac-2D ( ehm, I am three-dimensional from birth ) therefore I decided to ask anyone who knows it well.

There are to pricipal algorithms for treatment of the linear system and I would like to concentrate on the one based on the wave equation:

TREATMENT OF THE LINEAR SYSTEM = 2 // SOLSYS

Compared to the well known and widely appreciated algorithm of Casulli (also based on the wave equation), where the horizontal advection and diffusion operators of shallow water equatons are treated explicitely in order to simplify drastically the resulting linear equation system to a diagonally dominated one, which can be extremely efficiently solved iteratively -- in Telemac-2D, so far as I understand from propag.f, at least the diffusion operator is taken semi-implicitely, making the matrix to be solved more awkward. Even more, the default value of the implicitness factor (which most users do not touch) is

IMPLICITATION FOR DIFFUSION OF VELOCITY = 1. // TETAD

i.e. it is suggested to take velocity diffusion fully implicitely (i.e. in theory no stability constraints on the time step but the physics, however with maybe more numerical "diffusion of diffusion").

My question is, what are your experiences with semi-implicit or even explicit velocity diffusion treatment in the wave-equation-treatment branch? Is it applicable? And by the way, does this branch assume explicit advection?

Best regards,
jaj

Hello Jacek,

The answer is rather technical. The diffusion of the velocity is basically explicit, even if the user asks for implicitation, but we can keep at least the diagonal part of the matrix implicit, and this is more stable. The point is that if you keep only diagonal terms implicit (i.e. friction and diagonal of diffusion) you can still solve explicitely the momentum equation, to get (just by dividing by a diagonal) a predictor value of the velocity at time n+1, that can be plugged into the Saint-Venant continuity equation. In 3D it can be even more tricky and it leads to a linear system solved on every vertical, because we sum over the vertical a number of implicit (but diagonnal only) momentum equations.

I hope this clarifies a bit the matter, otherwise it is described in my book, and in more detail in the first release notes available on this website (a chapter called "revisiting the wave equation in 3D" or something like this, I am not at work...).

With best regards,

Jean-Michel

Dear Jean-Michel,

well, "the diffusion of the velocity is basically explicit, even if the user asks for implicitation" is a bit confusing to me, what is then the use of the implicitness factor TETAD? I presume -- with a difference to the original Casulli algorithm -- that you take the diffusion "pre-treated" in the velocities taken from the momentum equation to be introduced int the continuity equation (resulting with the wave equation).

The important question then is stability. With purely a explicit diffusion operator (treated separately) one has a "mild" stability limitation on the time step applied, something like timestep =< 1 / (factor * nu_max / resolution_min^2), whereby the factor is 4, if I can quote from Casulli's UnTRIM.

Well the users would tend to use large time steps with the wave equation treatment, even with fine meshes, and in 2D they would apply large viscosities due to the 2D averaging (momentum dispersion)...

What is with the stability then? Any estimates on the time step?

Best regards,
jaj

Hello Jacek,

I never investigated stability of the diffusion equation, actually I have seen only once in a user application a divergence due to an explicit diffusion, so the criterion timestep =< 1 / (factor * nu_max / resolution_min^2) is not really a problem, except in Telemac-3D on the vertical with tidal flats, where the resolution tends to 0...

Regards,

JMH

Dear Jean-Michel,

the fact you've never investigated the stability if the diffusion makes me slightly worry, but with fingers crossed... Anyway, let us remain by 2D, the sigma-mesh troubles by nearby the shoreline are well known, but let us treat them as irrelevant here.

Well, the problem is that people do use notoriously Elder dispersion or diffusion model for the momentum as well, and with the original constants untouched, and remember that they were found out empirically for tracer transport...! In this way, quite high "turbulent viscosity" coefficients result. And the finer the mesh, the more chances of violating the stability criterion. The values are much too high than "we like it" physically, due to the fact we have to deal with the 2D momentum dispersion being the artefact of the 2D-integration alone.

Please note that if the calibration requirement is to have constant dissipation parametrisation for a range of waterlevels -- as required for "prognose ability" in this range, and if one is forced to use small/constant depth-independent roughness coefficients due to morphology requirements, larger diffusion/dispersion coefficients are applied in the practice... I know it is ugly -- it makes no chance for horizontal turbulence models to affect anything -- but this is the dirty daily practice!

Well, people blame usually for the (in)famous Telemac-2D instability, besides the collocated linear elements, of course, the thresholds for this and that. Like on the tidal flats or by the shoreline: thresholds appearing not only for water depth, but also for roughness and/or diffusion!...

I was just wondering if we might have also a trouble with the diffusion stability as a trigger with all that instabilities people are fighting with.

Best regards,
jaj

Hello Jacek,

You cannot escape a threshold in 3D for diffusion, since the vertical gradients can become infinite with elements that have no height, thus with DZ=0 in derivative DF/DZ.

For the stability of explicit diffusion, I said in my previous post that the diagonal is kept implicit. If your case is a steady or quasi steady state, which is about always the case in river applications, you can consider trying an implicitation coefficient greater than 1, even equal to 2, it has a stabilising effect, and it should not change the result in steady state.

With best regards,

Jean-Michel

Dear all,

by an attempt to implement the semi-implicit treatment of the horizontal diffusion terms in the way described above in another code it has been proven, that in terms of the linear FD/FV method this seemingly attractive approach in which in the implicit part just only the diagonal terms treated implicitly and the extradiagonal ones explicitly, is unfortunately inconsistent although always stable.

Inconsistent means that the numerical scheme in reality approximates another equation as it is originally supposed to. Inconsistency in one term (diffusion) may affect other correctly approximated terms as well (advection, pressure gradients). For users it is very difficult to grasp inconsistency in a stable scheme if they do not make extensive sensitivity studies with the time and space resolution and/or verification examples for principal numerical properties are missing.

Because of the well known analogies between finite differences, volumes and elements method this fact cannot be ignored for Telemac as well, we have discussed the matter with Jean Michel anew. Now one has a recommendation, that one should use explicit treatment of the diffusion terms by applying the (pseudo) wave quation approach to solve the shallow water equations to be on the safe side of life. Please note that per default, the (suspected to be inconsistent) fully implicit treatment is set in Telemac-2D.

However, please be warned again, that for an explicit scheme a stability constraint for the time step is valid, see the posts before, in the form timestep =< 1 / (factor * nu_max / resolution_min^2).

For more information, please see the descriptions of the wave equation method in the book of J.-M. Hervouet "Hydrodynamics of free sufrace flows", paragraphs 4.12, 5.4 and especially 5.4.2 "With implicit diffusion?".

For even more information on the original (FD/FV) method see papers of Casulli et al. (numerous publications, for open access see e.g. http://dx.doi.org/10.1016/S0895-7177(02)00264-9 and the references there) with the characteristic implicit treatment of the vertical and explicit treatment of the horizontal diffusion/viscosity terms.

Best regards,
jaj

Dear all,

I have made an effort of checking the Telemac validation suite examples with the question if the mentioned above inconsistency in the formulation of the viscous terms by the "pseudo wave equation" treatment of the equations could be detected in the past. Please remember this is a very perverse sort of errors: everything is stable, processes take place, stationary solutions are not affected -- only the progress of the processes is different and this depending on the resolutions in time and space. Without a suspicion and additional efforts it is very difficult to grasp.

The most perverse hypothetical effect would be affecting not only the viscosity, but also the advective velocity and the wave speed according to the (local, global) resolution -- with e.g. artificial wave dispersion.

So I have searched through all the examples for v6p3 where
TREATMENT OF THE LINEAR SYSTEM : 2.

(1) Unfortunately by most (academical) examples where the effect could be quantitatively discovered, the "velocity diffusivity" was switched off ot the viscosity coefficient set to 0. No way to see anything.

(2) By all examples which concern the flows we are interested, stationary flows are studied -- with large viscosity coefficients -- where it is very difficult to grasp this kind of inconsistency.

(3) By examples in which the inconsistency would show up with the largest possibility -- waves, tides, estuaries, unfortunately no convergence studies were done. Additionally in the free surface filtering with FREE SURFACE GRADIENT COMPATIBILITY was applied (and values of even 0.5, 0.0) which would suppress also the above mentioned artifical wave dispersion.

(4) And in all other suitable cases: no convergency investigations.

In my conclusions: the error could go undetected for years. Additionally: there's no example applying IMPLICITATION FOR DIFFUSION OF VELOCITY = 0. (default 1.) which would be in theory free of inconsistency.

Question: What should we do now? Has anyone an idea of a good test example to test the inconsistency existence?

Best regards,
jaj

PS. I add the list of the validation examples I have looked at.

---

.../telemac/v6p3r2/examples

---

Defaults are

TREATMENT OF THE LINEAR SYSTEM = 1
DIFFUSION OF VELOCITY = YES
IMPLICITATION FOR DIFFUSION OF VELOCITY = 1.
VELOCITY DIFFUSIVITY = 1.E-6
TURBULENCE MODEL = 1
FREE SURFACE GRADIENT COMPATIBILITY = 1.

---

Bump case 3

./hydraulic_jump/t2d_hydraulic_jump_v1p0.cas
TREATMENT OF THE LINEAR SYSTEM : 2 SOLVER : 1

TURBULENCE MODEL = 1
VELOCITY DIFFUSIVITY = 0.

---

Solitary wave / conical island

./conical_island/t2d_island.cas
TREATMENT OF THE LINEAR SYSTEM = 2

TURBULENCE MODEL = 1
DIFFUSION OF VELOCITY = NO

---

Breakwater

./break/t2d_break.cas
TREATMENT OF THE LINEAR SYSTEM : 2

TURBULENCE MODEL : 1
VELOCITY DIFFUSIVITY = 0.00001
FREE SURFACE GRADIENT COMPATIBILITY = 0.9

---

River confluence

./confluence/t2d_confluence.cas
TREATMENT OF THE LINEAR SYSTEM = 2

TURBULENCE MODEL = 1
VELOCITY DIFFUSIVITY = 1.E-3

---

Oil spill in a flume

./riv_art/t2d_riv_art.cas
TREATMENT OF THE LINEAR SYSTEM : 2

TURBULENCE MODEL : 3
VELOCITY DIFFUSIVITY : 0.000001

---

"The Wesel river"

./wesel/t2d_wesel.cas
TREATMENT OF THE LINEAR SYSTEM:2 SOLVER:8

TURBULENCE MODEL = 1 VELOCITY DIFFUSIVITY = 2.
/TURBULENCE MODEL = 2

./wesel/t2d_wesel_pos.cas
TREATMENT OF THE LINEAR SYSTEM:2 SOLVER:2 PRECONDITIONING:11
TREATMENT OF THE LINEAR SYSTEM:2 SOLVER:8

TURBULENCE MODEL = 1 VELOCITY DIFFUSIVITY = 2.
/TURBULENCE MODEL = 2

---

Okada free surface tsunami displacement

./okada/t2d_okada.cas
TREATMENT OF THE LINEAR SYSTEM : 2

TURBULENCE MODEL : 1
VELOCITY DIFFUSIVITY = 1.0
FREE SURFACE GRADIENT COMPATIBILITY = 0.9

---

Culm

./culm/t2d_culm.cas
TREATMENT OF THE LINEAR SYSTEM : 2 SOLVER : 1 PRECONDITIONING : 2
VELOCITY DIFFUSIVITY = 2. TURBULENCE MODEL : 1

---

Solitary wave? Trangular shelf

./triangular_shelf/t2d_triangular_shelf.cas
TREATMENT OF THE LINEAR SYSTEM = 2

FREE SURFACE GRADIENT COMPATIBILITY = 0.9
TURBULENCE MODEL = 1
DIFFUSION OF VELOCITY = NO

---

Shoal?

./shoal/t2d_shoal.cas
TREATMENT OF THE LINEAR SYSTEM : 2

FREE SURFACE GRADIENT COMPATIBILITY = 0.9
VELOCITY DIFFUSIVITY = 0.001
TURBULENCE MODEL : 1

---

Gaussian open boubdarries / Thompson non-reflecting BCs

./thomson/t2d_thomson.cas
TREATMENT OF THE LINEAR SYSTEM = 2

TURBULENCE MODEL = 1
VELOCITY DIFFUSIVITY = 0.

---

Dam break over dry bed

./dambreak/t2d_dambreak_v2p0.cas
TREATMENT OF THE LINEAR SYSTEM : 2

TURBULENCE MODEL = 1
VELOCITY DIFFUSIVITY = 0.

---

Tide

./tide/t2d_tide-jmj_real_gen.cas
./tide/t2d_tide-jmj_real_gen.cas
./tide/t2d_tide-jmj_type.cas

TREATMENT OF THE LINEAR SYSTEM : 2

/TURBULENCE MODEL : 3
/DIFFUSION OF VELOCITY : YES
/VELOCITY DIFFUSIVITY : 1.E-6
FREE SURFACE GRADIENT COMPATIBILITY : 0.5

---

Gironde Estuary / Oil spill?

./estu_gir/t2d_estu_gir.cas
TREATMENT OF THE LINEAR SYSTEM : 2

TURBULENCE MODEL : 1
VELOCITY DIFFUSIVITY : 1.D0
FREE SURFACE GRADIENT COMPATIBILITY : 0.

---

Hello Jacek & Telemac Users,

I don't have a perfect idea but some hint.

I did some studies on the dependency of the scheme on the chosen time step size (see my TUC paper from 2014). When I remember correctly, there was a strong dependency on the time step size and steady state solutions.

If we think about implicit schemes and steady state solutions, my experience form other fully implicit schemes is, that there should be no dependency on the time step size when all terms are discretized correctly. (Yes, there are always some restrictions but they are often far away from the CFL condition).

So what I suggest is to use a simple test case (with a steady state solution) to investigate the dependency of the solution on the time step size. You can increase the complexity of the test case step by step.

Dear Lufia,

thinking positively:

(1) We have always the option of taking the viscosity explicitly -- setting IMPLICITATION FOR DIFFUSION OF VELOCITY = 0., but then being confronted with some "mild" limitation on the time step. (The trouble is that in twodimensional river applications people have to take the dispersion due to the vertical velocity averaging into account via "velocity diffusion", what results in huge viscosity coefficients, think of "Elder turbulence model"...)

(2) And, of course, there's the "good old" approach of TREATMENT OF THE LINEAR SYSTEM : 1, which should be free of inconsistencies, or at least from -this- inconsistency...

According to the great philosophical question if steady states reached with a consistent and an inconsistent scheme (for, let's say shallow water equations) will be exactly the same I must admit my ignorance.

Anyway, if you take the questioned diffusion equation only (no advection, no pressure gradients, no friction nor Coriolis) the same steady state will be reached in -this- particular case with both. However, in the inconsistent scheme, with the "numerical" viscosity coefficient depending on the time step and local spatial resolution, this steady state will be reached in the transient way (so as a development in time) in a different way than with the consistent scheme.

I might imagine troubles -- or I would feel insecure with -- when we add the mentioned above additional terms and the inconsistency spreading over these terms (which depends on the dicretisation), but let me not think about it, my mind is already set for the weekend...

Best regards,
jaj