Compressible atmospheric model

[ali-ramadhan]

Over at JULES.jl @thabbott, @RaphaelRR, and I have been building a compressible non-hydrostatic model on top of Oceananigans.jl for atmospheric simulations. It follows the numerics of the Advanced Research WRF Model as described in Klemp et al. (2007) and Skamarock et al. (2019).

With three qualitative verification experiments I think it's ready to be tagged as v0.1.

I'm opening this issue to discuss the possibility of merging JULES.jl into Oceananigans.jl which would add a CompressibleModel type (discussed a bit in #566). The two models share a lot of the infrastructure we've already built for Oceananigans, and in integrating the two we will readily get a GPU-ready research-grade compressible atmosphere model.

Verification experiments:

1. Hydrostatic adjustment in an isothermal atmosphere [very qualitative, see Bannon (1995)]
2. Dry rising thermal bubble: Dry rising thermal bubble verification experiment ali-ramadhan/Atmosfoolery.jl#31
3. Nonlinear density current: Nonlinear density current verification experiment ali-ramadhan/Atmosfoolery.jl#32

Dry rising thermal bubble can be compared with figure 5 of Wicker and Skamarock (1998), figure 7 of Jahn et al. (2015), and https://faculty.nps.edu/fxgirald/projects/mesoscale/rtb_movie.html.

Nonlinear density current can be compared with figure 1 of Straka et al. (1993) and https://faculty.nps.edu/fxgirald/projects/mesoscale/dc_movie.html.

