Please comment on the way I discover and fix the bug. Modified fortran code and .cas files are attached.
Program: telemac3d\particles
To trigger the bug,
(1) Change “NON-HYDROSTATIC VERSION” in t3d_particles.cas to “NO”
(2) Change “PRINTOUT PERIOD FOR DROGUES” in t3d_particles.cas to “1”
(3) Change “user_fortran\user_flot3d.f” to:
IF(LT.EQ.10) THEN
CALL ADD_PARTICLE(-220.D0,400.D0+LT/3.D0,259.D0+LT/100.D0,
& LT,NFLOT,
& NFLOT_MAX,XFLOT,YFLOT,ZFLOT,TAGFLO,
& SHPFLO,SHZFLO,ELTFLO,ETAFLO,MESH3D,NPLAN,
& 0.D0,0.D0,0.D0,0.D0,0,0)
ENDIF
Then run.
Expected result: Particle at X=-220.0, Y=403.33, Sigma=0.4423 at LT=10 (t=50 s)
Observed result: Particle at X=-218.79576450, Y=403.36784403, Sigma=0.38305969, Y=, Z= at LT=10 (t=50 s) and Sigma=0.00000 afterwards.
The discrepancy between expected and observed X and Y is understandable (there may be one-time-step advection before printing out), the discrepancy in is not. Further checking shows the size of WCONV=1 while the size of UCONV=20390, VCONV=20390. The number of points for the domain is 2039*10 horizontal layers=20390, WCONV=1 is run into variable W in subroutine DERIVE (v8r1\sources\utils\bief\derive.f) and give problematic result.
Correction: change WCONV to WSCONV (line 2129-2140) in the CALL DERIVE in telemac3d.f.
CALL DERIVE(UCONV%R,VCONV%R,
WSCONV%R,DT,AT,
& X,Y,Z,
& MESH2D%IKLE%I,MESH3D%IFABOR%I,LT,IELM3,UCONV%ELM,
& 3,3,
& NPOIN3,NPOIN2,NELEM2,MESH2D%NELMAX,
& MESH2D%SURDET%R,XFLOT%R,YFLOT%R,ZFLOT%R,
& SHPFLO%R,SHZFLO%R,TAGFLO%I,ELTFLO%I,ETAFLO%I,
& NFLOT,NFLOT_MAX,FLOPRD,MESH3D,T3D_FILES(T3DFLO)%LU,
& IT1%I,T3_01%R,T3_02%R,T3_03%R,IT2%I,
! NO STOCHASTIC DIFFUSION
& MTRA1%X%R,MTRA2%X%R,NPOIN3,0,SVIDE,
& NPLAN,ZCHAR%R,TRANSF)
To make the notation consistent, subroutine SCHAR41_STO (which is called by subroutine DERIVE via subroutine SCARACT) in streamline.f (line 3176-3182 and line 3530-3533) are modified as follows:
Line 3176-3182
DZ(IPLOT) = ((W(I1,IET )*SHP(1,IPLOT)
& + W(I2,IET )*SHP(2,IPLOT)
& + W(I3,IET )*SHP(3,IPLOT))*(1.D0-SHZ(IPLOT))
& +(W(I1,IET+1)*SHP(1,IPLOT)
& + W(I2,IET+1)*SHP(2,IPLOT)
& + W(I3,IET+1)*SHP(3,IPLOT))*SHZ(IPLOT) )
& * PAS
* (ZSTAR(IET+1)-ZSTAR(IET)) / DELTAZ
Line 3530-3533
DZ(IPLOT) = ( W(I1,IET2)*SHP(1,IPLOT)
& + W(I2,IET2)*SHP(2,IPLOT)
& + W(I3,IET2)*SHP(3,IPLOT) ) * PAS2
& * (ZSTAR(IET+1)-ZSTAR(IET)) / DELTAZ
MOC scheme “SCHEME FOR ADVECTION OF VELOCITIES = 1” is used in the original “.cas” file and the scheme seems to have relevance to the subroutine SCHAR41_SIGMA,
“SCHEME FOR ADVECTION OF VELOCITIES = 14” is used here.
Observed result after the modification:
Particle X=-218.79596354, Y=403.36717802, Sigma=0.44551731 at LT=10 (t=50 s) and further transportation seems making sense.