Troubleshooting

Sometimes the simulation can fail. It usually fails with an error about NaN values or negative pressure computed. The output looks like this:

NaN value or Negative AUX in Density_Pressure_Eqs, EXITING

In that case, the program will write out a file in the output folder, in the diagnostics subdirectory, that contains the grid position that failed. The output looks like this:

$ cat diagnostics/diag_density_pressure_eq_45053_1_0.txt
NAN_VALUE - idx: 94 - lev: 14
NAN_VALUE - idx: 94 - lev: 15                                          

There are two identified cases where this happens.

Low top of atmosphere pressure.

If the pressure at the top of the atmosphere gets to small (below millipascal), the simulation can fail in computing the next pressure step and ends with a negative pressure.

The simulation outputs the minimum pressure at each step and if it fails, you should see that the level where it fails is close to the top of the atmosphere.

To avoid this, either lower the top altitude of the grid, or increase the bottom pressure.

G_pm limiting and mu_star wiggle

In some cases, the different parameters in the computation of the G_plus and G_minus terms make it so that the denominator gets close to 0 and the G_plus and G_minus explode. This results in an unphysical step in the flux and an exploding Qheat value in one cell. This can then lead in some cases to negative pressures and a crash. This can be identified in the diagnostics file by seeing that usually some intermediate layers had bad values.

To mitigate this, the simplest is to slightly change the value of mu_star (i.e. change a bit the angle of the incoming beam), to avoid being in the zone of the discontinuity for the column that presents the issue. This is done by iteratively checking the G_± value and incrementing the incoming angle until the G_± value check passes.

There are some parameters that can be used.

The check can either be done in two ways:

• limit on the full G_± term, setting Alf_G_pm_limit_on_G_pm = true and Alf_G_pm_limit to something between 100 and 10000.
• limit on the denominator G_± term, setting Alf_G_pm_limit_on_G_pm = false and Alf_G_pm_limit to something around 1e-5.

The angle gets incremented by Alf_G_pm_mu_star_increment value (in degrees), usually something of 0.5 or 1.0 degree works. It increments the angle value and checks if the G_± criteria is met, and if not increments again, up to Alf_mu_star_iteration_max times, before giving up.

This shows up in this way in the output:

Hit G_pm denom limit, wiggling mu_star, loop counter 1
Hit G_pm denom limit, wiggling mu_star, loop counter 2
Hit G_pm denom limit, wiggling mu_star, loop counter 3
Hit G_pm denom limit, wiggling mu_star, loop counter 4
Hit G_pm denom limit, wiggling mu_star, loop counter 5
Hit G_pm denom limit, wiggling mu_star, loop counter 6
Hit G_pm denom limit, wiggling mu_star, loop counter 7
Hit G_pm denom limit, wiggling mu_star, loop counter 8
Hit G_pm denom limit, wiggling mu_star, loop counter 9
Hit maximum iteration of mu_star wiggle, bailing out

Sometimes, several levels in a column are close to the limit, so changing the angle might solve the issue on one cell but bring the neightbour cell close to the limit. To stringent criteria on G_± might make it so that it can’t really find a reasonable solution and changes the angle by a big value.

If it hits Alf_mu_star_iteration_max and gives up, it might still have a discontinuity and crash.

It seems that Alf_G_pm_limit_on_G_pm = false is less prone to crash, but can still lead to some discontinuities in the final flux, whereas Alf_G_pm_limit_on_G_pm = true with a small Alf_G_pm_limit value (below 300) can have difficulties finding a solution and still crash.

This needs a bit of experimenting to find the correct value.