Things we can hopefully share between IncompressibleModel and `CompressibleModel:

1. Operators
2. Grids
3. Coriolis
4. Forcing functions
5. Boundary conditions
6. Turbulence closures
7. Diagnostics
8. Output writers
9. Advection schemes
10. User interface (setting initial conditions, boundary conditions, time stepping, etc.)

Atmosphere-specific abstractions that would differentiate compressible from incompressible models:

1. Prognostic temperature (potential temperature, entropy, enthalpy, etc.)
2. Base state (dry adiabatic, isothermal, analytic sounding, form from initial conditions, etc.)
3. Acoustic time stepper (nothing, fully explicit for regular grids, vertically implicit for stretched grids)
4. Setting fields and initial conditions is more complicated because the density needs to be adjusted when setting the potential temperature. Also because the prognostic variables are ρ, ρu, ρv, ρw, ρθ, etc.
5. Other atmospheric things like cloud microphysics and radiation but for this there are packages that can be relied on, e.g. RRTMGP.jl and Cloudy.jl.

It would be ideal if we have a single Model type and just change the equation between the CompressibleModel and IncompressibleModel but the numerical methods and abstractions employed are sufficiently different that I think it makes sense to have two separate model types.

It might even make sense to keep a separate package for this CompressibleModel but I think for now it's better if we keep it all inside the Oceananigans.jl repository (see #284 for some previous discussion) until the code is more stable and also for extra visibility while the package is young.

References:
Bannon (1995): https://doi.org/10.1175/1520-0469(1995)052%3C1743:HALP%3E2.0.CO;2
Jahn et al. (2015): https://doi.org/10.5194/gmd-8-317-2015
Klemp et al. (2007): https://doi.org/10.1175/MWR3440.1
Skamarock et al. (2019): https://opensky.ucar.edu/islandora/object/opensky%3A2898
Straka et al. (1993): https://doi.org/10.1002/fld.1650170103
Wicker and Skamarock (1998): https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2

[johncmarshall54]

Excellent Ali. We could perhaps use the compressible model as a starting point to explore making it quasi-incompressible, moving it an ocean direction. This would lead to some interesting research papers. John

…

On Wed, Jan 29, 2020, 9:57 PM Ali Ramadhan ***@***.***> wrote: Over at JULES.jl <https://github.com/thabbott/JULES.jl> @thabbott <https://github.com/thabbott>, @RaphaelRR <https://github.com/RaphaelRR>, and I have been building a compressible non-hydrostatic model on top of Oceananigans.jl for atmospheric simulations. It follows the numerics of the Advanced Research WRF Model as described in Klemp et al. (2007) and Skamarock et al. (2019). With three qualitative verification experiments I think it's ready to be tagged as v0.1. I'm opening this issue to discuss the possibility of merging JULES.jl into Oceananigans.jl which would add a CompressibleModel type (discussed a bit in #566 <#566>). The two models share a lot of the infrastructure we've already built for Oceananigans, and in integrating the two we will readily get a GPU-ready research-grade compressible atmosphere model. Verification experiments: 1. Hydrostatic adjustment in an isothermal atmosphere [very qualitative, see Bannon (1995)] 2. Dry rising thermal bubble: ali-ramadhan/Atmosfoolery.jl#31 <ali-ramadhan/Atmosfoolery.jl#31> 3. Nonlinear density current (PR incoming) Dry rising thermal bubble can be compared with figure 5 of Wicker and Skamarock (1998), figure 7 of Jahn et al. (2015), and https://faculty.nps.edu/fxgirald/projects/mesoscale/rtb_movie.html. Nonlinear density current can be compared with figure 1 of Straka et al. (1993) and https://faculty.nps.edu/fxgirald/projects/mesoscale/dc_movie.html. Things we can hopefully share between IncompressibleModel and `CompressibleModel: 1. Operators 2. Grids 3. Coriolis 4. Forcing functions 5. Boundary conditions 6. Turbulence closures 7. Diagnostics 8. Output writers 9. Advection schemes 10. User interface (setting initial conditions, boundary conditions, time stepping, etc.) Atmosphere-specific abstractions that would differentiate compressible from incompressible models: 1. Prognostic temperature (potential temperature, entropy, enthalpy, etc.) 2. Base state (dry adiabatic, isothermal, analytic sounding, form from initial conditions, etc.) 3. Acoustic time stepper (nothing, fully explicit for regular grids, vertically implicit for stretched grids) 4. Other atmospheric things like cloud microphysics and radiation but for this there are packages that can be relied on, e.g. RRTMGP.jl <https://github.com/climate-machine/RRTMGP.jl> and Cloudy.jl <https://github.com/climate-machine/Cloudy.jl>. It would be ideal if we have a single Model type and just change the equation between the CompressibleModel and IncompressibleModel but the numerical methods and abstractions employed are sufficiently different that I think it makes sense to have two separate model types. It might even make sense to keep a separate package for this CompressibleModel but I think for now it's better if we keep it all inside the Oceananigans.jl repository (see #284 <#284> for some previous discussion) until the code is more stable and also for extra visibility while the package is young. References: Bannon (1995): https://doi.org/10.1175/1520-0469(1995)052%3C1743:HALP%3E2.0.CO;2 Jahn et al. (2015): https://doi.org/10.5194/gmd-8-317-2015 Klemp et al. (2007): https://doi.org/10.1175/MWR3440.1 Skamarock et al. (2019): https://opensky.ucar.edu/islandora/object/opensky%3A2898 Straka et al. (1993): https://doi.org/10.1002/fld.1650170103 Wicker and Skamarock (1998): https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#605>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKXUEQSTLZJ7MSNAXN7LJS3RAHUS3ANCNFSM4KNLLK5Q> .

[ali-ramadhan]

Yup it should be flexible enough to adapt it for ocean modeling. Might be worth thinking whether we need any additional abstractions.

One important point is that the numerical method used does not time step pressure, it is instead diagnosed from potential temperature or entropy. May require some thought if we want to modify the speed of sound.

[glwagner]

I think its a fine idea to merge, though it is interesting that the code will still be called Oceananigans :-D

I'm not sure I understand point 4 completely. The ocean model also has the concept of a model "state"; for example, ghost cell values and nonlinear diffusivities depend on the velocity and tracer fields, in general. We currently handle this only during time-stepping --- there is no guarantee, for example, that the diffusivities or ghost cell values are consistent with the velocity field if the velocity field is set by the function set! (in fact, the diffusivities are always inconsistent with the current velocity field due to the fact that they are calculated prior to taking a time-step). If we want to guarantee such a consistency, we can develop the concept of a model "state" and apply it to both ocean and atmospheric models (eg a function update_state!(model), or something similar).

What does the jargon "sounding" mean?

Merger means we need to think carefully about how to reduce boilerplate to minimize (within reason) the cost of maintaining two parallel models as we develop their shared subcomponents.

