ali-ramadhan/atmosfoolery.jl

Compressible non-hydrostatic model built on top of Oceananigans.jl so it runs on CPUs and GPUs

I'm an Applied Scientist at Afresh working on reducing food waste, a huge contributor to climate change, by using machine learning to help grocery stores optimize their inventory management and ordering decisions.

Before that I received my PhD in Computational Earth, Atmospheric, and Planetary Sciences from MIT where I worked at the intersection of climate science and scientific machine learning. As part of the Climate Modeling Alliance I developed Oceananigans.jl, a fast and flexible next-generation ocean model written in Julia that runs on GPUs, and used it to train machine learning models of geophysical turbulence and simulate all kinds of fluid dynamics.

And before that I received my MS and BS in Physics from the University of Waterloo where I used ultrashort pulse lasers and synchrotrons to make movies of molecular dynamics.

You can find out more about me and read my blog/ramblings on my personal website: http://aliramadhan.me/

atmosfoolery.jl's Issues

Do we really need halos of size 2?

Overriding Base.getindex(A::OffsetArray{T,N}, I::Vararg{Int,N}) where {T,N} to print a stacktrace if it's in the second halo point,

function Base.getindex(A::OffsetArray{T,N}, I::Vararg{Int,N}) where {T,N}
    any(I .< 0) && display(stacktrace())
    @boundscheck checkbounds(A, I...)
    J = map(parentindex, axes(A), I)
    @inbounds parent(A)[J...]
end

I get three results which makes sense but can we rewrite so that the strain rate tensor doesn't dip into the second halo point?

 getindex(::OffsetArray{Float64,3,Array{Float64,3}}, ::Int64, ::Int64, ::Int64) at REPL[74]:2
 getindex at abstract_field.jl:172 [inlined]
 ℑxᶠᵃᵃ at interpolation_operators.jl:8 [inlined]
 U_over_ρ at compressible_operators.jl:15 [inlined]
 ∂xᶜᵃᵃ at derivative_operators.jl:18 [inlined]
 strain_rate_tensor_ux at compressible_operators.jl:171 [inlined]
 viscous_flux_ux at compressible_operators.jl:182 [inlined]
 δxᶠᵃᵃ at difference_operators.jl:21 [inlined]
 ∂ⱼτ₁ⱼ at compressible_operators.jl:193 [inlined]
 FU at right_hand_sides.jl:8 [inlined]

 getindex(::OffsetArray{Float64,3,Array{Float64,3}}, ::Int64, ::Int64, ::Int64) at REPL[74]:2
 getindex at abstract_field.jl:172 [inlined]
 ℑyᵃᶠᵃ at interpolation_operators.jl:11 [inlined]
 V_over_ρ at compressible_operators.jl:16 [inlined]
 ∂yᵃᶜᵃ at derivative_operators.jl:21 [inlined]
 strain_rate_tensor_vy at compressible_operators.jl:172 [inlined]
 viscous_flux_vy at compressible_operators.jl:183 [inlined]
 δyᵃᶠᵃ at difference_operators.jl:24 [inlined]
 ∂ⱼτ₂ⱼ at compressible_operators.jl:199 [inlined]
 FV at right_hand_sides.jl:13 [inlined]

 getindex(::OffsetArray{Float64,3,Array{Float64,3}}, ::Int64, ::Int64, ::Int64) at REPL[74]:2
 getindex at abstract_field.jl:172 [inlined]
 ℑzᵃᵃᶠ at interpolation_operators.jl:14 [inlined]
 W_over_ρ at compressible_operators.jl:17 [inlined]
 ∂zᵃᵃᶜ at derivative_operators.jl:24 [inlined]
 strain_rate_tensor_wz at compressible_operators.jl:173 [inlined]
 viscous_flux_wz at compressible_operators.jl:184 [inlined]
 δzᵃᵃᶠ at difference_operators.jl:27 [inlined]
 ∂ⱼτ₃ⱼ at compressible_operators.jl:205 [inlined]
 FW at right_hand_sides.jl:17 [inlined]

JULES name already taken

@milankl pointed out Hydro-JULES: https://hydro-jules.org/

There's also https://github.com/milankl/Juls.jl which is pretty close.

But probably don't need to pick a name until we want to register the package or something.

