Bizarre blow-up issue, seemingly due to tracers

[gli353]

Hi,
I have run into a situation not similar to any posted problems on this forum (as far as I can see). The model blows up after about 2,000 steps, with the reason being PEchar = NAN; the potential energy, hence the tracers, or the depth, have gone crazy (while KEchar is a valid number; please see the attached log file, for the model log and my cpp flags). I know that others in this forum have seen PEchar = NAN as well, yet they concurrently also had KEchar = NAN. Given the mode-split nature of ROMS, I wonder if this simply means that the problem occurred for the stepping of the tracers from step n to n+1 (the <<ubar,vbar>> step in the diagram : https://www.myroms.org/wiki/Numerical_S ... _Technique)
I have done an autopsy of the very last snapshot stored in the rst file, and all to my surprise, here is not a single water grid cell that has a NaN tracer value! I have paid particular attention to the tracer values near my river input, yet they seem ok in the rst file as well, no NaN there either! As for the water depth, here is not any NaN in water grid either. (I did the above checking in matlab using the command
find(maskr(:)==1 & isnan(salt(:))), where maskr is the rho-mask in the grid file, not the wet-and-dry mask hence the wetting-and-drying grids are also included in this checking, and salt is the salt value from a particular sigma layer. Other variables are checked in the same way).

What is the possible reason for this blowup? Thanks in advance.

Kind Regards,
Gaoyang

[wilkin]

You should look for NaNs irrespective of the masks because ROMS computes on all points, and then applies the mask. A NaN could have arisen in the land prior to masking.

You don't have any explicit horizontal viscosity activated (UV_VIS2 is undefined) yet you have VISC_3DCOEFF, so I suspect those viscosities are not applied. I also notice some odd characters printed in coefficient ranges report.

Code: Select all

 Minimum horizontal diffusion coefficient =  0.00000000E+00 ､路Q
 Maximum horizontal diffusion coefficient =  0.00000000E+00 ､路Q

 Minimum horizontal viscosity coefficient =  0.00000000E+00 ､路Q
 Maximum horizontal viscosity coefficient =  0.00000000E+00 ､路Q

Not sure what to make of that.

Your river inflow salinities are in the range 16 to 22, not zero. Is that intentional?

The report on bathymetry shows you have tall topography on land ... h = -244 m. I don't have compelling evidence on this, but I would be wary of having negative depths in the permanent land adjacent to points that might wet or dry. It's kind of asking for trouble.

Other things to try are increasing the minimum water depth (you have 1 cm) for the wet/dry switch.

[gli353]

Hi Wilkin,
Thanks for the prompt reply!
1. I totally forgot about this. Yet, from the very initial salinity field, NaN values have populated the land points. Below is what I mean, the snapshot of the initial salinity field, as it is, without any mask being applied (so the NaNs have been there from the very beginning). And so is the initial velocity field stored in the first history file (the first time layer in roms_his_0001.nc): seems like the stored versions of the variables have been masked before being written to nc file? Not sure if I've got anything wrong here, and not quite certain how I may check (without entering the debug mode) if any land points' values have turned NaN during the simulation.

InitialSalt.JPG (26.47 KiB) Viewed 9568 times

2. That odd character there began to appear after an update of the HPC centre of an update of my WinSCP software. It came out of no where, yet was not in the log file a while ago. However, indeed I did not switch on UV_VIS2 in this simulation. It didn't improve the model, and tended to crash the model (even though dt*nu/(dx)^2 < 1/2 was not violated) so I undefined it.

3. Yes, the river salinity was set intentionally like that, the salinity is decreased gradually from these numbers down to 0. The initial salinity field has salinities between 16 to 22 at the headwaters, and so do the river flow injection there.

4. The water depth is derived from several DEM products... I thought I had set the very negative numbers to something sensible, though I will ensure this is done properly. I think this might also be the reason why applying harmonic viscosity would cause blowup (point 2) even though the numeric stability criterion dt*nu/(dx)^2 < 1/2 was not violated, that this could be due to applying mixing near a high dry land point, as discussed in many other posts?

5. Dcrit is actually 0.1m, not 1cm, as shown in the log

" 1.6000E-02 Zob Bottom roughness (m).
0.0000E+00 Zos Surface roughness (m).
1.0000E-01 Dcrit Minimum depth for wetting and drying (m)."

Thanks a lot for these suggestions!
Gaoyang

[gli353]

Also here, a snapshot of the bathymetry, without being masked:

Bathy_NoMask.JPG (44.67 KiB) Viewed 9563 times

The region-of-interest is in the black circle, and the very negative numbers are in the red circle, the two are very far-apart. The bathymetry in the black circle is based on nautical chart(s) and some LiDAR dataset, and I have made sure that the mountain range (Waitākere) has been smoothed out in the model. I will try if further removing the negative values in the lower right corner could really help...

Thanks,
Gaoyang

[wilkin]

I recommend that at all points in permanent land mask the depth be given a modest positive value so that they play nicely with adjacent wet point at the coast.

Within ROMS internal calculations you don't ever want NaNs anywhere. Is it possible that after some time wetting/drying dipped into a region where your initial conditions had NaN? But you have #define ANA_INITIAL so where did they come from?

The _FillValue in output files converts the land point values to NaN when you read them into e.g. Matlab, but they are not NaN inside ROMS. So, yes, it can be hard to track. You might be able to circumvent the conversion to NaN by reading files with the lower level netcdf.open routines in Matlab, before _FillVale, add_offset and scale_factor are applied.

[gli353]

1. I am doing this; and I have just further smoothed the grid and opted for VStretching = 4 (originally was 3) and see how it goes...

2. Aha, the story is like this. I decided to use the initial salinity field (representative of the April Waitematā) derived from some old reference and our regular monitoring. Yet the initial zeta and velocities are still set to zero (so I started out with the zero velocities and the estimate initial salinity field being the initial condition and did a warm-up run, with TS_FIXED; after a simulation of 10 days, when the residual current structure has been spun up for a spring-to-neap transition condition, I restart the model there). I did not want to prepare a large Initial file whose only spatially-varying variable is salt, so I hacked the ROMS code a bit, to make the ANA_INITIAL.h set initial zeta and velocities to zero, yet read the initial salinity field from file (the attached figure shows that the warm-up model reads in an extra variable Sini from my grid file). I thought here could be some bizarre values in my Sini variable and the wetting and drying, yet after checking I have confirmed this is not the case.

3. I will try these low-level functions and update what I find there in the rst file. Unfortunately, the original rst file storing the information related to this blowup incident was accidentally overwritten. I will later replicate the incident and try netcdf open. I believe here must be a NaN in a certain land point... At least, I have used these functions to view my Initial Fields: no NaNs there, not in initial velocities, nor in tracers, nor in zeta. Only the desired values and the _FillValue show up there.

Thanks!
Gaoyang

[gli353]

Thanks to Dr. Wilkin's suggestions, I have solved this issue. It was indeed due to some errors in my initial salinity field. After interpolation, I forgot to set the land point values to _FillValue, which later on caused troubles at some land points.... Case closed

I will now investigate how to add some horizontal tracer diffusion into the model; previously it would give me a blow-up after a short run of merely tens of steps; not sure if it was related to the wrong initial field or to the wetting-and-drying.