From c77f94b2ca40228b2503bc418da5e540dc50e628 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Chi-Tu=C3=A2n=20Pham?= <chi-tuan.pham@edf.fr>
Date: Thu, 26 Dec 2024 17:38:00 +0100
Subject: [PATCH] [VnV][telemac3d] Update pildepon example with better choices
 for a few keywords - TREATMENT OF NEGATIVE DEPTHS = 2 (new default value
 since release 9.0) to improve mass conservation of water (already chosen for
 NH example) - use of default advection scheme for velocities (LIPS): better
 parallelisation than with method of characteristics - use of gradient
 conjugate to solve PPE (little bit more efficient)

Add comparison of the advection of tracers with PC1/PC2/LIPS

Same time step for non-hydrostatic and hydrostatic versions (now = 0.4 s as NH
before)
---
 examples/telemac3d/pildepon/doc/pildepon.tex  | 119 +++++++++++++++++-
 .../pildepon/f3d_piledepon-nonhydro.slf       |   4 +-
 examples/telemac3d/pildepon/f3d_piledepon.slf |   4 +-
 .../pildepon/t3d_piledepon-nonhydro.cas       |  41 +++---
 examples/telemac3d/pildepon/t3d_piledepon.cas |  31 +----
 examples/telemac3d/pildepon/vnv_pildepon.py   |  18 ++-
 6 files changed, 155 insertions(+), 62 deletions(-)

diff --git a/examples/telemac3d/pildepon/doc/pildepon.tex b/examples/telemac3d/pildepon/doc/pildepon.tex
index 5d2e65f3b0..5e21a35dde 100644
--- a/examples/telemac3d/pildepon/doc/pildepon.tex
+++ b/examples/telemac3d/pildepon/doc/pildepon.tex
@@ -6,6 +6,7 @@ This example demonstrates the ability of \telemac{3D} to represent
 the impact of an obstacle on a channel flow.
 It also demonstrates the capability to represent unsteady eddies in a
 model with steady state boundary.
+Finally, it compares the advection of tracers with several methods.
 
 \section{Description}
 
@@ -62,15 +63,44 @@ is made of 4,304 triangular elements (2,280 nodes).
 2 computations are run: one with the hydrostatic hypothesis and one with
 non-hydrostatic version.
 
-To solve the advection, the method of characteristics is used for velocities.
+To solve the advection, the LIPS scheme is used for velocities (default choice).
+It used to be the method of characteristics but there were differences between
+sequential and parallel runs with this method for this example.
+
+GMRES is used for solving the propagation (option 7), conjugate gradient for the
+diffusion of velocities (option 1).
+Conjugate gradient is used to solve the Poisson Pressure Equation because
+a little bit more efficient for this example.
 
-GMRES is used for solving the propagation and diffusion of velocities (option 7).
 The implicitation coefficients for depth and velocities are respectively equal
 to 0.6 and 1.
 
-The time step is 0.1~s for the hydrostatic version, 0.4~s for the
-non-hydrostatic version.
+The time step is 0.4~s for both hydrostatic and non-hydrostatic versions.
 The simulated period is 80~s for both cases.
+\\
+
+Default value since release 9.0 for keyword \telkey{TREATMENT OF NEGATIVE DEPTHS}
+(= 2, i.e. flux control) is chosen to improve the mass-conservation of water
+(up to machine precision).
+\\
+
+In the non-hydrostatic example, 10 tracers are released from the inlet
+(initial value = 0 for every tracer).
+Horizontal coefficient for diffusion of tracers is equal to 0.01~m$^2$/s
+whereas vertical coefficient for diffusion of tracers is equal to 0.1~m$^2$/s.
+Every tracer is advected with a different advection scheme:
+\begin{itemize}
+\item Method of characteristics (= 1),
+\item SUPG,
+\item Leo Postma,
+\item MURD N-type,
+\item MURD PSI-Type,
+\item NERD \#13,
+\item NERD \#14,
+\item Predictor-Correct order 1,
+\item Predictor-Correct order 2,
+\item LIPS.
+\end{itemize}
 
 \subsection{Physical parameters}
 
@@ -128,6 +158,87 @@ computations.
   \label{t3d:pildepon:FreeSurfNH}
 \end{figure}
 