I wonder if splitting off the output writers and diagnostics into a Simulation type that wraps AbstractModel may help. In this paradigm, a "Model" is reduced to numerics + physics specification. This would be easy to implement (while I think an Equation abstraction would be relatively difficult due to myriad difficult design problems, including the abstraction of tuples of terms with heterogeneous function signatures, numerical aspects of the equation, implicit vs explicit treatment of terms, etc).
The Simulation type can then be used to "run" simulations over multiple time steps, eg

simulation = Simulation(model, Δt=1.0, end_time=8day, output=output_writers, diagnostics=diagnostics)
run!(simulation)

and is tasked with managing things like output writing, diagnostics calculation, adaptive time-stepping, and progress logging. Each Model then simply needs to define a function that performs a single time-step to interface with Simulation. This is discussed further in #447. Such a orthogonalization of the code means we can develop the Simulation abstraction without having to worry about updating each AbstractModel individually.

[johncmarshall54]

Agree we need another name, perhaps a variant of shenanigans.

…

On Thu, Jan 30, 2020, 3:31 PM Gregory L. Wagner ***@***.***> wrote: I think its a fine idea to merge, though it is interesting that the code will still be called Oceananigans :-D I'm not sure I understand point 4 completely. The ocean model also has the concept of a model "state"; for example, ghost cell values and nonlinear diffusivities depend on the velocity and tracer fields, in general. We currently handle this only during time-stepping --- there is no guarantee, for example, that the diffusivities or ghost cell values are consistent with the velocity field if the velocity field is set by the function set! (in fact, the diffusivities are *always* inconsistent with the current velocity field due to the fact that they are calculated *prior* to taking a time-step). If we want to guarantee such a consistency, we can develop the concept of a model "state" and apply it to both ocean and atmospheric models (eg a function update_state!(model), or something similar). What does the jargon "sounding" mean? Merger means we need to think carefully about how to reduce boilerplate to minimize (within reason) the cost of maintaining two parallel models as we develop their shared subcomponents. I wonder if splitting off the output writers and diagnostics into a Simulation type that wraps AbstractModel may help. In this paradigm, a "Model" is reduced to numerics + physics specification. This would be easy to implement (while I think an Equation abstraction would be relatively difficult due to myriad difficult design problems, including the abstraction of tuples of terms with heterogeneous function signatures numerical aspects of the equation, implicit vs explicit treatment of terms, etc). The Simulation type can then be used to "run" simulations over multiple time steps, eg simulation = Simulation(model, Δt=1.0, end_time=8day, output=output_writers, diagnostics=diagnostics) run!(simulation) and is tasked with managing things like output writing, diagnostics calculation, adaptive time-stepping, and progress logging. Each Model then simply needs to define a function that performs a single time-step to interface with Simulation. This is discussed further in #447 <#447>. Such a orthogonalization of the code means we can develop the Simulation abstraction without having to worry about updating each AbstractModel individually. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#605>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKXUEQSOSWBJGZJKTTSGCZ3RALQCRANCNFSM4KNLLK5Q> .

[thabbott]

