open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi all,
I am making a 3D simulation of Bodensee (Lake Constance) with a constant wind field above, and with thermal stratification as an initial condition. I use K-eps as vertical and horizontal turbulence model. My problem is that after several iterrations (40 to 60) the step "CALCULATION OF VERTICAL VELOCITY STEP" takes more and more time until it hangs up. This happens both for short (30sec) or long (300sec) timesteps. Does anyone have an idea how can I identify what causes the problem and how to solve it?

Thanks in advance
Marcin Kawka

Hello,

If you have such a message it means that you are using the hydrostatic option, you should perhaps try the non hydrostatic option. In case of wind a critical parameter is the THRESHOLD DEPTH FOR WIND which is by default 1. Try also a case without k-epsilon, with a mixing length model for example. Depending on your mesh size, it is not sure that k-epsilon can see the real velocity gradients that would allow a correct behaviour of this model.

You can get inspiration in test case called Berre which is also a lake.

With best regards,

Jean-Michel Hervouet

Thank you for your answer Jean-Michel. I have increased the THRESHOLD DEPTH FOR WIND to higher values, and tried NON-Hydrostatic option and it made my model hang-up after 100-150 iterations. After increasing K-EPS ACCURACY to 1.E-9 model runs more than 300 iterations, however it still hangs up after 310-320. With constant viscosity it works fine, however I would like to avoid using it. Do you think that refinement of the grid can help in my case?
Where can I find it? Looks like it's not included in validation/training packages, or it has a different name?

Thank you in advance,
Marcin Kawka

Hello,

Yes, this test case was not included, though it is a perfect example of a very complex case, including wind, salinity and temperature. I copy below the steering file (but in French...). In this case some elements are totally crushed because we impose that some planes are really horizontal, to follow stratifications. This is the reason why we treat it like a tidal flats case, though there are none. An important keyword could be :

HYDROSTATIC INCONSISTENCY FILTER : YES

to cancel diffusion in places where elements are too much distorted.

Try also if not already the case :

MASS-LUMPING FOR DEPTH : 1.
MASS-LUMPING FOR DIFFUSION : 1.
FREE SURFACE GRADIENT COMPATIBILITY : 0.9

In this Berre lake case we find k-epsilon useful for taking wind into account. Note that in subroutine soukep.f, where most of the k-epsilon is, you have two options, option 2 being hardcoded. You can try also option 1 (simplest k-epsilon implementation) though you will see that in option 2, the wind stress is put directly into the turbulence production, which avoids refining too much at the free surface.

I hope this helps,

Jean-Michel Hervouet






PROCESSEURS PARALLELES : 8
REGIME DE TURBULENCE POUR LE FOND : 2
LOI DE FROTTEMENT SUR LE FOND : 5
COEFFICIENT DE FROTTEMENT POUR LE FOND : 0.001
/LOI DE FROTTEMENT SUR LE FOND : 3
/COEFFICIENT DE FROTTEMENT POUR LE FOND : 35.
/
/DEBUGGER : 1
/**********************************************************************/
/ FICHIER DE DECLARATION DES MOTS CLES DU CODE /
/ TELEMAC3D /
/**********************************************************************/
/
/ POINT SA3 : X=24363 Y=28336
/ CORRESPOND A LAMBERT II ETENDU : 824362,660 E
/ 1828335,516 N
/
/ AJOUTE PAR JMH : PLANS ECRASES +
/
PRESSION DYNAMIQUE DANS L'EQUATION D'ONDE : NON
VITESSE PROJETEE SUR LES PAROIS LATERALES SOLIDES : OUI
VITESSE PROJETEE SUR LE FOND : OUI