+Mass-conservation of tracer is not ensured with the method of characteristics.
+With SUPG, it is not as good as other distributive schemes.
+In the distributive schemes family, mass conservation of tracer with PSI and
+NERD \#14 is less well preserved compared to the other advection schemes
+(Leo Postma, N, NERD \#13, NERD Leo Postma for Tidal Flats, Predictor-Corrector
+order 1 and 2 + LIPS).
+
+\begin{lstlisting}[language=TelFortran]
+                FINAL MASS BALANCE
+T =          80.0000
+
+--- WATER ---
+INITIAL VOLUME                      :     1637.919
+FINAL VOLUME                        :     1649.809
+VOLUME EXITING (BOUNDARY OR SOURCE) :    -11.88995
+TOTAL VOLUME LOST                   :   -0.5378809E-11
+
+--- TRACER 1: TRACEUR 1       , UNIT : CARACTERISTIQUES* M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     122.0862
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -161.2828
+TOTAL QUANTITY OF TRACER LOST       :     39.19659
+
+--- TRACER 2: TRACEUR 2       , UNIT : SUPG            * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     112.9911
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -112.9831
+TOTAL QUANTITY OF TRACER LOST       :   -0.7956231E-02
+
+--- TRACER 3: TRACEUR 3       , UNIT : LPO             * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     131.6895
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -131.6895
+TOTAL QUANTITY OF TRACER LOST       :    0.3009006E-09
+
+--- TRACER 4: TRACEUR 4       , UNIT : NSC             * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     131.6393
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -131.6393
+TOTAL QUANTITY OF TRACER LOST       :    0.3065850E-09
+
+--- TRACER 5: TRACEUR 5       , UNIT : PSI             * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     124.5492
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -124.5492
+TOTAL QUANTITY OF TRACER LOST       :   -0.1365517E-05
+
+--- TRACER 6: TRACEUR 6       , UNIT : LPO_TF          * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     131.7259
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -131.7259
+TOTAL QUANTITY OF TRACER LOST       :   -0.7301537E-09
+
+--- TRACER 7: TRACEUR 7       , UNIT : NSC_TF          * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     131.6728
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -131.6728
+TOTAL QUANTITY OF TRACER LOST       :    0.8511336E-07
+
+--- TRACER 8: TRACEUR 8       , UNIT : PC1             * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     123.0018
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -123.0018
+TOTAL QUANTITY OF TRACER LOST       :    0.7421335E-09
+
+--- TRACER 9: TRACEUR 9       , UNIT : PC2             * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     122.9811
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -122.9811
+TOTAL QUANTITY OF TRACER LOST       :    0.5133813E-09
+
+--- TRACER10: TRACEUR 10      , UNIT : LIPS            * M3)
+INITIAL QUANTITY OF TRACER          :     0.000000
+FINAL QUANTITY OF TRACER            :     128.0557
+QUANTITY EXITING (BOUNDARY/SOURCE)  :    -128.0557
+TOTAL QUANTITY OF TRACER LOST       :    0.4271499E-09
+\end{lstlisting}
+
+SUPG and LIPS for the advection of tracers are not the best parallelised
+advection schemes for this example, to be investigated.
+
 \section{Conclusion}
 
 \telemac{3D} can be used to study the hydrodynamic impact of engineering
diff --git a/examples/telemac3d/pildepon/f3d_piledepon-nonhydro.slf b/examples/telemac3d/pildepon/f3d_piledepon-nonhydro.slf
index 3690798b2b..023e250049 100644
--- a/examples/telemac3d/pildepon/f3d_piledepon-nonhydro.slf
+++ b/examples/telemac3d/pildepon/f3d_piledepon-nonhydro.slf
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:0ec845fe6fddfde859dc29194b8bbcd54b215de94f3e4796aa8bce082aa7e9da
-size 1885361
+oid sha256:855f5e4e48551af9b80560cab14c5c631050d327f5f2dddec3c8169492b62b23
+size 1447644
diff --git a/examples/telemac3d/pildepon/f3d_piledepon.slf b/examples/telemac3d/pildepon/f3d_piledepon.slf
index fd17f31399..70a1e9a40d 100644
--- a/examples/telemac3d/pildepon/f3d_piledepon.slf
+++ b/examples/telemac3d/pildepon/f3d_piledepon.slf
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:d8486c6413ff1632130cc72d1e5a85d82cc3fd27e4a00c68082fa8331b85c0bb
-size 2651325
+oid sha256:ba48b659117f662ae21e00d991f212b513b96c8ba11416ef0c693be68a9c41f1
+size 899964
diff --git a/examples/telemac3d/pildepon/t3d_piledepon-nonhydro.cas b/examples/telemac3d/pildepon/t3d_piledepon-nonhydro.cas
index de4499b28d..0ddec245b5 100644
--- a/examples/telemac3d/pildepon/t3d_piledepon-nonhydro.cas
+++ b/examples/telemac3d/pildepon/t3d_piledepon-nonhydro.cas
@@ -13,20 +13,15 @@
 /  SUR HPC3700  : version 6.1 : 285 s  21/01/11
 /  SUR HPC3700  : version 6.2 : 205 s  30/05/12  Calibre 7 : 37 s (1 proc)
 /