Halo regions from the beginning?

Do we want to support halo regions from the start?

If so OffsetArrays.jl might be useful.

Ignoring the messy question of how to implement boundary conditions, we found using halo regions significantly simplifies the spatial operators as the x, y, z operators can be made to all look the same.

Does it make sense to support non-zero bulk viscosities?

Somewhat of a followup from #64 but does it make sense to support a non-zero volume/bulk viscosity μᵥ?

It's suggested that μᵥ = 0 by the Stokes assumption is found to be accurate in many applications. So I'm guessing this effect is probably negligible for the atmosphere as it's approximately incompressible?

But it seems like μᵥ can be quite large for ideal gases: https://aip.scitation.org/doi/full/10.1063/1.4729611

Note: In #58 I think I mistakenly refer to the ⅓μ∇(∇⋅u) term as the bulk viscosity when it is not.

---

Some helpful material from Kundu et al. 6th edition:

Update viscosity (and diffusivity?) values in verification experiments.

Straka et al. (1993) use a constant isotropic diffusivity of 75 m²/s for the density current test case (no word on viscosity).

Can't find the value for the dry rising thermal bubble but will update with the value when I find it.

Bug in diffusion operators?

I think the diffusive fluxes in compressible_operators.jl:

@inline diffusive_flux_x(i, j, k, grid, κᶠᶜᶜ, ρ, C) = κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * δxᶠᵃᵃ(i, j, k, grid, C_over_ρ, C, ρ)
@inline diffusive_flux_y(i, j, k, grid, κᶜᶠᶜ, ρ, C) = κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, grid, C_over_ρ, C, ρ)
@inline diffusive_flux_z(i, j, k, grid, κᶜᶜᶠ, ρ, C) = κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * δzᵃᵃᶠ(i, j, k, grid, C_over_ρ, C, ρ)

@inline function ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, κˣ, κʸ, κᶻ, ρᵈ, C)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, diffusive_flux_x, κˣ, ρᵈ, C) +
                                    δyᵃᶜᵃ(i, j, k, grid, diffusive_flux_y, κʸ, ρᵈ, C) +
                                    δzᵃᵃᶜ(i, j, k, grid, diffusive_flux_z, κᶻ, ρᵈ, C))
end

are missing a division by the grid spacing. (Without it the units don't work out.) This might explain when we've been finding that the verification experiments require a smaller diffusivity than is typically reported in publications.

Implement flux calculation

(based on meeting on 9/16)

Some comments, suggestions, and typos on the algorithm notes

Been reading the algorithm notes and had some comments. Feel free to ignore any super pedantic points.

I'm looking at algorithm.pdf at this particular commit: https://github.com/thabbott/JULES.jl/blob/b884c353830a60d85283570da76c66bb10dcd9f7/notes/algorithm.pdf

1. The prognostic variable is not momentum, but momentum per unit volume right?

2. The index gymnastics identities given on page 1 aren't really used, should we still keep them?

3. Should we add Coriolis terms (2Ω×u or f×u) to equations (2)-(4)?

4. Could probably condense equations (2)-(4) into a single equation using u_i, \tau^{(i)} (or tensor \tau), and adding a unit vector to the buoyancy term if we want to. This would also get rid of the Cartesian-looking derivatives so the equation remains coordinate-independent with summation notation. Also, do we want to include forcing terms here?

5. Equation (5) is an empty equation.

6. Might be good to define specific entropy as the entropy per unit mass? Also what s_0, c_v, T_0, R, and \rho_0 are. Maybe this isn't the right time to make the documentation more accessible though.

7. Should we include the equation for arbitrary tracers in section 1, equations (6)-(10)? Presumably they satisfy an advection-diffusion equation like the conservation law defined at the beginning of section 1 plus a forcing term? Is it common to run with arbitrary tracers in the atmosphere for e.g. gas concentrations? I guess in the ocean it's mostly used for biogeochemical tracers.

8. More of a question but are J, τ, and ϵ all given by all SGS turbulence closures? I'm only familiar with LES closures but they usually employ eddy viscosity and eddy diffusivity models and give τ via an eddy viscosity and presumably J through an eddy diffusivity. Do they give ϵ? That seems more something that comes out of a k-ϵ turbulence closure.

9. In section 2, seems confusing to use R to refer to the right-hand-side as R is already used as the specific gas constant in the definition of specific entropy. We also have the same problem... MITgcm uses G to refer to the RHS, which seems a little cryptic...

10. Equation set integrated on the acoustic time steps (page 3) needs minus signs on ∂y(p′) and ∂z(p′).

11. δτ for small acoustic time step might be confusing as it implies infinitesimal to me. Why not use Δτ like Klemp et al. (2007)? I also like their superscript primes, but it's definitely tricky notation.

12. Equation (12) should have a -∂y(p′) not -∂x(p′).

13. I wonder if superscript n is better than superscript t to indicate a discrete time. t is usually used for continuous time. I guess Klemp et al. (2007) use superscript t. Depends if we want to consistent with their notation.

14. On page 1, specific entropy is defined with an R, which changes to a c_p in section 2.1. Defining R_specific = c_p - c_v might make this clearer.

15. Seems to me that linearization is not neccessary to solve for p' analytically. Just curious why it's being done.

AppVeyor code coverage

The AppVeyor build keeps failing while trying to detect coverage statistics---see e.g. this build. Since the Travis builds successfully report coverage should we just have the AppVeyor build run tests but not detect and report coverage?

Stress tensor divergence terms are incomplete?

Maybe I'm misunderstanding how the stress tensor divergence is being computed but looking at the stress tensor components written out explicitly in matrix form (see third posted image), and assuming the volume viscosity μᵥ = 0 and the viscosity is constant, then the x-component of the divergence of the stress tensor should be ∂ⱼτ₁ⱼ = μ∇²u + ⅓μ∂x(∇·u⃗) but the following operators seem to be missing terms (e.g. shouldn't viscous_flux_ux include a ∇·u⃗ term and not sure how the factor of 1/3 ended up in the strain rate tensor).