Nice work with the test cases @ali-ramadhan! No objections to merging from me, but I think we should have a discussion (in person or in a separate issue) about ways to separate aspects of the model that make it an atmosphere vs. an ocean model from things that make it a compressible vs. incompressible model. (The Boussinesq solver in Oceananigans is already useful for atmospheric simulations with depths small relative to a scale height and could be great for cloud modeling if there's any easy way to convert it to an anelastic solver.)

[ali-ramadhan]

I think its a fine idea to merge, though it is interesting that the code will still be called Oceananigans :-D

Agree we need another name, perhaps a variant of shenanigans.

True. The window of opportunity to change names is closing though. It might not be good to change package name after the JOSS paper comes out.

I'm open to a name change if we can all agree on a new name.

I'm not sure I understand point 4 completely. The ocean model also has the concept of a model "state"; for example, ghost cell values and nonlinear diffusivities depend on the velocity and tracer fields, in general.

Right. But I guess users don't have to worry about ghost cell values and nonlinear diffusivities as they'll get filled in or calculated as needed.

But in JULES.jl, if you want to add a potential temperature perturbation and modify θ you need to modify both the density ρ and the ρθ field. Example:

    θ = Θ[i, 1, k] / ρ[i, 1, k] + θ′(xC[i], 0, zC[k])
    π = Π(i, 1, k, grid, gas, Θ)


    ρ[i, 1, k] = pₛ / (Rᵈ*θ) * π^(cᵥ/Rᵈ)
    Θ[i, 1, k] = ρ[i, 1, k] * θ

From a user's perspective this is much more complicated than using Oceananigans' set! function. To further complicate things, how ρ is computed will change if you're running with entropy as a prognostic variable instead of potential temperature.

If we want to guarantee such a consistency, we can develop the concept of a model "state" and apply it to both ocean and atmospheric models (eg a function update_state!(model), or something similar).

This sounds like a good idea and would be of great benefit to users.

What does the jargon "sounding" mean?

Ah sorry I'm referring to atmospheric soundings. From talking to @thabbott apparently it's pretty common to initialize a model using observational data obtained from e.g. a weather balloon.

Apparently some models get pretty fancy and construct a hydrostatic base state from the sounding to avoid transient acoustic waves when running the model.

I wonder if splitting off the output writers and diagnostics into a Simulation type that wraps AbstractModel may help.

Hmmm, I wasn't super into the idea of a new Simulation type before but with multiple model types it makes more sense and would simplify the Model type by moving out output writers and diagnostics.

I still think we don't need a new type and can simply define a new function like

run_model!(model, Δt=1.0, end_time=8day, output=output_writers, diagnostics=diagnostics)

which just needs time_step! to be defined.

Nice work with the test cases @ali-ramadhan! No objections to merging from me, but I think we should have a discussion (in person or in a separate issue) about ways to separate aspects of the model that make it an atmosphere vs. an ocean model from things that make it a compressible vs. incompressible model.

Thanks! Hmmm, that's a good point that currently the IncompressibleModel is really more of an OceanModel as evidence by the surface_waves functionality and the SeawaterBuoyancy default.

And the proposed CompressibleModel is already more of an AtmosphericModel as tracers representing mixing ratios like ql, qv, qi, etc. are given special treatment (they contribute to the moist density).

The Boussinesq solver in Oceananigans is already useful for atmospheric simulations with depths small relative to a scale height and could be great for cloud modeling if there's any easy way to convert it to an anelastic solver.

My understanding (maybe @jm-c or someone else can confirm for the MITgcm) is that the current Boussinesq solver can be repurposed for anelastic simulations but it would have to be in pressure coordinates.

Either way, you raise a good point that solvers may be repurposed to model different systems/fluids.

We can meet in person to discuss (with @glwagner if he's also interested). I'll also open a new issue to discuss this on GitHub as we should also document these design choices.

[glwagner]

But I guess users don't have to worry about ghost cell values and nonlinear diffusivities as they'll get filled in or calculated as needed.

That's not quite true for the ordinary definition of "as needed". Halos and nonlinear diffusivities are only filled or calculated immediately prior to a time-step! If you set! the model state (or even if you time-step the model forward), you cannot actually calculate, diagnose, or plot 1) the gradients of fields across boundaries or 2) the nonlinear diffusivities.

In fact, it hadn't occurred to me previously but with our current paradigm the LES diffusivity field that is output is actually inconsistent with the velocity field even if they are output at the same model.clock.time. So at the moment this is a legitimate concern.

Are the "atmospherigans" being discussed here of the same form? Are some of these fields related diagnostically to one another, or are they only related diagnostically along a model boundary? (Or something --- its hard for me to generalize from the example given above.)

We may want to develop the concept of a model "state". We can provide a function update_state! that updates all fields (eg halos + fields that depend diagnostically on others like nonlinear diffusivities and possibly others), and use this both in set! as needed, during time-stepping, and prior to output. We can also create a system for 'tagging' components of the model that have update_state! methods attached to them. For example, some turbulence closures require the computation of perturbations from a horizontal average (eg, the buoyancy perturbation). Such fields would also need to be updated during an update_state! call. The lack of such an abstraction is currently a barrier to the implementation of those turbulent closures.

[glwagner]

What's the status of this issue?

[ali-ramadhan]

I'll close this issue since it's essentially resolved with PR #1079 which we decided to not merge as it will live in a different repo: https://github.com/ali-ramadhan/Atmosfoolery.jl is it's current foster home.