MetOffice / MetOffice/lfric_apps
Unphysical vertical wave amplification in non-Earth isothermal CRM setup
- Dominant language
- Fortran
- Stars
- 31
- Forks
- 118
- Avg merge
- 2d 4h
- Merged PRs (30d)
- 8
Description
### Version
main
### Are there any linked Issues or Pull Requests?
_No response_
### What happened?
LFRic version: 2.2
I am attempting to run CRM simulations of a dry H2-dominated atmosphere with an isothermal initialisation. Since there is no isothermal profile option for planar domain initialisation (see #211 ), I am using specified_profiles to set the initial profile manually as isothermal with variable gravity. I have turned off the shallow atmosphere approximation as something in that initialisation case seems to adjust the profile away from isothermal.
When changing the composition of the atmosphere to something lighter (i.e. increasing $R_d$ and $c_p$), a vertical wave forms that increases in amplitude. I have found the point where this subsequently causes the model to crash to be roughly around $R_d \sim 1000$ J/(kg K), $c_p \sim 3500$ J/(kg K) - but it also depends on other things such as uniform vs stretched extrusion methods. It seems that the greater $R_d$ & $c_p$, the quicker this wave appears, amplifies, and crashes the model.
Below I have plotted the central column of the domain for a run with $R_d=975.92$ J/(kg K), $c_p=3535.20$ J/(kg K). The extrusion method is `um_L140_122t_18s_40km` and I have set alpha=1 for fully implicit timestepping - something @Bshipway suggested.
For the same run, I have included the log for the timestep where the model crashes in the 'relevant log output' section.
My understanding is that LFRic should be invariant to changes to the gas parameters. I'd greatly appreciate some help figuring this out. I am using suite [u-dv345](https://code.metoffice.gov.uk/trac/roses-u/browser/d/v/3/4/5/trunk).
### Relevant log output
timestep 44
```shell
20260123173130.734+0000:P000:S44:INFO : Start of timestep 44
20260123173130.735+0000:P000:S44:DEBUG: Gungho: computing si_operators
20260123173130.736+0000:P000:S44:DEBUG: si_ops: mg level=1
20260123173130.787+0000:P000:S44:DEBUG: si_ops:MG restrict L1 to L2
20260123173130.787+0000:P000:S44:DEBUG: si_ops: mg level=2
20260123173130.804+0000:P000:S44:DEBUG: si_ops:MG restrict L2 to L3
20260123173130.804+0000:P000:S44:DEBUG: si_ops: mg level=3
20260123173130.811+0000:P000:S44:DEBUG: si_ops:MG restrict L3 to L4
20260123173130.811+0000:P000:S44:DEBUG: si_ops: mg level=4
20260123173130.819+0000:P000:S44:DEBUG: slow_physics: Running explicit JULES layer
20260123173130.823+0000:P000:S44:DEBUG: slow_physics: Running explicit Boundary layer
20260123173130.845+0000:P000:S44:INFO : Field Min/max dtheta slow = 0.00000000 0.00000000
20260123173130.856+0000:P000:S44:DEBUG: Mass matrix solve for W2 space.
20260123173130.856+0000:P000:S44:DEBUG: Created mass matrix operator
20260123173130.857+0000:P000:S44:DEBUG: Created diagonal pre-conditioner
20260123173130.858+0000:P000:S44:DEBUG: Gungho: mass matrix solve:
20260123173130.902+0000:P000:S44:INFO : Field Min/max lipschitz_3d = -490.36114502 160.70668030
20260123173130.902+0000:P000:S44:INFO : excess count = 1594 1455
20260123173130.902+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.902+0000:P000:S44:INFO : at level = 139 139
20260123173130.902+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.902+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.902+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.902+0000:P000:S44:INFO : Field Min/max lipschitz_2d = -490.29998779 160.70898438
20260123173130.902+0000:P000:S44:INFO : excess count = 1568 1469
20260123173130.902+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.902+0000:P000:S44:INFO : at level = 139 139
20260123173130.902+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.902+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.902+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.903+0000:P000:S44:INFO : Field Min/max lipschitz_x = -245.91545105 80.87847137
20260123173130.903+0000:P000:S44:INFO : excess count = 727 811
20260123173130.903+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.903+0000:P000:S44:INFO : at level = 139 139
20260123173130.903+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.903+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.903+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.903+0000:P000:S44:INFO : Field Min/max lipschitz_y = -244.38453674 81.62170410
20260123173130.903+0000:P000:S44:INFO : excess count = 1179 1141
20260123173130.903+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.903+0000:P000:S44:INFO : at level = 139 139
20260123173130.903+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.903+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.903+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.903+0000:P000:S44:INFO : Field Min/max lipschitz_z = -0.06114225 0.06581653
20260123173130.904+0000:P000:S44:INFO : Field Min/max lipschitz_sz = -0.03057112 0.03290826
20260123173130.904+0000:P000:S44:INFO : Field Min/max lipschitz_sx = -245.94602217 80.87732098
20260123173130.904+0000:P000:S44:INFO : excess count = 731 803
20260123173130.904+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.904+0000:P000:S44:INFO : at level = 139 139
20260123173130.904+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.904+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.904+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.904+0000:P000:S44:INFO : Field Min/max lipschitz_sy = -244.41510787 81.59958322
20260123173130.904+0000:P000:S44:INFO : excess count = 1185 1136
20260123173130.904+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.904+0000:P000:S44:INFO : at level = 139 139
20260123173130.904+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.904+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.904+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.904+0000:P000:S44:INFO : Field Min/max lipschitz_s2 = -490.33055892 160.70783398
20260123173130.904+0000:P000:S44:INFO : excess count = 1576 1463
20260123173130.904+0000:P000:S44:INFO : degeneracy = 1 1
20260123173130.904+0000:P000:S44:INFO : at level = 139 139
20260123173130.904+0000:P000:S44:INFO : at height = 3.85307780E+04 3.85307780E+04
20260123173130.904+0000:P000:S44:INFO : at latitude = 45.00000000 45.00000000
20260123173130.904+0000:P000:S44:INFO : at longitude = 0.00000000 0.00000000
20260123173130.905+0000:P000:S44:DEBUG: Min/max unity dep_dist_z, step 1 = -0.36197167E-01 0.90064573E+00
20260123173130.906+0000:P000:S44:DEBUG: Min/max unity dep_dist_xy_inner, step 2 = -0.40000000E+01 0.40000000E+01
20260123173130.906+0000:P000:S44:INFO : Transport: max horizontal FFSL Courant 0.40000000E+01
20260123173130.906+0000:P000:S44:INFO : Transport: max vertical FFSL Courant 0.18012915E+01
20260123173130.906+0000:P000:S44:INFO : Transporting density...
20260123173130.907+0000:P000:S44:DEBUG: Min/max field_np1: end_of_conservative_step = -0.13452171E+02 0.55662571E+02
20260123173130.912+0000:P000:S44:DEBUG: Min/max unity dep_dist_xy_outer, step 2 = -0.60000000E+01 0.60000000E+01
20260123173130.913+0000:P000:S44:DEBUG: Min/max field_np1: end_of_conservative_step = -0.37405038E+03 0.51859863E+03
20260123173130.914+0000:P000:S44:DEBUG: Min/max unity dep_dist_z, step 3 = -0.26128577E+03 0.17302687E+03
20260123173130.914+0000:P000:S44:DEBUG: Min/max field_np1: end_of_conservative_step = -Infinity Infinity
20260123173130.915+0000:P000:S44:INFO : Transporting wind...
20260123173130.917+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.20521445E+04 0.17954231E+04
20260123173130.920+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.26113070E+06 0.46123445E+05
20260123173130.920+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.28404100E+05 0.75427992E+05
20260123173130.921+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.18084467E+04 0.18266996E+04
20260123173130.923+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.37357473E+05 0.16078830E+06
20260123173130.923+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.36980902E+05 0.75643922E+05
20260123173130.924+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.98903418E+01 0.90326027E+02
20260123173130.926+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.38761514E+04 0.28332629E+04
20260123173130.926+0000:P000:S44:DEBUG: Min/max field_np1: end_of_advective_step = -0.98805801E+04 0.74775337E+04
20260123173130.930+0000:P000:S44:INFO : Transporting potential temperature...
20260123173130.931+0000:P000:S44:ERROR: The minimum of field ref_mass, step 1 is -0.38437868E+08 which is below the tolerance of 0.11920929E-06 which makes it impossible to find consistent departure distances. Most likely the transporting wind has a maximum Lipschitz number that exceeds 1. Stopping model rather than continuing
```
Contributor guide
No contributing guide indexed for this repository
Research direction
Reproduce the CRM case from suite u-dv345 using specified_profiles, the non-Earth isothermal setup, and the reported R_d and c_p values. Start with the timestep 44 log, especially the transport and Lipschitz diagnostics, and compare runs across extrusion methods and gas parameters. Done means identifying the cause of the growing vertical wave and preventing the simulation from producing infinities or crashing.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- fortran
- Domain
- hpc
- Issue type
- Bug
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Stale
- Clarity
- Needs clarification
- Newbie friendliness
- 25/100