https://github.com/thabbott/JULES.jl/blob/3dac4ab6bc7d7f3a9f84c6240aa90dd365936ff5/src/Operators/compressible_operators.jl#L170-L196

---

Some helpful material from Kundu et al. 6th edition:

Adding support for base states

Not sure how important this is given that we're not super worried about truncation or rounding errors in the vertical momentum equation right now. From my understanding, the errors are more severe when running with terrain-following coordinates so perhaps we don't even base states to do most of what we want to do.

A big reason to add base states though is to provide hydrostatic states in which we can initialize simulations without worrying about transient acoustic waves (see issue #42).

It would also reduce the amount of boilerplate code by a lot. Currently the verification experiments spent a non-trivial amount of code setting up a dry adiabatic atmosphere and running it to steady state. Users cannot be expected to do this.

Integrate `Oceananigans.OutputWriters` with JULES.jl

This will help a lot with verification experiments (#33) and integration.

Milestones until we can submit a JOSS paper?

X-Ref: CliMA/Oceananigans.jl#605

Opening this issue to keep track of what needs to be done before we think JULES.jl and its CompressibleModel can be merged into Oceananigans.jl.

We just released Oceananigans v0.38.0 yesterday whiles JULES.jl is relying on Oceananigans v0.24.1 so I'll do some updating first then see what needs to be done.

---

Roughly in order:

• Get JULES.jl up to date with the latest version of Oceananigans.jl: #71
• Clean up the code by using a time stepper abstraction: #72
• Clean up verification experiments by making use of output writers and analysis scripts: #33
• Get a high-order advection scheme from Oceananigans.jl to work with JULES.jl: #75
• Get an LES closure from Oceananigans.jl to work with JULES.jl: CliMA/Oceananigans.jl#654
• Ensure you can run a compressible non-atmospheric verification test case: #69
• Rewrite JULES.jl time-stepping kernels to use KernelAbstractions.jl: #74
• Write some documentation: #73

Probably want JULES to be running on the GPU before merging.

Implement acoustic time steps

(based on meeting on 9/16)

Inertia-gravity waves verification experiment

Seems like a popular test and an easy one to implement.

See:

1. https://www2.mmm.ucar.edu/projects/srnwp_tests/IG_waves/ig_wave.html
2. https://faculty.nps.edu/fxgirald/projects/mesoscale/igw_movie.html

Abstracting away p′ and the equation of state

@thabbott suggested I take a look at section 2 and see if we can abstract away the calculation of p′ and the equation of state as they might depend on terms we don't know.

Just starting to think about this, and it seem like it'll be easy for the dry ideal gas law, but had some questions as I'm not familiar with atmospheric thermodynamics.

1. For atmospheric modeling, you always want to get an expression for the specific entropy? Presumably if you cannot obtain such an expression, you cannot use that equation of state.

2. The entropy of an ideal gas is given by the Sackur-Tetrode equation (probably not useful here) while the specific entropy defined in the notes is the change in entropy for an ideal gas undergoing some process. Would it be useful to denote it as Δs? Maybe not if it's universally understood that we always talk about entropy changes.

3. Will p′ ever depend on any other prognostic variables besides ρ′ and (ρs)′? I suppose for a moist atmosphere, water vapor might be a prognostic variable too?

4. Seems that c_p and c_v are introduced in the definition of specific entropy. Do we imagine having other equations of state that do not have a c_p or c_v, or have a wildly different definition of specific entropy? If so, we'll have to abstract away a lot of terms, mostly p′ probably.

Might be helpful to at least consider the dry ideal gas and moist ideal gas cases to see what abstraction would work best.

Are there other equations of state we should consider when designing this abstraction? @thabbott suggested maybe ideal gas for 2+ condensible components, but also maybe an equation of state where the ideal gas law breaks down?

Maybe we should just stick to the dry and moist ideal gases for now, and just refactor to support exotic equations of state when we want to implement them.

Still any interest?

Hey @thabbott and @RaphaelRR, been thinking about this project recently as there's no atmospheric LES model in Julia (or any GPU-enabled one although I may be wrong here) so I think it's a chance to make something easy to use and very useful. As well as tackle some cool research problems that would be within reach with such a code.

MIT is getting a pretty big GPU cluster at the end of the month so we should have plenty of local resources.

It's definitely a time-consuming project, and may be difficult to carry out as a side project though. Might be good to be realistic if the time commitment is too large. But we should discuss if there's still any interest, I have some ideas.

Integrate `Oceananigans.BoundaryConditions` with JULES.jl

Right now models are hard-coded to be horizontally periodic with no-flux and no-normal flow at the top and bottom.

Would be good to integrate Oceananigans' boundary conditions framework with CompressibleModel.

I personally think it would be cool if fields just carried around their boundary conditions (CliMA/Oceananigans.jl#606) but that's a non-trivial refactor and may or may not happen.

How to initialize an atmosphere in (discrete) hydrostatic balance?

I tried solving for pressure and density profiles that would satisfy the discrete form of hydrostatic balance but I still get transient acoustic waves when running the model.

Not super high priority but might be good to figure out how to do this to avoid transient acoustic waves.

Currently I just put a strong damping sponge layer at the top and all the acoustic waves are damped after ~500 seconds then I turn off the sponge layer.

Sandbox script: https://github.com/thabbott/JULES.jl/blob/master/sandbox/solve_for_hydrostatic_profiles.jl

Non-atmospheric test case: compressible Kelvin–Helmholtz

Might be a good test case to demonstrate that JULES.jl can do more than just atmospheric setups

Public vs. private repo?

Just curious on what everyone's thoughts are on making the repository public vs. private. I noticed the repo is currently private but not sure if that was on purpose.

Viscosity μ (Pa·s) vs. kinematic viscosity ν (m²/s)?

Just to confirm, we want to use kinematic viscosity ν (m²/s)?

it seems that the atmospheric modeling community always means kinemtic viscosity (m²/s) when they say "viscosity".

Account for bulk viscosity

We should probably be including the effects of the bulk viscosity in the momentum equation: the ⅓μ∇(∇⋅u) term.

See: https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Compressible_flow

Probably doesn't matter for the simpler verification tests we've carried out. Also not sure if it's included (or accounted for) with LES. Will have to investigate more.

Implement time stepper

(based on meeting on 9/16)

General questions

Took another attempt at reading through the algorithm notes and Klemp et al. (2007) [K07] but had some questions as I'm thinking about how we can code this all up.

Apologies for the spam of questions! Not expecting a reply on GitHub!

1. I think we agreed to use height coordinates. I assume we want to start with flat-bottom Cartesian coordinates, i.e. ζ(x, y, z) = z following K07?

2. I probably just misunderstood things but the advection term in K07's Eq. (4)-(9) seems to suggest that advection terms are made up of three terms, while the CM1 manual suggests six terms in their Eq. (6)-(7).

3. Seems that you can't run a model with just U, V, W (and no buoyancy) as a temperature-like variable such as entropy s or modified virtual temperature θₘ is required to diagnose the pressure and evolve U, V, W?

4. What is the minimal set of prognostic variables you can evolve? If a temperature-like variable is required then it seems to be U, V, W, S, and ρ?

5. Does it make sense to design the code such that the prognostic temperature-like variable can be changed in the future? We've chosen entropy s but maybe there's a reason to choose modified virtual temperature or enthalpy?

6. Setting initial conditions seems to be non-trivial because users usually want to specify velocities u, v, w and temperature T. But these are not the prognostic variables. Am I right in assuming that velocities get converted into momenta U, V, W and the temperature gets converted to entropy s, after which the initial density is computed?

7. A "time-invariant hydrostatically balanced dry reference state" is used to define perturbation variables as deviations from this reference state. Is there only one way to define this base state? The base state actually enters in the acoustic time steps so presumably changing the base state will change the dynamics?

8. Is the reference surface pressure p₀ user-defined? Is it a constant or is it p₀ = p₀(x, y)? It seems like changing p₀ would change the dynamics. Not sure if a 2D field makes sense as it implies topography which is taken care of by ζ(x, y, z). Any changes in surface pressure would be accounted for in p′.

9. While entropy s and the mixing ratios qv, ql, qi, etc. are all tracers, it seems they should be separated into difference structs as the temperature-like prognostic variable effects how both pressure and buoyancy are computed, while your choice of moist species affects your cloud microphysics and you can run without any moist species. Seems to make sense that passive tracers (e.g. ) should be further separated into a different struct.

10. Am I right in assuming that a cloud microphysics scheme just picks a number of moist species to evolve (expressed as mixing ratios) and adds terms to the right-hand-side of their evolution equations to specify how they interact?

11. Do moist species contribute to buoyancy solely through the moist density ρ_m?

12. If we wanted to evolve the acoustic modes explicitly, would this make sense as the smallest modification on the vertically implicit method? Not sure if it's more stable/accurate if we calculate W″(τ+Δτ) first, or if we calculate S″(τ+Δτ) and ρd″(τ+Δτ) first. Maybe they're both bad.

    1. Calculate U″(τ+Δτ) and V″(τ+Δτ) using variables from τ.
    2. Calculate W″(τ+Δτ) using S″(τ) and ρd″(τ).
    3. Calculate S″(τ+Δτ) and ρd″(τ+Δτ) using U″(τ+Δτ), V″(τ+Δτ), and W″(τ+Δτ).

13. How would user-defined forcing terms work? I see the potential for three different kinds of forcing terms:

    1. "slow" non-RK3 physics evaluated as part of the F terms in Eq. (4)-(8) of K07. These apparently "contain the Coriolis and diffusion terms as well as diabatic effects and any parameterized physics that arise in a full atmospheric prediction model". Seems like an LES diffusivity would be calculated here?
    2. "fast" RK3 physics evaluated as part of the R terms in Eq. (13)-(16) of K07.
    3. Not in K07 but you could potentially add a forcing term that is evaluated at every acoustic time step, i.e. on the left hand side of K07's Eq. (13)-(16).

14. I'm confused by Figure 3.1 of the WRF4 manual between steps (1) and (9). Step (1) is described as computing "physics tendencies including mixing" during the first RK step while step (9) is "compute non-RK3 physics (currently microphysics)". So it seems that microphysics enjoys a special status as it's computed separately at the end of every large time-step. Seems to me they can be combined?

15. Do we envision supporting topography one day? Would that involve switching to the terrain following coordinates ζ(x, y, z)? We've been thinking of an immersed boundary method for Oceananigans where the boundary is accounted for by adding forcing terms to the momentum equations so no refactoring is required and the pressure solver doesn't have to change.

Integrate `Oceananigans.TurbulenceClosures` with JULES.jl

Right now Laplacian diffusion and viscosity is hard-coded in with μ = κ = 0.5 which is really bad to say the least.

Should integrate with Oceananigans.TurbulenceClosures to be able to use the basic constant diffusivity closures as well as all the LES closures: Smagorinsky-Lilly, AMD, etc.

Turbulence closures in Oceananigans return kinematic viscosities while I think the compressible equations of motion generally use regular viscosity. I'm sure they can be meshed together, just gotta keep this difference in mind.

Cataloging different atmospheric test cases

Just throwing together some ideas.

Dry two-dimensional

• Gravity current or negatively buoyant bubble: Straka et al. (1993)

Moist two-dimensional

• Saturated buoyant bubble: Bryan & Fristch (2002)

Three-dimensional moist LES

• Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II RF01): Stevens et al. (2005)
• Barbados Oceanographic and Meteorological Experiment (BOMEX): Siebesma et al. (2003)
• Deep convection? Hohenegger and Bretherton (2011)

Other ideas

• A standard test set for nonhydrostatic dynamical cores of NWP Models, Bill Skamarock

The name `JULES.jl` may be too generic

Just a note --- this repository is difficult to find with a google search. It might be a good idea to find a more specific and informative name.

Stop hard coding in gravity [and allow for g(z)]

Right now g = 9.80665 is hard coded in which is bad.

In Oceananigans the gravitational acceleration is stored inside the buoyancy struct, but here I think I'll make it a separate model property which will also allow us to add a height-varying g(z) in the future if we want it for crazy stuff like simulating Jupiter's atmosphere.

Should we use PkgTemplates.jl to boostrap a Julia package?

https://github.com/invenia/PkgTemplates.jl

It'll generate all the files we need and set up YAML configuration files for CI testing (Travis, Appveyor, etc.) that should come in handy.

Worth retaining modified potential temperature as a prognostic variable?

It's convenient if you want to work with dry air + dilute water vapor but doesn't extend in a particularly straightforward way to other situations, so it's unclear to me whether it's worth keeping it in JULES (as opposed to just developing the model with entropy for now). Thoughts?

Reference pressures

Should the reference surface pressure (e.g. used in the definition of potential temperature) and the reference pressure in the equation of state always be the same?

If so, then it would make sense to keep track of the reference surface pressure in the IdealGas struct.

Thermodynamic reference states need re-organizing

The way we're diagnosing temperature and pressure from entropy assumes that all massive tracers have an entropy defined with the same reference temperature, but the way that reference states are stored (one for each massive tracer) doesn't enforce that, since each massive tracer has its own reference temperature. (This is something we could address when we add support for massive tracers that aren't non-condensible ideal gases, since this will require some substantial additions to the model thermodynamics.)

Simplest working model for v0.1?

Been thinking about the simplest working model we can use to demonstrate a working 3D compressible atmospheric, i.e. something we can call v0.1.

What do you guys think?

1. Uniform grid (already have)
2. Second-order operators (already have)
3. RK3 time stepping with explicit acoustic sub-stepping (@thabbott seems to almost have this done?)

The other nice features, e.g. vertically stretched grid, WENO advection schemes, implicit vertical solvers, etc., would all have to be developed and tested which might be a lot to do for v0.1 with our time budgeet.

Need a `RungeKutta3` time stepper abstraction

Right now the CompressibleModel struct contains four collections of fields that clutter things up. Seems like it would be good to combine them into a time stepper struct.

             slow_forcings :: S
          right_hand_sides :: R
    intermediate_variables :: I
     acoustic_time_stepper :: W

Questions about setting up DYCOMS verification experiment

The basic state for RF01 was compiled from all of the measurements and is idealized as a quasi-two-layer structure in liquid water potential temperature θₗ and total-water specific humidity qₜ...

1. For the n-density model I'm guessing we want to add a second gas for water vapor with water vapor density ρₗ and water vapor total energy eₗ?

We specify geostrophic winds of U = 7 m/s and V = 5.5 m/s, which produce winds within the PBL near 6 and 4.25 m/s"

2. How do we best impose this? Presumably you specify the geostrophic winds and expect a small shear to develop? Two ideas from CliMA.jl include geostrophic forcing and Rayleigh damping.

For the sea surface temperature we specify a value of 292.5 K, which is 2.1 K warmer than the surface air temperature. Given a bulk aerodynamic drag coefficient, CD = CH =
CQ = 0.0011, this should correspond to a surface sensible heat flux near 15 W/m² and a surface latent heat flux of approximately 115 W/m².

To test the degree to which the assumption of fixed fluxes masked differences among the simulations, additional simulations were performed by most groups for which surface temperatures were held fixed and surface fluxes were computed interactively. However, these simulations did not differ substantially from those with surface fluxes fixed at the above values, and so the results reported on below are from simulations with specified fluxes.

3. Should be easy enough to enforce a couple of fluxes but I'm not super sure how to enforce them. Is the latent heat flux converted into a water vapor energy flux? And the sensible heat flux is converted into a total energy flux?

4. The radiative cooling is converted into an appropriate energy flux?

Generalized split-explicit time-stepping

I've seen people talking about Knoth & Wensch (2014) recently as a nice alternative to Wicker-Skamarock time-stepping. They report being able to take time steps 4 times as large which seems pretty impressive (although not sure if they came up with this for turbulent 3D flows).

I think they just replace they generalize the method so you can pick something else for the fast integrator, and they suggest a bunch of multirate infinitesimal step (MIS) methods.

Probably not something we want to implement until we have a working model and time to explore, but might be a good idea to keep the code flexible enough to switch between fast integrators.

Abstraction for velocity fields without extra memory allocations

Right now getting velocity is a little clunky as you have to do something like

model.momenta.ρu[i, j, k] / model.density[i, j, k]

Would be nice to have an abstraction for velocity fields like

struct LazyVelocityField
    momentum:: Field
    density :: Field
end

getindex(A::LazyVelocityField, inds...) = A.momentum[i, j, k] / A.density[i, j, k]

then we can just write model.velocities.u[i, j, k] and it'll lazily calculate the velocity.

The powerful feature (besides zero memory allocations) is that we can pass this lazy velocity field around to operators and other functions just like a regular field.

Could also use it to tracers (instead of mass-weighted tracer densities [?] if that's fields like ρθ are called) and it could replace the use of functions like U_over_ρ and C_over_ρ.

Some concerns about the Klemp (2007) equation set and suggestions for transitioning to a moist model

I want to discuss what I see as a potential problem with following the Klemp (2007) paper exactly (which is where JULES.jl currently seems to be heading) as we work on converting from a dry to a moist model, which is that the way Klemp (2007) uses dry density prevents its equation set from working in a limit where the entire atmosphere is made of water. This is admittedly a pretty academic problem, but for the sake of flexibility I would rather use an equation set that doesn't break down for extremely condensible-rich atmospheres.

My suggestions is that, moving forward, we use the Klemp (2007) paper as a useful guide but base our model a slightly different set of equations. Something like equations 24-27 and 31-32 from Satoh (2003) (plus an appropriate thermodynamic equation) would be a more robust set of governing equations, since they still make sense even with zero dry density. The key differences from Klemp (2007) are that Satoh (2003)

1. tracks the mass of water species as concentrations relative to total mass rather than mixing ratios relative to dry mass, and
2. defines ρu as the total mass flux, not just the flux of dry mass

One intermediate goal we could work toward that would give us a chance to start playing with equation sets with multiple massive species but without all of the complexity of moist thermodynamics would be to work on writing code that integrates a set of equations for two non-condensible ideal gases (call them a and b):

(I'm using the notation from the Satoh paper---so the Ds are diffusive fluxes, and the F is viscous dissipation. I've put question marks on the RHS of the entropy equation since I need to think more about what should go there---presumably a source term corresponding the momentum sink and probably a mixing terms of the two gases aren't identical.)

This would give us a chance to check the robustness of the equation set (e.g., does everything still work with ρ_a = 0?) and play with how best to implement a model with multiple massive species in JULES, all without having to deal with a ton of complicated thermodynamic bookkeeping. If you think this sounds like a reasonable next step I can work on closing this equation set (and replacing entropy with some other thermodynamic variable, if you'd prefer) to give us something you can actually work on implementing.

Add higher-order advections schemes, e.g. WENO-5

WENO-5 is implemented in Oceananigans.jl v0.38.0 and seems to work well so we should see if we can get JULES.jl to use the higher-order advection schemes from Oceananigans.jl

Some high-level documentation of the numerical methods

I believe the numerics of the n-density model implemented by @thabbott is described by Satoh (2003). Would be good to write some documentation about this.

Satoh, Conservative Scheme for a Compressible Nonhydrostatic Model with Moist Processes,
Monthly Weather Review (2003) 131 (6): 1033–1050: https://journals.ametsoc.org/mwr/article/131/6/1033/66966

Mixing ratios and moist density

Right now we've mostly been running in dry mode so it hasn't mattered but mixing ratios contribute to the moist density so the model needs to know what's a mixing ratio.

Two approaches we can take:

1. Mixing ratios are just tracers but depending on the "microphysics scheme" being used, specific mixing ratios will contribute to moist density. This is similar to how we currently treat active tracers and passive tracers (they're both lumped into model.tracers).
2. Add a new property model.mixing_ratios and sum over all mixing ratios when calculating moist density.

I'm gonna go with 1 for now to avoid further complicating the model but we may want to revisit this issue in the future.

CFL condition utilities

Will be useful especially to quantify and verify speedups with acoustic time stepping.

To calculate CFL condition due to acoustic waves, should probably use https://en.wikipedia.org/wiki/Euler_equations_(fluid_dynamics)#Ideal_gas

Use KernelAbstractions.jl

Should be easy to get on the GPU then. Maybe we should limit our use of fancy abstractions until we're running on the GPU (where we can make sure said fancy abstractions can be made to work efficiently on GPUs).

Add analysis scripts to verification experiments

Once we integrate output writers with JULES.jl then the verification experiments should be split into two scripts:

1. A simulation script that runs the verification experiment.
2. An analysis script that makes figures and movies.

Example: https://github.com/climate-machine/Oceananigans.jl/tree/master/verification/stratified_couette_flow

Differentiate between potential temperature and modified potential temperature?

So far we haven't really done any moist simulations so the prognostic temperature called ModifiedPotentialTemperature is really just potential temperature.

With moisture added Klemp et al. (2007) define a "modified potential temperature" that's a little different from virtual temperature.

I wonder if we should have both PotentialTemperature (used in dry atmospheres) and a ModifiedPotentialTemperature (used when moist species are present) that the user should pick from.

Batched tridiagonal solver?

@thabbott suggested looking into tridiagonal solvers for the implicit (ρw)′ solve.

I imagine a tridiagonal system needs to be solved for at every column? If so might be good to develop a batched solver.

It'll be trivial on the CPU (just loop over all columns) but might be good to model our batched solver on https://gist.github.com/maleadt/1ec91b3b12ede9898958c95596cabe8b so we can easily extend it to the GPU when the time comes.

In that gist @maleadt was able to develop some really fast batched tridiagonal solvers on the GPU. I think he was able to perform 256^2 tridiagonal solves (256 vertical levels each) in 1.5 ms. Might have been because cyclic reduction is really good for problem sizes that are powers of 2.

---

Aside showing my ignorance:

Is the vertical treated implicitly because of vertically propagating waves? I thought the acoustic sub-stepping takes care of acoustic waves in all three dimensions.

Not clear to me why an implicit solve is required for high-resolution LES where the grid cell aspect ratio is close to 1, unless this is not the use care being planned for.

Adding pressure as a prognostic thermodynamic variable

Using pressure as a prognostic thermodynamic variable has fallen out of fashion in mesoscale atmospheric models (you generally won't conserve energy to numerical precision if you don't use it as a prognostic variable) but apparently it's necessarily in compressible ocean models since the equation of state is generally too complicated to diagnose pressure from a different prognostic thermodynamic variable. Depending on the level of interest in eventually using the dynamic core from this model for ocean simulations it might be worth adding it as an option, even if it's not the best choice for atmospheric simulations.

Update JULES.jl to work with the latest version of Oceananigans.jl (v0.40.0)

Lots of good stuff we can use here.

Override getproperty to access thermodynamic variable

Would be nice to override getproprty(CompressibleMode, ...) to be able to write things like model.energy, model.entropy, etc. based on the thermodynamic variable in use.