Adaptive time step constraints not stict enough

[fjansson]

In tstep.f90, tstep_update() determines the largest dt that is safe to use.
In the model, there are two diffusion rates, ekm for momentum and ekh for other quantities (thl, qt). W.r.t. diffusion, the time step is currently limited only based on ekm. However, ekh is larger than ekm, and thus the time step limit is not always strict enough.
Some of our model crashes showed checkerboard-like disturbances in qt, of the kind that show up when a diffusion problem is solved with too large time steps. When I added ehk to the time step criterion, those problems disappeared.

peclettotl=max(peclettotl,maxval(ekm(2:i1,2:j1,k))*rdt/minval((/dzh(k),dx,dy/))**2) 
peclettotl=max(peclettotl,maxval(ekh(2:i1,2:j1,k))*rdt/minval((/dzh(k),dx,dy/))**2)   ! propose adding this

Note there are two places in the code where this should be done. This of course comes with the cost of smaller time steps, so the simulations take longer.

In another sense, the diffusion time step limit is sometimes too strict. It considers the smallest of (dx, dy, dz). When dz << dx, dy, this is stricter than necessary.

The current condition is

ekh * dt * (1/min(dx,dy,dz)**2) < 1/6

which could, I think, be relaxed to

ekh * dt * (1/dx**2 + 1/dy**2 + 1/dz**2) < 1/2

The constant 1/6 is actually the parameter peclet.

[Chiil]

My experience is that the 1/2 that comes out of the linear stability analysis is too permissive for LES, because of the large variations in the diffusivity in short distances. I agree with your statement, so maybe the combination of the both gives the best solution. In MicroHH we use: std::max(dnmul, std::abs(tPrfac*evisc[ijk]*(dxidxi + dyidyi + dzi[k]*dzi[k]))); Which is pretty much what you are suggesting. 1/2 leads to instabilities in the diffusion, for the course resolution simulations, especially those with clouds, a 1/6 value is already at the edge with centered advection schemes.

…

On 05 Jan 2018, at 13:49, Fredrik Jansson ***@***.***> wrote: In tstep.f90, tstep_update() determines the largest dt that is safe to use. In the model, there are two diffusion rates, ekm for momentum and ekh for other quantities (thl, qt). W.r.t. diffusion, the time step is currently limited only based on ekm. However, ekh is larger than ekm, and thus the time step limit is not always strict enough. Some of our model crashes showed checkerboard-like disturbances in qt, of the kind that show up when a diffusion problem is solved with too large time steps. When I added ehk to the time step criterion, those problems disappeared. peclettotl=max(peclettotl,maxval(ekm(2:i1,2:j1,k))*rdt/minval((/dzh(k),dx,dy/))**2) peclettotl=max(peclettotl,maxval(ekh(2:i1,2:j1,k))*rdt/minval((/dzh(k),dx,dy/))**2) ! propose adding this Note there are two places in the code where this should be done. This of course comes with the cost of smaller time steps, so the simulations take longer. In another sense, the diffusion time step limit is sometimes too strict. It considers the smallest of (dx, dy, dz). When dz << dx, dy, this is stricter than necessary. The current condition is ekh * dt * (1/min(dx,dy,dz)**2) < 1/6 which could, I think, be relaxed to ekh * dt * (1/dx**2 + 1/dy**2 + 1/dz**2) < 1/2 The constant 1/6 is actually the parameter peclet. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#29>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAGu6eikS6tNNshT0_jRoMBuUwgCbQDnks5tHhpAgaJpZM4RUaCc>.

[fjansson]

The issue of ekh and ekm is fixed in version 4.2, commit e9257ef.
The question of considering all of dx,dy,dz instead of the smallest is still open (i.e. the limit may now be stricter than necessary).