open TELEMAC-MASCARET The mathematically superior suite of solvers

Hello,

I am upgrading a 3D model in a Mediterranean location from v6p3 to v7p2. This model comprises thermal exchanges at the air-sea interface which have been hardcoded in BORD3D as boundary conditions (so it is not exactly equivalent to the WAQTEL model n°5 where solar ray penetration is modelled as source term). Yet using identical parameterisations for both v6p3 and v7p2 models, the v7p2 simulation crashes halfway and gives NaN values for temperature, while v6p3 simulation runs fine (results are similar between both versions until the NaN values arise). I would like to understand why...

I attached a sequence of snapshots (dt=5s) showing surface temperature a few time steps before NaN appear in v7p2; water height for one time step is also shown. What we observe is the sudden apparition of "singular spots" where temperature gets either very large or very small, located in zones where water height is quite low (H<1m). In the following time step(s), the solver crashes and gives NaN.

I suspect diffusion is somehow responsible for this but not sure...

The precision for tracer diffusion is 10^-8 and the related maximum number of iterations is 200 (which is regularly exceeded before NaN the problem occurs).

Does anybody have any suggestion for understanding what is actually going on and fixing that issue?

Thank you very much,
JBS.

Here is the figure:

After some investigation, I have found this old discussion lead by JM Hervouet.

The problem might be caused by the relative numerical instability of the boundary condition:

nu_T*dT/dz = ATABOS*T+BTABOS

If this is the same cause, I would be happy to know more about this instability and the way it can be overcome (is it possible to make the formulation more stable?).

Thank you.
JBS.

Hello JBS

Before diving in the numerical aspects, can you make small preliminary investigations?
if your model was running well with v6p3 with hardcoded radiation effects, the problem could be linked to that part of the code.
You can for instance, change the option for radiation with the keyword FORMULA OF ATMOSPHERIC RADIATION (you have all implemented fomula in the module exchange_with_atmosphere.f in sources/waqtel).

The idea is to find the guilty part of source term (or boundary condition) that causes the crash. You can see that these terms are strongly non linear and could cause numerical issues depending on how they are tackled in time (implicitly or explicitly). By the way, you can also try to decrease the time step to see if it can overcome the problem.

with my best regards

Riadh

Thank you Riadh,

Sorry for the delay in replying.
It seems reducing the time step from 5s to 2.5s was successful, which is the most important!
The formulation for the boundary condition at the surface (highly non-linear indeed) is implicit/explicit (ie. ATABOS<>0 and BTABOS<>0) - the net vertical heat flux had been linearized in T to this end. I have tried a fully explicit formulation (ie. ATABOS=0, BTABOS<>0) but the simulation immediately crushed, whatever the time step.
I still don't know why 5s were ok in v6p3 and no longer in v7p2, but this might not be a great concern.

Thanks again,
JBS.