CONTOURNEMENT DES VOLUMES NULS : OUI
BANCS DECOUVRANTS : OUI
TRAITEMENT DES HAUTEURS NEGATIVES : 2
MASS-LUMPING POUR LA HAUTEUR : 1.
PRECISION POUR LA DIFFUSION DES TRACEURS : 1.E-10
SOLVEUR POUR LA PROPAGATION : 2
PRECISION POUR LA PROPAGATION : 1.E-6
COMPATIBILITE DU GRADIENT DE SURFACE LIBRE = 0.9
NOMBRE DE SOUS ITERATIONS POUR LES NON LINEARITES : 1
/
/
/
/ Temps sur HP Linux Z600 comp. Intel Version 6.0 1 procs : 158428 s
/ optimisation mt14pp : 1 procs : 143614 s
/ 2 procs : 90339 s
/ 4 procs : 59675 s
/ optimisation mt14pp : 8 procs : 33076 s
/ ancien mt14pp + bypass : 8 procs : 26319 s
/ schema 5 8 procs : 50465 s
/ nouveau mt14pp + bypass : 8 procs : 22740-23800 s
/ caracteristiques + nouveau mt14pp + bypass : 8 procs : 26737 s !!
/
/ avec k-epsilon et nouveaux schemas : 8 procs : 34287 s (10/01/2011)
/ avec k-epsilon et nouveaux schemas : 8 procs : 29140 s (20/01/2011)
/ avec k-epsilon et nouveaux schemas : 8 procs : 33013 s (03/02/2011) mis en ref
/
/
/
/ essai sur 20 jours avec schema 5 : 48633 s (regime de turbulence sur le fond : 1)
/ essai sur 20 jours avec schema 5 : 40800 s (regime de turbulence sur le fond : 2)
/ essai sur 20 jours avec schema 14 : 22698 s (regime de turbulence sur le fond : 1)
/ essai sur 20 jours avec schema 5 : 35649 s (regime de turbulence sur le fond : 2)
/ et Nikuradse 0.2
/ essai sur 20 jours avec schema 5 : 43890 s (regime de turbulence sur le fond : 2)
/ et Nikuradse 0.01 + k-epsilon
/ essai sur 20 jours avec schema 5 : 24001 s (regime de turbulence sur le fond : 2)
/ et Nikuradse 0.01 + Longueur de melange
/
/ essai sur 20 jours avec schema 5 : 22105 s (regime de turbulence sur le fond : 2)
/ et Nikuradse 0.01 + k-epsilon
/ + DZ_MIN = 0.01
/
/ Le 06/12/2011
/ essai sur 20 jours avec schema 5 : 19109 s (regime de turbulence sur le fond : 2)
/ et Nikuradse 0.001 + k-epsilon
/ + DZ_MIN = 0.01
/ Idem le 03/01/2012 : 13199 s !!!!!
/
TITRE = 'ETANG DE BERRE'
/
/

---

/
/ OPTIONS GENERALES /
/

---

/
/
FICHIER FORTRAN : princi.f
FICHIER DES CONDITIONS AUX LIMITES : cl_B
FICHIER DE GEOMETRIE : geo_B
FICHIER DE REFERENCE : ref3d_14_nikuprism VALIDATION : OUI
/
FICHIER DES RESULTATS 3D : res3d_5_nikuprism
FICHIER DES RESULTATS 2D : res2d_5_nikuprism
/
FICHIER DE DONNEES FORMATE 1 : vent_010906.txt
FICHIER DE DONNEES FORMATE 2 : thermique_010906.txt
FICHIER DES FRONTIERES LIQUIDES : qsl_010906.txt
/
/

---

/
/ OPTIONS DE CALCUL dec 2007 /
/

---

/
/
VARIABLES POUR LES SORTIES GRAPHIQUES 2D : U,V,B,S,H
VARIABLES POUR LES SORTIES GRAPHIQUES 3D : Z,U,V,W,TA1,TA2,NUZ
/1er septembre 2006 au 28 fevrier 2007 = 181 jours
NOMBRE DE PAS DE TEMPS : 86400 /30 jours (mois de septembre)
NOMBRE DE PAS DE TEMPS : 57600 /20 jours (mois de septembre)
PAS DE TEMPS : 30.
PERIODE POUR LES SORTIES LISTING : 360 /une sortie toutes les 3 heures
PERIODE POUR LES SORTIES LISTING : 240 /une sortie toutes les 2 heures

PERIODE POUR LES SORTIES GRAPHIQUES : 360 /Une sortie toutes les 3 heures
PERIODE POUR LES SORTIES GRAPHIQUES : 180 /Une sortie toutes les 1,5 heure
/PERIODE POUR LES SORTIES LISTING : 1
/ POUR FABRIQUER LA REFERENCE
/PERIODE POUR LES SORTIES GRAPHIQUES : 57600

