n-density dry model

[thabbott]

I have an implementation of an n-density dry model that seems to be working reasonably well using energy as the prognostic thermodynamic variable ...

Energy, one gas
Energy, three gases

... but that doesn't work particularly well with entropy ...

Entropy, one gas
Entropy, three gases

I think the issue is representing entropy production by mixing, which is hard to compute accurately. The rate of entropy production by mixing depends non-linearly on mass mixing rates and the entropy of a specific constituent blow up at low densities, and those two things together make me suspect that it will be hard to model mass diffusion when using entropy as the thermodynamic variable without producing spurious energy/pressure/temperature sources. (It's also possible there's a bug in the code I wrote...)

I'm labeling this as work in progress (and think we should hold off on merging for now) because the code is getting a little clunky. Some additions to the code I think would be useful include:

1. a more robust approach to tracking which model tracers are the thermodynamic variable and the densities (right now the code assumes that they're at certain indices of model.tracers)
2. a more streamlined way to include diagnostic fields in operators. Abstraction for velocity fields without extra memory allocations #59 would address this.
3. some way to enforce boundary conditions on diagnostic fields (e.g. pressure, diffusive fluxes). Right now this can be done indirectly by setting boundary conditions on prognostic fields, but having a more direct way to apply them would simplify the code (and let us avoid hacks like the one in diffusive_pressure_flux_z in compressible_operators.jl, which is currently there because I couldn't figure out any other way to ensure the flux was zero at the top and bottom boundaries). Abstraction for velocity fields without extra memory allocations #59 could also address this.

@ali-ramadhan, let me know if you're interesting in working on any of these additions. (Might be worth doing some pair programming so we make sure we're both happy with the resulting code.)

It also might be worth adding some regression tests before we start working on any of the above so that we have a way of making sure we're not changing model behavior (apart from fixes to any bugs we find) while refactoring.

Resolves #49
Resolves #52
Resolves #54
Resolves #55
Resolves #56
Resolves #57
Resolves #58
Resolves #59

[thabbott]

Also want to note that this branch has or will have changes that address #49 (but need to update the notes from that issue), #52, #54, #55, #56, #57, #58, and #59.

[ali-ramadhan]

Pretty epic PR! Super cool to see the three gas thermal bubble and the difference between energy and entropy.

For entropy, I'm not familiar with the spurious mixing issue but it seems really bad with the three gas thermal bubble. Even with the single gas case (both energy and entropy) you get numerical artifacts but that could just be the advection scheme. Quick test would be to see if increasing the resolution reduces those numerical artifacts.

But in general, good advection scheme seems pretty important.

Just skimmed through the changes but I won't review yet as it's a WIP. If there's a bug with entropy hopefully we'll find it during review or pair programming.

1. a more robust approach to tracking which model tracers are the thermodynamic variable and the densities (right now the code assumes that they're at certain indices of model.tracers)

Hmmm, yeah this might warrant us taking a slightly different approach than for incompressible models. I think a big benefit of having them be part of model.tracers is that they can be processed and time-stepped together.

One solution could be to define model.thermodynamic_variables and model.densities to be a named tuple of pointers to specific tracers in model.tracers. So we can treat them all similarly when calculating right hand sides or time stepping, but can pass in just model.thermodynamic_variables when needed, for example in a pressure calculation or microphysics scheme.

2. a more streamlined way to include diagnostic fields in operators. Abstraction for velocity fields without extra memory allocations #59 would address this.

True. Maybe model.tracers can hold the prognostic fields while properties like model.thermodynamic_variables and model.velocities can hold lazy primitive fields.

3. some way to enforce boundary conditions on diagnostic fields (e.g. pressure, diffusive fluxes). Right now this can be done indirectly by setting boundary conditions on prognostic fields, but having a more direct way to apply them would simplify the code (and let us avoid hacks like the one in diffusive_pressure_flux_z in compressible_operators.jl, which is currently there because I couldn't figure out any other way to ensure the flux was zero at the top and bottom boundaries). Abstraction for velocity fields without extra memory allocations #59 could also address this.

Ah sorry yeah the boundary conditions in master are still hacked in as I was hoping to include each field's boundary conditions in the Field struct (CliMA/Oceananigans.jl#606).

But yeah we definitely want to avoid hard-coding in boundary conditions like that!

Hmmm, this one seems like it could be a bit tricky if we really need boundary conditions and halo filling on diagnostic fields. Filling halos on the fly or on demand for lazy computations could be expensive.

It would be nice if we can get away with just imposing boundary conditions on prognostic fields so that lazy evaluations of the diagnostic variables will never have to worry about boundary conditions.

Hmmm, looking into diffusive_pressure_flux_z, I think it should be able to diagnose a zero pressure gradient if the thermodynamic variable has no-flux boundary condition (zero gradient or no-flux at the boundary), which should produce p/ρ = constant at the boundary. Then ∂zᵃᵃᶠ(i, j, k, grid, p_over_ρ_diag, tvar, g, Ũ, ρ, ρ̃, C̃) = 0. But that must have not worked as you had to hard code in the no-flux boundary conditions. Looks like halo regions are filled correctly during time stepping too.

Worse case even if we had to "fill in halo regions" for lazy diagnostic fields I think we can do it generally and there will be a performance hit but it should be rather negligible.

It also might be worth adding some regression tests before we start working on any of the above so that we have a way of making sure we're not changing model behavior (apart from fixes to any bugs we find) while refactoring.

True. I guess we don't have easy output writers built in yet (#38) but should be easy to output/load a few fields for the purposes of regression testing.

Been thinking of moving diagnostics and output writers out of the model and into a Simulation type (CliMA/Oceananigans.jl#447) so I haven't added an output_writers property field to CompressibleModel yet.

@ali-ramadhan, let me know if you're interesting in working on any of these additions. (Might be worth doing some pair programming so we make sure we're both happy with the resulting code.)

Pair programming could be a good idea, especially to settle on design choices.

I'll try to push through a couple of upstream issues: CliMA/Oceananigans.jl#606 so it'll be much easier to add boundary conditions on all fields here, and CliMA/Oceananigans.jl#447 so we can start using diagnostics and output writers (and get rid of the time stepping loop).

But maybe more importantly, I can work on adding lazy fields (#59) which seem like they might be a much needed abstraction for a CompressibleModel.

[ali-ramadhan]

Also want to note that this branch has or will have changes that address

I think if you include a list like

Resolves #49 
Resolves #52 
...

(but without the code block) in your original comment then those issues will be automatically closed when this PR is merged.

[thabbott]

Even with the single gas case (both energy and entropy) you get numerical artifacts but that could just be the advection scheme. Quick test would be to see if increasing the resolution reduces those numerical artifacts.

I think the solutions in this branch might be more oscillatory than master in part because diffusion is less strong---master has κ = 0.5 m^2/s, but with a bug in the diffusion operators that makes this equivalent to κ = 500 m^2/s in this branch, whereas the videos I posted use κ = 75 m^2/s. I'm also using half the spatial resolution. So we could make the verification experiments look nicer pretty easily (but having better advection operators eventually would be good too).

One solution could be to define model.thermodynamic_variables and model.densities to be a named tuple of pointers to specific tracers in model.tracers. So we can treat them all similarly when calculating right hand sides or time stepping, but can pass in just model.thermodynamic_variables when needed, for example in a pressure calculation or microphysics scheme.

Something like this seems good to me---I agree that keeping the thermodynamic variable and densities in the list of tracers is convenient, and this would give a straightforward and relatively low-cost way to access them when we need them individually.

True. Maybe model.tracers can hold the prognostic fields while properties like model.thermodynamic_variables and model.velocities can hold lazy primitive fields.

It might make sense to store all the lazy diagnostic fields in another struct called something like model.diagnostic_variables so that we clearly distinguish them from variables that get stepped forward in time. But we can think more about the best way to organize things.

Hmmm, looking into diffusive_pressure_flux_z, I think it should be able to diagnose a zero pressure gradient if the thermodynamic variable has no-flux boundary condition (zero gradient or no-flux at the boundary), which should produce p/ρ = constant at the boundary. Then ∂zᵃᵃᶠ(i, j, k, grid, p_over_ρ_diag, tvar, g, Ũ, ρ, ρ̃, C̃) = 0. But that must have not worked as you had to hard code in the no-flux boundary conditions. Looks like halo regions are filled correctly during time stepping too.

The problem comes from the fact that energy includes kinetic energy and halos for the vertical velocity field aren't filled the same way as other variables at the top and bottom boundaries. So even though the energy gradient is zero at boundaries the temperature and pressure gradients aren't. Calculating diffusive fluxes using lazy field with its own boundary condition might let us do something like

struct LazyGradient{B}
    scalar :: Field
    bc :: B
end

getindex(A::LazyGradient, inds...) = gradient(i, j, k, A.scalar.grid, A.scalar, A.bc)

function gradient(i, j, k, grid, scalar, bc::AllPeriodic)
   xcomp = ∂x(i, j, k, grid, scalar)
   ycomp = ∂y(i, j, k, grid, scalar)
   zcomp = ∂z(i, j, k, grid, scalar)
   return (xcomp, ycomp, zcomp)
end

function gradient(i, j, k, grid, scalar, bc::NoGradientZ)
   xcomp = ∂x(i, j, k, grid, scalar)
   ycomp = ∂y(i, j, k, grid, scalar)
   zcomp = (k == 1 || k == grid.Nz) ? 0.0 : ∂z(i, j, k, grid, scalar)
   return (xcomp, ycomp, zcomp)
end

This doesn't require extra halo filling but also handles the case where halo filling doesn't do what we'd like it to.

I'll try to push through a couple of upstream issues: CliMA/Oceananigans.jl#606 so it'll be much easier to add boundary conditions on all fields here, and CliMA/Oceananigans.jl#447 so we can start using diagnostics and output writers (and get rid of the time stepping loop).

But maybe more importantly, I can work on adding lazy fields (#59) which seem like they might be a much needed abstraction for a CompressibleModel.

Sounds good to me!

[ali-ramadhan]

Of course it took longer than expected but I've made some changes to Oceananigans.jl that will benefit JULES so I'll update JULES to use the latest version and start with building up some lazy fields.

A quick summary:

1. Grids now have a topology which help inform boundary conditions (Grids now have a topology CliMA/Oceananigans.jl#614, Create field boundary conditions using information from grid topology CliMA/Oceananigans.jl#620)
2. time_step! now takes a single time step. A high-level Simulation abstraction handles time stepping, diagnostics, and output writing (High-level Simulation type to manage time stepping CliMA/Oceananigans.jl#621).
3. Oceananigans.jl now implements an IncompressibleModel so there's room for a CompressibleModel (Rename Model to IncompressibleModel CliMA/Oceananigans.jl#626).
4. Boundary conditions are now a field property so JULES can readily support arbitrary boundary conditions (Upgrade fields to store boundary conditions CliMA/Oceananigans.jl#631).

[ali-ramadhan]

Resolves #49
Resolves #52
Resolves #54
Resolves #55
Resolves #56
Resolves #57
Resolves #58
Resolves #59

[codecov]

Codecov Report

Merging #60 into master will not change coverage by %.
The diff coverage is n/a.

@@           Coverage Diff           @@
##           master      #60   +/-   ##
=======================================
  Coverage   94.91%   94.91%           
=======================================
  Files          11       11           
  Lines         295      295           
=======================================
  Hits          280      280           
  Misses         15       15           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a96c8f8...a96c8f8. Read the comment docs.