+TITLE = 'TELEMAC 3D : RIVIERE AVEC DEUX PILES DE PONT$'
+/
 BOUNDARY CONDITIONS FILE : geo_piledepon.cli
 GEOMETRY FILE            : geo_piledepon.slf
 LIQUID BOUNDARIES FILE   : t3d_piledepon.qsl
 3D RESULT FILE           : r3d_piledepon-nonhydro.slf
 2D RESULT FILE           : r2d_piledepon-nonhydro.slf
-/----------------------------------------------------------------------/
-FREE SURFACE GRADIENT COMPATIBILITY=0.75
-MASS-BALANCE : YES
-TIDAL FLATS : YES
-TREATMENT OF NEGATIVE DEPTHS : 2
-TREATMENT OF FLUXES AT THE BOUNDARIES : 2;2
-NUMBER OF HORIZONTAL LEVELS = 6
 /
-TITLE = 'TELEMAC 3D : RIVIERE AVEC DEUX PILES DE PONT$'
+NUMBER OF HORIZONTAL LEVELS = 6
 /
 / NOTE : IL FAUDRAIT 400 PAS DE TEMPS POUR VOIR LA PERIODICITE
 /
@@ -40,6 +35,7 @@ GRAPHIC PRINTOUT PERIOD = 200
 LISTING PRINTOUT PERIOD = 20
 VARIABLES FOR 2D GRAPHIC PRINTOUTS = 'U,V,H,B,S,TA*'
 VARIABLES FOR 3D GRAPHIC PRINTOUTS = 'Z,U,V,W,TA*'
+/
 LAW OF BOTTOM FRICTION = 3
 FRICTION COEFFICIENT FOR THE BOTTOM = 40.
 /
@@ -50,6 +46,7 @@ COEFFICIENT FOR VERTICAL DIFFUSION OF VELOCITIES   = 1.D-6
 MIXING LENGTH MODEL         = 1   / OLD DEFAULT VALUE UNTIL V8P5
 /
 IMPLICITATION FOR DEPTH = 0.6
+FREE SURFACE GRADIENT COMPATIBILITY = 0.75
 INITIAL CONDITIONS : 'ZERO ELEVATION'
 PRESCRIBED FLOWRATES  : 0.;62.
 PRESCRIBED ELEVATIONS : 0.;0.
@@ -58,32 +55,37 @@ PRESCRIBED ELEVATIONS : 0.;0.
 /
 NON-HYDROSTATIC VERSION : YES
 /
+MASS-BALANCE : YES
+/
 MASS-LUMPING FOR DEPTH : 1.
 /  SYSTEMES LINEAIRES DE L'OPTION NON-HYDROSTATIQUE
-ACCURACY FOR PPE = 1.E-5
+SOLVER FOR PPE = 1
 /
-SCHEME FOR ADVECTION OF VELOCITIES : 1
+TREATMENT OF FLUXES AT THE BOUNDARIES : 2;2
 /
-NUMBER OF TRACERS : 7
+NUMBER OF TRACERS : 10
 NAMES OF TRACERS : 'TRACEUR 1       CARACTERISTIQUES';
 'TRACEUR 2       SUPG            ';
 'TRACEUR 3       LPO             ';
 'TRACEUR 4       NSC             ';
 'TRACEUR 5       PSI             ';
 'TRACEUR 6       LPO_TF          ';
-'TRACEUR 7       NSC_TF          '
-INITIAL VALUES OF TRACERS : 0.;0.;0.;0.;0.;0.;0.
+'TRACEUR 7       NSC_TF          ';
+'TRACEUR 8       PC1             ';
+'TRACEUR 9       PC2             ';
+'TRACEUR 10      LIPS            '
+INITIAL VALUES OF TRACERS : 0.;0.;0.;0.;0.;0.;0.;0.;0.;0.
 /
 /VALEURS MISES DANS CLI.TXT
 / GMRES IS CHOSEN FOR SUPG SCHEME (AND DOES NOT WORK WITH PRECONDITIONING 17)
