Wrong dycoms results with 1D IMEX

[smarras79]

Dycoms output is incorrect for 1D IMEX with either ARK2GiraldoKellyConstantinescu or AtmosAcousticGravityLinearModel using the value of the explicit stable dt.

The IMEX solution is here:

versus the EXPLICIT solution (closer to what is expected) is here:

Notice: for ease of comparison, just compare theta in the top right plots.

@fxgiraldo @mwarusz I used the following driver

[tapios]

The fluxes (e.g., ) look off. Do the fluxes include the SGS (Smagorinsky) contribution? If not, could you include it? And could you please also plot <w \theta_l>?

[tapios]

And could you run this with the exact same time step for the explicit and implicit method (i.e., not exploiting the computational efficiency IMEX affords)? I want to see whether stability limitations for the diffusion are the problem here.

[smarras79]

hi @tapios / @fxgiraldo running it with a fixed dt=0.01 (the max stable for explicit) using IMEX the solution is what we expect to have. I am currently at 2000s into the simulation and the stats look good. It seems that the issue above is tied to larger dt only. More tests are needed.

[smarras79]

At 14400 we see big differences that may be tied to either the flux bc or mixing. A plot is attached and we can discuss options tomorrow.

The plot of <w'theta_l> follows:

[tapios]

At least your subcloud-layer stays well mixed now (in q_t and theta_l). It seems to me the rest of the differences can plausibly attributed to the numerical differences, without explicit (physical or coding) errors. Or do you see differences that point to outright errors in the code?

[smarras79]

At least your subcloud-layer stays well mixed now (in q_t and theta_l). It seems to me the rest of the differences can plausibly attributed to the numerical differences, without explicit (physical or coding) errors. Or do you see differences that point to outright errors in the code?

Hard to tell. I am going through the code (b.c., turbulence, sources) of the current master that I used to run this to make sure everything is still consistent. I can't tell why a larger dt would cause mixing issues with IMEX though.

[smarras79]

I found what causes the problem described in this open issue. The wrong results occur when AtmosModel.jl does the following check:

if m.ref_state isa HydrostaticState
    flux.ρu += (p-aux.ref_state.p)*I
  else
    flux.ρu += p*I
  end

Without this check (i.e. simply usingflux.ρu += p*I) the IMEX solution is correct.

This issue is related to issue #591 .

[tapios]

Good point. The hydrostatic state should not be subtracted again when we use the linearization.

[smarras79]

Good point. The hydrostatic state should not be subtracted again when we use the linearization.

You say "again"; where is it currently being subtracted other than at the lines mentioned above?

[tapios]

Good point. The hydrostatic state should not be subtracted again when we use the linearization.

You say "again"; where is it currently being subtracted other than at the lines mentioned above?

We are linearizing about the hydrostatic state in the linear model. I'm unclear which version is currently used to compute the remainder model. Could you point me to it? The error likely occurs there.

[smarras79]

Good point. The hydrostatic state should not be subtracted again when we use the linearization.

You say "again"; where is it currently being subtracted other than at the lines mentioned above?

We are linearizing about the hydrostatic state in the linear model. I'm unclear which version is currently used to compute the remainder model. Could you point me to it? The error likely occurs there.

I am not sure. This may be more a question for @fxgiraldo and/or @simonbyrne But this file may be what you are looking for?

https://github.com/climate-machine/CLIMA/blob/master/src/Atmos/Model/remainder.jl

[tapios]

I'm not sure this is being used. Will check with @simonbyrne.

[tapios]

This is obsolete, isn't it? @smarras79 can we close it?

[smarras79]

Yep, obsolete. Closing it.