NOMBRE DE LISSAGES DU FOND : 0
BILAN DE MASSE : OUI
INFORMATION SUR LE BILAN DE MASSE A CHAQUE SORTIE LISTING : OUI
/
NOMBRE DE PLANS HORIZONTAUX : 31
TRANSFORMATION DU MAILLAGE : 3
/
NOMBRE DE TRACEURS : 2
NOMS DES TRACEURS : 'TEMPERATURE °C ';
'SALINITE G/L '
LOI DE DENSITE : 3
/
VERSION NON-HYDROSTATIQUE = OUI
/
CONDITIONS INITIALES : 'COTE CONSTANTE'
COTE INITIALE : 0.42
COTES IMPOSEES : 0.;0.;0.;0.
DEBITS IMPOSES : 0.;0.;0.;0.
TRAITEMENT DES FLUX AUX FRONTIERES : 2;2;2;2
VALEURS INITIALES DES TRACEURS : 20.1 ; 28.
VALEURS IMPOSEES DES TRACEURS : 14.4; 38.;20.;0.;20.;0.;20.;0.
/
/

---

/
/ TERMES SOURCES /
/

---

/

LOI DE FROTTEMENT SUR LES PAROIS LATERALES : 0

SEDIMENT : NON
VENT : OUI
/
FILTRE LES INCONSISTANCES HYDROSTATIQUES : OUI
PLUIE OU EVAPORATION : NON
/
/

---

/
/ ETAPE DE CONVECTION /
/

---

/
/
SCHEMA POUR LA CONVECTION DES VITESSES : 5
SCHEMA POUR LA CONVECTION DES TRACEURS : 5
SCHEMA POUR LA CONVECTION DU K-EPSILON : 5
/
/

---

/
/ ETAPE DE DIFFUSION /
/

---

/
/
SCHEMA POUR LA DIFFUSION DES TRACEURS : 1
MODELE DE TURBULENCE VERTICAL : 3
MODELE DE TURBULENCE HORIZONTAL : 3
/
PRECONDITIONNEMENT POUR LA DIFFUSION DU K-EPSILON : 34
/
COEFFICIENT DE DIFFUSION HORIZONTAL DES VITESSES : 1.E-6
COEFFICIENT DE DIFFUSION VERTICAL DES VITESSES : 1.E-6
COEFFICIENT DE DIFFUSION HORIZONTAL DES TRACEURS : 1.E-6
COEFFICIENT DE DIFFUSION VERTICAL DES TRACEURS : 1.E-6
/
PRECISION POUR LA DIFFUSION DES VITESSES : 1.E-6
PRECONDITIONNEMENT POUR LA DIFFUSION DES VITESSES : 34
PRECISION POUR PPE : 1.E-6
/
/ LONGUEUR DE MELANGE
/
/MODELE DE TURBULENCE VERTICAL : 2
/MODELE DE LONGUEUR DE MELANGE : 3
/FONCTION D'AMORTISSEMENT : 3
/MODELE DE TURBULENCE HORIZONTAL : 1
/COEFFICIENT DE DIFFUSION HORIZONTAL DES VITESSES : 0.01
/COEFFICIENT DE DIFFUSION HORIZONTAL DES TRACEURS : 0.01
/
/

---

/
/ PROPAGATION /
/

---

/
/
ORDRE DU TIR INITIAL POUR LA HAUTEUR : 1
IMPLICITATION POUR LA HAUTEUR : 1.
IMPLICITATION POUR LES VITESSES : 1.
/
IMPLICITATION POUR LA DIFFUSION : 1.
PRECONDITIONNEMENT POUR LA DIFFUSION DES TRACEURS : 34
MASS-LUMPING POUR LES VITESSES : 1.
MASS-LUMPING POUR LA DIFFUSION : 1.
PRECONDITIONNEMENT POUR LA PROPAGATION : 2
PRECONDITIONNEMENT POUR PPE : 34
STOCKAGE DES MATRICES : 3
/
PROFILS DE VITESSE : 1;1;1;3
/
/&ETA
&FIN

Thank you Jean-Michel, that helped,
Now my model is ruining. Friction velocity caused by wind seem to be correct, but I'm getting too much turbulent kinetic energy in lake waters and this causes total mixing and termocline destruction. I will try manual fixed layer positioning.

Regards,
Marcin Kawka

Hello,

It looks much like our Berre lake study. However in our case we have much less turbulence with k-epsilon than with mixing length. If you used option 2 of k-epsilon in subroutine soukep.f, make sure that the wind formula that is hardcoded here corresponds to what you have in subroutine bord3d. You can also try option 1, then the turbulence due to wind will be only deduced from velocity gradients, which may give smoother results. So now you are on the cutting edge!!

Another thing to check is the average slope of the free surface, compared with measurements, to validate the friction coefficient of wind.

With best regards,

Jean-Michel Hervouet