-SCHEME FOR ADVECTION OF TRACERS          : 1;2; 3; 4; 5;13;14
-SOLVER FOR DIFFUSION OF TRACERS          : 1;7; 1; 1; 1; 1; 1
-PRECONDITIONING FOR DIFFUSION OF TRACERS :34;2;34;34;34;34;34
+SCHEME FOR ADVECTION OF TRACERS          : 1;2; 3; 4; 5;13;14; 5; 5; 5
+SCHEME OPTION FOR ADVECTION OF TRACERS   : 1;1; 1; 1; 1; 1; 1; 2; 3; 4
+SOLVER FOR DIFFUSION OF TRACERS          : 1;7; 1; 1; 1; 1; 1; 1; 1; 1
+PRECONDITIONING FOR DIFFUSION OF TRACERS :34;2;34;34;34;34;34;34;34;34
 SUPG OPTION : 1
 ACCURACY FOR DIFFUSION OF TRACERS : 1.D-12
 COEFFICIENT FOR HORIZONTAL DIFFUSION OF TRACERS : 0.01
 COEFFICIENT FOR VERTICAL DIFFUSION OF TRACERS   : 0.1
-ACCURACY FOR PROPAGATION : 1.E-10
 /
 / POUR ESSAI DE L'ELEMENT 51
 /
@@ -99,8 +101,3 @@ ACCURACY FOR PROPAGATION : 1.E-10
 /SCHEMA POUR LA CONVECTION DES VITESSES : 5
 /SCHEMA POUR LA CONVECTION DES TRACEURS : 2;3;4;5;13;14
 /SOLVEUR POUR LA DIFFUSION DES TRACEURS : 7;1;1;1; 1; 1
-/
-/ DEFAULT VALUES UNTIL V8P0 KEPT FOR NON REGRESSION
-IMPLICITATION FOR VELOCITIES = 1.
-SCHEME OPTION FOR ADVECTION OF TRACERS   : 1;1; 1; 1; 1; 1; 1
-OPTION OF SOLVER FOR DIFFUSION OF TRACERS = 3
diff --git a/examples/telemac3d/pildepon/t3d_piledepon.cas b/examples/telemac3d/pildepon/t3d_piledepon.cas
index 3fa9f7b518..01223e101e 100644
--- a/examples/telemac3d/pildepon/t3d_piledepon.cas
+++ b/examples/telemac3d/pildepon/t3d_piledepon.cas
@@ -1,16 +1,3 @@
-/  HP C3700 : 314 s en version 5.3
-/  HP C3700 : 271 s en version 5.4
-/  HP C3700 : 128 s en version 5.5  16/12/2004
-/  HP C3700 : 119 s en version 5.5  04/03/2005 (nouvelle diffusion)
-/  HP C3700 : 113 s en version 5.6  24/08/2005
-/  HP C3700 : 101 s en version 5.6  25/08/2005 (BANDEC=.NO.)
-/  HP C3700 : 107 s en version 5.7  13/02/2007 (Nag:268 s,Dell :61 s)
-/  HP C3700 :  95 s en version 5.8  19/12/2007 (Nag:239 s,Dell :50 s)
-/  HP C3700 :  97 s en version 5.9  08/08/2008 (Nag:267 s,Dell :50 s)
-/  HP C3700 :  94 s en version 6.0  04/12/2009 (Nag:251 s,Intel:38 s)
-/  HP C3700 :  99 s en version 6.1  21/01/2011
-/  HP C3700 : 100 s en version 6.2  30/05/2012
-/
 TITLE = 'TELEMAC 3D : RIVIERE AVEC DEUX PILES DE PONT$'
 /
 BOUNDARY CONDITIONS FILE : geo_piledepon.cli
@@ -21,10 +8,10 @@ LIQUID BOUNDARIES FILE   : t3d_piledepon.qsl
 /
 NUMBER OF HORIZONTAL LEVELS = 6
 /
-TIME STEP = 0.1
-NUMBER OF TIME STEPS = 800
+TIME STEP = 0.4
+NUMBER OF TIME STEPS = 200
 GRAPHIC PRINTOUT PERIOD = 100
-LISTING PRINTOUT PERIOD = 50
+LISTING PRINTOUT PERIOD = 20
 VARIABLES FOR 2D GRAPHIC PRINTOUTS = 'U,V,H,B,S'
 VARIABLES FOR 3D GRAPHIC PRINTOUTS = Z,U,V,W
 /
