Reproducible NaNs in TELEMAC-3D Salinity modelisation

System Modules TELEMAC-3D
1

Hello,

I am investigating recurrent [name removed] in a TELEMAC-3D v9.1.0 simulation of a long, shallow, inverse estuary. The equivalent TELEMAC-2D simulations completed successfully, while the 3D model becomes unstable after a variable number of simulated days.

The model uses:

hydrostatic TELEMAC-3D;

5 horizontal levels;

salinity and temperature as tracers;

density law 3;

WAQTEL thermal exchange;

tidal flats;

FES tidal forcing at the ocean boundary;

river inflows, wind, atmospheric pressure, rainfall and evaporation;

a 5 s time step;

MPI on 128 cores.

The most robust configuration tested so far includes:

SCHEME FOR ADVECTION OF VELOCITIES : 1

SCHEME FOR ADVECTION OF TRACERS : 13;13

SCHEME OPTION FOR ADVECTION OF TRACERS : 1;1

SCHEME FOR DIFFUSION OF TRACERS : 1;1

TIDAL FLATS : YES

OPTION FOR THE TREATMENT OF TIDAL FLATS : 1

TREATMENT OF NEGATIVE DEPTHS : 2

BYPASS VOID VOLUMES : YES

THRESHOLD HEIGHT BEFORE CRUSHED ELEMENTS : 0.01

THRESHOLD FOR VISCOSITY CORRECTION ON TIDAL FLATS : 0.2

OPTION FOR LIQUID BOUNDARIES : 2;2;2;2;2;2;2;2;2

TREATMENT OF FLUXES AT THE BOUNDARIES : 2;2;2;2;2;2;2;2;2

NON-HYDROSTATIC VERSION : NO

IMPLICITATION FOR VELOCITIES : 1.

SOLVER FOR PROPAGATION : 7

OPTION OF SOLVER FOR PROPAGATION : 5

ACCURACY FOR PROPAGATION : 1.E-8

In most failed runs, the first visible error is reported during tracer diffusion:

GRACJG EXCEEDING MAXIMUM ITERATIONS

RELATIVE PRECISION: NaN

The solver normally requires only a few iterations before the failure. Therefore, increasing the maximum number of iterations does not seem relevant: the matrix, right-hand side or tracer field appears to contain a non-finite value already.

2

Immediately afterwards:

salinity and temperature become NaN over the domain;

velocity and free surface are subsequently contaminated;

TELEMAC may continue running and writing invalid results instead of stopping at the first NaN.

The instability is reproducible at the same simulated time on different compute nodes. In the best run, all restart records are valid through day 129, while [name removed] appear between days 129 and 130. The first contaminated 2D field contains 73 NaN nodes centred approximately near the ocean boundary.

This location is on or close to the ocean boundary and is similar to the location found in earlier failures. However, the ocean boundary is not located on a drying flat: its bottom elevations range approximately from -16.75 to [name removed] m, and no obviously degenerate adjacent elements were identified.

casa_3d_salinite_2024_new_config1.cas (6.34 KB)

3

Sensitivity tests produced the following results:

Tracer advection scheme 13 is more robust than scheme 14 in the full configuration.

TREATMENT OF FLUXES AT THE BOUNDARIES : 2;… delays the instability and appears stabilising.

Using SCHEME FOR DIFFUSION OF TRACERS : 2;2 causes a rapid temperature increase and a much earlier failure, around day 13.

The result is sensitive to MPI partitioning: some runs using 256 cores failed considerably earlier than comparable runs using 128 cores.

Restarting from a valid SERAFIND restart file can sometimes pass the time at which the continuous run failed.

The last observation suggests that a non-restarted internal field or accumulated numerical state may contribute to the failure.

Could this behaviour be caused by an interaction between:

tracer diffusion and the tidal open boundary;

inflow/outflow switching during tidal reversals;

the vertical boundary structure;

crushed vertical elements away from the boundary;

WAQTEL thermal forcing;

or density law 3 outside its usual salinity range?

Salinity exceeds 42 g/L during the dry season and can become much higher upstream. Maximum simulated temperatures also exceed 40 C, so I am checking whether these values may affect the density calculation or the baroclinic terms.

Which variables or TELEMAC routines would you recommend instrumenting to locate the first non-finite value before it reaches GRACJG ? Is there an available debug option that stops the computation immediately when the first NaN is created?

I can provide the steering file, and .out file.

Thank you for your help and sorry for the 3 posts

4

Hi

Just to gives you some hints for you tests, you could use a compiler option to stop on Nan. On gfortran, the option is

-ffpe-trap=invalid

This could be completed with

,zero,overflow

to check related floating exceptions…

To investigate deeper and in an easier way, I would recommend to fix the number of proc and the domain decomposition and only run the computation with those fixed pre-computation steps.

There is also an option to ensure a perfect restart and then avoid the problem you mention when restarts pass the time when the previous run failed…

Hope this helps

1 Like
5

Hello [name removed],

Thank you very much for your answer, it is very useful.

I will rebuild TELEMAC with -ffpe-trap=invalid,zero,overflow -g -fbacktrace and run the continuous simulation again. This way, the code will stop at the first invalid operation and show me the name of the routine. This is exactly what I needed.

I have one question about the perfect restart. In my configuration, I already use:

  • RESTART FILE FORMAT : SERAFIND

  • PREVIOUS COMPUTATION FILE FORMAT : SERAFIND

  • a restart from the RESTART FILE (not from the result file), with RECORD NUMBER FOR RESTART

  • the same number of processors (128)

Is there another keyword in TELEMAC-3D v9.1 to make the restart perfectly identical (bit for bit)? Or is this configuration already a perfect restart?

If it is already perfect, then the fact that my restarts pass the time where the continuous run fails could mean that my case is at the limit of stability and very sensitive to round-off. This would be consistent with the partitioning sensitivity you mentioned (128 versus 256 cores).

From now on, I will keep the same decomposition (128 processors, same split reused) for all comparisons.

I will report back when I have the name of the trapped routine.

[name removed]

6

Hi

you should also use RESTART MODE = TRUE

I don’t have v9.1 installed but I see on my side there is:

- RECORD NUMBER IN RESTART FILE


- RECORD NUMBER FOR RESTART

the function is different so have a look to the dictionnary (search for RESTART) to be sure to use the right one…

Regards