open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi all,

I am trying to implement the Parker (1990) bedload transport equation in GAIA. As it is the basis for the Wilcock and Crowe (2003) equation, I thought I would start from the BEDLOAD_WILCOCK_CROWE_GAIA subroutine instead of the USER_BEDLOAD_QB subroutine.

After a few compilation errors, I found myself with a simulation that runs, but almost immediately creates NaN values over the whole domain for most output variables, such as BOTTOM, VELOCITY, BEDLOAD, etc. Listing printouts are also mostly composed of NaN, and I keep getting "GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS: 50 RELATIVE PRECISION: NaN". Interestingly, it's "almost" immediately as in the early timesteps of one simulation I tried, only part of the BOTTOM variables was NaN values and the rest was as it should be.

A summary of how the subroutine should work, past other elements that were already present for Wilcock and Crowe:

1. The geometric mean diameter D_sg is calculated by summing the product of the ratios and diameters at the given node (grabbed from DECLARATIONS_GAIA). There is a natural log() during the summation and then exp() at the end.
2. The arithmetic standard deviation in the phi scale σ_φ is likewise calculated by summing the product of the ratios and diameters, but including a division by D_sg, hence the need for a second loop. Natural log() and sqrt() are used throughout and at the end.
3. The non-dimensional reference shear stress φ_sg0 is calculated using the bed shear stress and the reference shear stress.
4. The reference straining factor ω₀ and the reference standard deviation in the phi scale σ_φ0 are approximated using a linearized parametric curve. I have checked the logic on this one and at least in Excel, it seems to work.
5. The straining factor ω is calculated using the reference values ω₀ and σ_φ0 calculated above, as well as the arithmetic standard deviation in the phi scale.
6. The "dummy variable" φ is calculated using the straining factor ω, the non-dimensional shear stress φ_sg0, and the ratio between the current diameter D_i and geometric mean diameter D_sg.
7. The function G(φ) is calculated using the 3-part equation in the paper.
8. The bedload QSC%R(I) is calculated using W_si^* which is itself a function of G(φ) and a constant.

I've tried this subroutine with 2 different geometries, and both give the NaN problem, so it should be simulation-independent. One is a cold start and the other is a hydrodynamic hot start. In the steering file, you can specify A,G,L,O (PRIVE_1-4) as output variables to see computed values for D_sg, σ_φ, σ_φ0, and ω₀ (meanings explained below). These values should be independent of the sediment class at study. In testing, you may also want to try outputting φ_sg0 or φ or G(φ) (although these will depend on which sediment class is being iterated over).

So, where's it going wrong? No clue! Honestly, I don't know where a NaN could be introduced, nor do I know how to find what's introducing it. Part of me thinks the EX_ variables could cause issues, but they don't seem to be, given the graphic output PRIVE that seems to work. I'm also concerned that bed shear stress τ seems to go to 0 or NaN very quickly; I'm not sure what could cause this.

If it can help, I'm testing the bedload of particle diameters between 0.0063 m and 0.112 m, in case those values trigger something strange.

I know this is a big ask, but any help in debugging this subroutine is more than welcome! Again, it compiles successfully, but gives NaN errors across the domain quickly but not immediately.
As far as I can tell, it matches the implementation given in Parker (1990) with the exception of the hacky parametric curve for ω₀ and σ_φ0, but if any mistakes on that front are evident, please let me know.

Kind regards,
André Renault

Hello again André,

As told in other posts, when getting NaN in your listing (Not a Number), you should use a debug configuration with debug options to investigate.
See e.g. the S11.gfortran.debug configuration in the $HOMETEL/configs/systel.edf.cfg configuration file, in particular the flag fflags_debug_gfo for gfortran compiler.
fflags_debug_gfo: -g -Wall -fcheck=all -fbacktrace -fbounds-check -finit-integer=-1 -finit-real=nan -ffpe-trap=invalid,zero,overflow

It will show in which subroutine the issue/nan occurs, the first suspicious line and may help you to change something in your computation.

It may be a division by 0, a square root of a negative float, a logarithm of a negative float...

Hope this helps,

Chi-Tuan