Better organization of rotation-related and buoyancy-related constants and parameters

[glwagner]

As discussed with @jm-c, the organization of physical constants and parameters is somewhat confusing.

Currently, constants are stored in three places:

1. PlanetaryConstants, which stores a rotation rate, gravitational acceleration, and a Coriolis parameter used in an f-plane approximation
2. ModelConfiguration, which stores anisotropic (potentially turbulent) viscosities and diffusivities
3. EquationOfState, of which there is only one kind: LinearEquationOfState, which stores both parameters associated with the equation of state in addition to a reference density

I see a few problems:

• f is not a property of a planet.
• 'Model configuration' is an obscure name for turbulent or molecular transport coefficients.
• A reference density is not a parameter in an equation of state.

I propose that we consolidate these three types into two, removing the reference density from EquationOfState and define a new type containing f, g, ρ0, ν, and κ.

I'm not sure what to call the new type. One possibility is FluidParameters or PhysicalParameters or PhysicalConstants.

I also propose that we cease support for anisotropic transport coefficients as parameters, defined generally, at least for the moment. We can support constant anisotropic transport coefficients as a type of LES closure in the future.

[johncmarshall54]

Greg, I think that rho_o belongs in the EOS. And we will have to run many problems with anisotropic diffusion etc. So we should support from the outset. John

…

On Thu, May 9, 2019, 5:55 PM Gregory L. Wagner ***@***.***> wrote: As discussed with @jm-c <https://github.com/jm-c>, the organization of physical constants and parameters is somewhat confusing. Currently, constants are stored in three places: 1. PlanetaryConstants, which stores a rotation rate, gravitational acceleration, and a Coriolis parameter used in an f-plane approximation 2. ModelConfiguration, which stores anisotropic (potentially turbulent) viscosities and diffusivities 3. EquationOfState, of which there is only one kind: LinearEquationOfState, which stores both parameters associated with the equation of state in addition to a reference density I see a few problems: - f is not a property of a planet. - 'Model configuration' is an obscure name for turbulent or molecular transport coefficients. - A reference density is not a parameter in an equation of state. I propose that we consolidate these three types into two, removing the reference density from EquationOfState and define a new type containing f, g, ρ0, ν, and κ. I'm not sure what to call the new type. One possibility is FluidParameters or PhysicalParameters or PhysicalConstants. I also propose that we cease support for anisotropic transport coefficients as parameters, defined generally, at least for the moment. We can support constant anisotropic transport coefficients as a type of LES closure in the future. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#217>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AKXUEQW3GJUNC2CYLVWD44TPUSMUVANCNFSM4HL6JVQA> .

[ali-ramadhan]

Thanks for bringing these up! This stuff certainly needs some work, and we should probably make it clear what's e.g. a free parameter.

I see a few problems:

• f is not a property of a planet.

Agreed. Perhaps now that we're also thinking of channel models on a β-plane, we should also build a "rotation" abstraction to choose between f-plane, β-plane, and Coriolis force (possible with cosine term).

• 'Model configuration' is an obscure name for turbulent or molecular transport coefficients.

Yes. This should be addressed by incorporating a TurbulentDiffusivity struct like the one you proposed in #120. I agree with John that when we do this, isotropic and anisotropic diffusion should be options.

• A reference density is not a parameter in an equation of state.
  Hmmm, but if ρ₀ is needed to calculate ρ then isn't it a parameter of the EOS?

Sounds like this issue is worth discussing and strategizing about. We could maybe get some ideas and inspiration from how CliMA.jl is handling parameters?

[glwagner]

I support a Rotation abstract type with subtypes like FPlane, BetaPlane, and FullCoriolis (you are referring to including the horizontal components of the rotation vector within the tangent plane approximation?)

Do you think that we should not have the concept of a molecular transport coefficient at all? I feel we should have something like MolecularTransportCoefficients for things like molecular viscosity and different molecular diffusivities for each tracer (so that we support the most basic and easiest simulation case: DNS). TurbulenceClosure can handle all the cases in which the transport coefficients are non-molecular, as in the case of a constant anisotropic viscosity.

