Abnormally small time step on finest level

[q86091483]

Dear AMReX-Combustion community:

We are running a reactive jet-in-cross-flow case with PeleLMeX. After running the case with coarse levels for several characteristic time scales, we are ready to run with the finest level. The problem we met is that: on the finest level, the time step ($\mathrm{d}t$) decreases to abnormally small value (about $10^{-10}~ \mathrm{[s]}$ ), much smaller than what CFL condition would approximately give (about $10^{-8}~ \mathrm{[s]}$). This does not occur in coarser level runs. In our case the $\mathrm{d}t$ is limited mainly by advective CFL condition instead of CVODE. The CFL number is 0.8.

The following is the $\mathrm{d}t$ at each time step. The finest level is turned on after step 9400. After that $\mathrm{d}t$ can increase to $2\times 10^{-8}~ \mathrm{[s]}$. This value is consistent with the estimation based on CFL=0.8 and jet exit velocity. But evetually $\mathrm{d}t$ keeps decreasing to $10^{-10}~ \mathrm{[s]}$.

My questions are:

1. For each time step, how to locate and print the position and the magnitude of domain-maximum velocity? It looks like to me that the MultiFab class can return the domain-max/min value by ReduceMax/Min method but how to find the corresponding index, something like argmax(x_velocity) which returns $(x,y,z)$?

2. I noticed that there are different methods for calculating $\mathrm{d}t$, as specified by peleLM.divu_dt_method. By default only 1 method is adopted. What are their differences and why these different methods are needed?

3. We observed that for some cases when density decreases close to peleLM.divu_dt_rhoMin, the $\mathrm{d}t$ would randomly decrease to abnormally small values. What's the role played by peleLM.divu_dt_rhoMin in the algorithm? Can we avoid the small $\mathrm{d}t$ by decreasing peleLM.divu_dt_rhoMin?

4. Do you think the abnormally small $\mathrm{d}t$ occurs due to refinement? If so, is it possible to improve refinement criteria to avoid the abnormally small $\mathrm{d}t$?

thanks,
Zisen

[esclapez]

Hello Zisen,

Thank you for reaching out.

I'll try to answer you specific questions, but will first circle back on the dt estimate in LMeX.
There are two main dt limits in LMeX: 1) a CFL based and 2) an expansion (or divU) based. No Fourier limit since all the viscous terms are solved implicitly.

1. Within the non-subcyling approach of LMeX, the CFL limit will reduce dt by a factor two for every new level of refinement (assuming the velocity maximas do not change).

2. The divU-based estimate limits the dt based on ~1/divU times how far the density is from a "min" value peleLM.divu_dt_rhoMin, defaulted to 0.1. If the density drops below peleLM.divu_dt_rhoMin, a baseline ~1/divU of dt is used. The intent here is to trigger a dt reduction when you approach a min value of density that is deemed unreasonable in your case, e.g. resulting from numerical articfacts or else. Shortening the dt is expected to reduce the numerical artifacts. hopefully returning to a normal operation. Practically, this min value should be adapted to your case: if your case operates in a low pressure environment for instance, you might want to decrease peleLM.divu_dt_rhoMin. The effect of refinement on this dt estimate is not straightforward: refinement can introduce small scale structures with locally large negative divU that were "averaged" out previsously by the cell-averaged representation of the solution within a finite-volume scheme.

In practice, I have seen divU-based time step restrictions for low velocity (say Mach < 0.01) cases. If it pops-up on a high velocity case (which your case might be with a CFL dt around 2e-8), it is likely the result of some numerical artifact.

Can you give us a bit more details on your case ? If you use hydrogen with pure hydrogen jets, the density of hydrogen is low such that reducing peleLM.divu_dt_rhoMin might be advisable (what is the density at the center of the jet, away from mixing) ? Additionally, having large density gradients can lead to local extremas with LMeX default advection scheme (peleLM.advection_scheme=Godunov_PLM). Adding refinement can unfortunalety increase these local extremas (stronger Gibbs oscillations). If you are not using EB, can you try peleLM.advection_scheme=Godunov_BDS ? This scheme was specifically designed to prevent introducing maximas, but is a bit more computationnally expensive.

On your specific questions:
Q1 - There is nothing built-in, but it might be possible to add it.
Q2 - Method is the default, the second method was inherited from PeleLM. Maybe @drummerdoc can recall case where method 2 was more appropriate ?
Q3 - See discussion above
Q4 - See discussion above

Cheers,

Lucas

[drummerdoc]

The divU_dt limit can be derived easily as a time step limit coming from a density sink due to expansion, as Lucas suggests.

drho/dt + r.DivU = 0, discretizting, we can find the target deltat allowable before density drops below rhomin:

(rhomin - rho^o)/deltat = -0.5*(rhomin + rho^o)(divU^o + 0.5deltat*d(divU)/dt)

This leads to a quadratic equation for deltat,and only comes into play when the flow is expanding. It is the largest step you can take before losing too much mass from the cell due to flow leaving the through the cell walls.

When the solution is fairly coarse, the divU_dt limit typically kicks in just before something bad happens...the classic example is when a cell somehow ends up with lots of fuel and high temperature so that the reaction step creates too much heat release and things go bad. The other case, like Lucas points out is when the flow is slow, often laminar, and the CFL condition is just not limiting things enough to prevent this heat release issue. The code really like strong turbulence because of this, and when you're working with laminar systems, you often have to apply an artificial max dt to prevent these bad things.

[q86091483]

Hi Lucas and Mark,

Sorry for the late update.

the problem was solved by reducing cfl from 0.8 to 0.5. As for the peleLM.divu_dt_rhoMin issue, the nominal minimum density in my case is ~0.3 kg /m3, which is near the threshold.

thanks,
Zisen