@@ -33,9 +20,9 @@ FRICTION COEFFICIENT FOR THE BOTTOM = 40.
 /
 COEFFICIENT FOR HORIZONTAL DIFFUSION OF VELOCITIES = 0.005
 HORIZONTAL TURBULENCE MODEL = 1
-VERTICAL TURBULENCE MODEL = 2
+VERTICAL TURBULENCE MODEL   = 2
 COEFFICIENT FOR VERTICAL DIFFUSION OF VELOCITIES   = 1.D-6
-MIXING LENGTH MODEL       = 1   / OLD DEFAULT VALUE UNTIL V8P5
+MIXING LENGTH MODEL         = 1   / OLD DEFAULT VALUE UNTIL V8P5
 /
 IMPLICITATION FOR DEPTH = 0.6
 FREE SURFACE GRADIENT COMPATIBILITY = 0.75
@@ -46,11 +33,3 @@ PRESCRIBED ELEVATIONS : 0.;0.
 NON-HYDROSTATIC VERSION : NO
 /
 MASS-BALANCE : YES
-/
-/ DEFAULT VALUES UNTIL V8P0 KEPT FOR NON REGRESSION
-SCHEME FOR ADVECTION OF VELOCITIES = 1
-IMPLICITATION FOR VELOCITIES = 1.
-ACCURACY FOR PROPAGATION = 1.E-6
-/
-/ DEFAULT VALUE UNTIL V8P5 KEPT FOR NON REGRESSION
-TREATMENT OF NEGATIVE DEPTHS = 1
diff --git a/examples/telemac3d/pildepon/vnv_pildepon.py b/examples/telemac3d/pildepon/vnv_pildepon.py
index 622063a88c..6ea6834562 100644
--- a/examples/telemac3d/pildepon/vnv_pildepon.py
+++ b/examples/telemac3d/pildepon/vnv_pildepon.py
@@ -77,32 +77,38 @@ class VnvStudy(AbstractVnvStudy):
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_1:T3DRES',
                             'f3d_piledepon.slf',
-                            eps=[1.E-1, 2.E0, 1.5E0, 1.E0])
+                            eps=[1.e-9])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_2:T3DRES',
                             'f3d_piledepon.slf',
-                            eps=[1.E-1, 2.E0, 1.5E0, 1.1E0])
+                            eps=[1.e-9])
 
         # Comparison between sequential and parallel run.
         self.check_epsilons('vnv_1:T3DRES',
                             'vnv_2:T3DRES',
-                            eps=[1.E-1, 2.E0, 1.5E0, 1.1E0])
+                            eps=[1.e-9])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_3:T3DRES',
                             'f3d_piledepon-nonhydro.slf',
-                            eps=[0.04, 1.5, 1., 0.4, 0.2, 0.3, 0.1, 0.2, 0.2, 0.1, 0.1])
+                            eps=[1.e-9, 1.e-9, 1.e-9, 1.e-9, 3.e-4, 1.e-10,\
+                                 1.e-10, 1.e-10, 1.e-10, 1.e-10, 1.e-10,\
+                                 1.e-10, 1.e-10, 1.e-10])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_4:T3DRES',
                             'f3d_piledepon-nonhydro.slf',
-                            eps=[0.2, 1.7, 1.1, 0.4, 0.3, 0.6, 0.3, 0.3, 0.4, 0.3, 0.3])
+                            eps=[1.e-9, 1.e-9, 1.e-9, 1.e-9, 6.e-4, 1.e-10,\
+                                 1.e-10, 1.e-10, 1.e-10, 1.e-10, 1.e-10,\
+                                 1.e-10, 1.e-10, 3.e-4])
 
         # Comparison between sequential and parallel run.
         self.check_epsilons('vnv_3:T3DRES',
                             'vnv_4:T3DRES',
-                            eps=[0.2, 1.7, 1.1, 0.4, 0.3, 0.6, 0.3, 0.3, 0.4, 0.3, 0.3])
+                            eps=[1.e-9, 1.e-9, 1.e-9, 1.e-9, 6.e-4, 1.e-10,\
+                                 1.e-10, 1.e-10, 1.e-10, 1.e-10, 1.e-10,\
+                                 1.e-10, 1.e-10, 3.e-4])
 
 
     def _post(self):
-- 
GitLab