My comment about the reference density was actually my attempt to channel @jm-c and I think I mangled it. The point is that density is not a variable in a Boussinesq code and the value of the reference density has no effect on the dynamics. Buoyancy is the Boussinesq variable. It gets confusing when you add nonlinear equations of state, because you will still need to store a reference density somewhere in order to calculate boundary fluxes of temperature and velocity associated with, for example, heat fluxes and momentum fluxes. For that reason it makes sense to divorce the reference density from the equation of state and consider it separately. As long as you have a Boussinesq code you will not need to calculate density. You'll get to this point in a hurry if you try to implement nonlinear equations of state for your Europa simulation. Maybe @jm-c can weigh in so we can get these thoughts from the source.

[ali-ramadhan]

you are referring to including the horizontal components of the rotation vector within the tangent plane approximation?

Ah I should have been clearer. For the Cartesian domains we might want FullCoriolis but including the cosine term would only make sense if we're running on a sphere.

Do you think that we should not have the concept of a molecular transport coefficient at all?

I think it makes sense to have this, especially so that the default is a DNS. Not sure if this is the most common use case, in which case maybe we can spit out a warning to new users so they know what they're running.

The point is that density is not a variable in a Boussinesq code and the value of the reference density has no effect on the dynamics.

This makes sense to me, but then I'm wondering if ρ₀ is important to know to calculate diagnostics etc. where should we store it? EOS makes the most sense to me as a storage place even though it doesn't influence the dynamics.

[glwagner]

including the cosine term would only make sense if we're running on a sphere.

I'm not sure what you mean by cosine term --- I am referring to including the full rotation vector, which projects onto the horizontal as well as the vertical in a tangent plane approximation. Keeping only the vertical component of the projection is known as the 'traditional' approximation. For weakly stratified non-hydrostatic scenarios there is evidence that the full Coriolis force should be used. See:

• https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/nearinertial-waves-in-the-ocean-beyond-the-traditional-approximation/8E5E79284F842F15A79E1570D1B5A2A7

• https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2015GL066472

For the sphere, a different abstraction for Rotation should probably be used, because there is no longer the concept of the projection of the rotation vector into a coordinate system that is locally tangent to the Earth's surface at a particular latitude (which is what gives rise to the parameters f).

I think if this code is successful, DNS will be a common use case. There are no julia DNS CFD codes, and the are very few GPU DNS codes in general.

This makes sense to me, but then I'm wondering if ρ₀ is important to know to calculate diagnostics etc. where should we store it? EOS makes the most sense to me as a storage place even though it doesn't influence the dynamics.

This is the question. If the EOS actually calculates buoyancy rather than density, I feel that it is an error to associate ρ₀ with the EOS. Including ρ₀ in the EOS contributes to the misconception that density is a variable in the code, which it is not.

[glwagner]

After talking a bit with @ali-ramadhan, I think we might have settled on the following solution:

• Introduce a new type called PhysicalParameters that holds g and ρ0
• Introduce a new abstract type called Rotation (or some such) that encodes information about the background rotation rate of the model --- TangentPlane with f and β, and possibly non-traditional Coriolis parameters, another type for the sphere, etc.
• Group viscous and diffusive transport coefficients into the upcoming TurbulenceClosure type, allowing for isotropic constant transport coefficients, anisotropic constant coefficients, or nonlinear closure schemes.

These changes will also require us to compute buoyancy rather than a density perturbation, and may motivate us to simplify the code by defining variables in pressures as having units of the 'kinematic pressure', which is the ordinary pressure divided by ρ0.

[glwagner]

This is partially solved by #234 and #245.

We still may want to introduce a new type for rotation. I was thinking that Coriolis and a new model field model.coriolis is actually more descriptive and specific than Rotation. Then we can write things like f = coriolis.f, which is utterly clear, and not have to worry about names like fCor.

[glwagner]

My latest thoughts:

We design a rotation functionality that resembles what we do for turbulence closures. The types could look something like

abstract type AbstractRotation end

struct FPlane{T} <: AbstractRotation
    f :: T
end

struct BetaPlane{T} <: AbstractRotation
    f0 :: T
    β :: T
end

struct TiltedFPlane{T} <: AbstractRotation
    fz :: T
    fy :: T
end

These types would then be associated with functions that add the associated terms to the u, v, and w momentum equations; for example:

@inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v)
@inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u)
@inline z_f_cross_U(i, j, k, grid::Grid{T}, rotation::FPlane, U) where T = zero(T)

I propose we add the field model.rotation to model, and put this new functionality in a new file called rotation.jl.