Available Potential Energy

[francispoulin]

Do we want to be able to compute Available Potential Energy (APE) to get an idea as to the energtics of a flow?

[ali-ramadhan]

Just pasting this definition (approximation?) of APE from Vallis' textbook (2nd edition):

Should be possible to compute/diagnose with a KernelComputedField (#1293).

Also found this paper on oceanic convective available potential energy (CAPE) but it looks pretty complicated and not sure if it's immediately useful: https://doi.org/10.1175/JPO-D-14-0155.1

cc @sandreza since APE might be a useful diagnostic for your mesoscale simulations?

[jklymak]

There was a recent JPO paper on APE, I want to say by Timor Radko, but not 100% sure

[navidcy]

The definition in Winters et al 2006 is very different (at least at first glance)
https://doi.org/10.1017/S002211209500125X

[jklymak]

Maybe: https://www.annualreviews.org/doi/abs/10.1146/annurev-fluid-011212-140620?casa_token=8OILtxnTBfEAAAAA:KI9M5PyS38-6W0vhjUG_4iJ5g0FgHwXc7UwdljGg-gMncp-QOuws3gyNg3qhPsK7lSEMee9r7wt4

[francispoulin]

Thanks @jklymak .

I looked at this paper and the APE for the incompressible Boussinesq model is defined explicitly in section 2.4. The APE is defined in equation 3.6 as

APE     =     \int_V E_a dV 
        =     \int_V \rho g (z - z_R) dV,                                                               (1)

where E_a is the APE density. This is exact and to compute this one needs to find the reference depth, z_R, which is I presume where the sorting comes in.

Then, it says that one can show that the APE density can be written as

E_a       =     - \int_0^\zeta g \tilde \zeta \rho_R' (z - \tilde \zeta) d \tilde \zeta               (2)
   
E_a   \approx   \rho_0 N_r^2 \zeta^2/2                                                                (3)

where N_R^2 = - g/\rho_0 \rho_R' is the Brunt-Vaisala frequency based on the reference density.

Equations (1) and (2) are exact and (3) is approximate, but I believe (3) is exact in the case of constant stratifiation.

I wonder whether using equation (1) is the simplest?

(If anyone can tell me how to make equations in markdown I am happy to change this to make it prettier.)

[glwagner]

I think it's probably a good idea to (cautiously) add diagnostics that are important yet difficult to write, such as this one. Note that this would be the first!

As a historical note, we have resisted adding diagnostics to the code base so far because we believed that a more generic solution would make it "easy" for users to define their own diagnostics (eg, as simple as writing the mathematical expression), and because long lists of diagnostics impose a maintenance burden.

But I think there are some examples (such as APE) for which we probably can't define simply via abstractions like those provided by AbstractOperations. And as the code becomes more stable (and we have more contributors), maintaining a list of difficult-to-code but commonly-used diagnostics is more feasible.

[francispoulin]

I agree that most diagnostics like total tracer and total kinetic energy are pretty easy to compute and if it is a one liner, then we probably don't need to do anything special.

Computing the APE could be a oneliner as well, after we compute the reference depth and density. So maybe that's something we could add somewhere?

I have not done this before but I know it's pretty standard. I presume the MITgcm does and we could maybe borrow from their weather of knowledge?

[jklymak]

I'm not aware of an APE calculation in the MITgcm. Its actually pretty rare because it is quite expensive - you need to do a global sort on the density field!

[francispoulin]

Ah, thanks @jklymak for pointing that out.

Kevin Lamb has a paper that uses a slightly different way of computing APE that looks like

APE = \int_V (\rho - \bar \rho_L) g z dV

where \rho_L is the adiabatically arranged density field.

His paper also has some links to other works that do sorting to compute \rho_L.

[tomchor]

As a historical note, we have resisted adding diagnostics to the code base so far because we believed that a more generic solution would make it "easy" for users to define their own diagnostics (eg, as simple as writing the mathematical expression), and because long lists of diagnostics impose a maintenance burden.

I agree that some diagnostics would be important, but I also agree with the above philosophy of Oceananigans.

In cases like this wouldn't it be good to consider adding a companion repo to oceananigas that has all these diagnostics rather than putting them in Oceananigans by default? Much like LESbrary.jl does, it would be someting that (hopefully) can use KernelComputedFields, but it wouldn't be dedicated only to LES (so I guess it can't be LESbrary.jl).

[navidcy]

Adiabatically arranged to me reads as a global resorting. It’s a tuff computation but it is what it is. What can we do?

…

On 13 Jan 2021, at 04:53, Francis J. Poulin ***@***.***> wrote:  Ah, thanks @jklymak for pointing that out. Kevin Lamb has a paper that uses a slightly different way of computing APE that looks like APE = \int_V (\rho - \bar \rho_L) g z dV where \rho_L is the adiabatically arranged density field. His paper also has some links to other works that do sorting to compute \rho_L. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[francispoulin]

From what I understand, we need to go through the entire domain and find out how much volume (or area in 2D) of each density we have. Then find the basic state where we have the different densities layed out flat, with the heaviest on the bottom. It sounds like the hard part is looking through the density profile to see how much of each density we have.

I have read, and believe, that it gets more complicated if you have geostrophic balance, which we typically would, but starting off with a 2D non-rotating example would be a reasonable first step. Is this what other people understand?

During non-covid times, I could walk down two doors and ask Kevin Lamb for help on how to do this. Even though I can't ask him in person I can email him questions, when they arise.

[glwagner]

In cases like this wouldn't it be good to consider adding a companion repo to oceananigas that has all these diagnostics rather than putting them in Oceananigans by default? Much like LESbrary.jl does, it would be someting that (hopefully) can use KernelComputedFields, but it wouldn't be dedicated only to LES (so I guess it can't be LESbrary.jl).

I think this is a nice approach. I think it could also encourage more community contributions.

[jklymak]

I have read, and believe, that it gets more complicated if you have geostrophic balance, which we typically would, but starting off with a 2D non-rotating example would be a reasonable first step. Is this what other people understand?

I think APE can still be unambiguously defined, even in the case of a geostrophic balance. If you want to decompose your flows in various ways to identify energy reservoirs associated with different classes of motion, that is even more specialized. I'm sure you can do it, but I expect the global sort is still necessary.

This all said, I think a reasonable approximation could be made to calculate the APE on a much coarser timestep to the full physics of the model. If you put such an expensive piece of code in, you would probably want a knob to twiddle to allow it to be coarsened, and then some tuning would be necessary.

[francispoulin]

I agree with you @jklymak , you can certainly compute APE in a rotating fluid, which will require a global sort. I wonder if the issue is how to pick the angle of the pycnoclines in the state of least potential energy? I am sure there are ways of doing these, and we can learn how when needed.

It would be good to give the user an option as to whether they want to compute it and how often. Given what is set up in the time stepping wizard, I don't imagine this will be a concern, and would be nice to give the user the power to control this.

[glwagner]

This all said, I think a reasonable approximation could be made to calculate the APE on a much coarser timestep to the full physics of the model. If you put such an expensive piece of code in, you would probably want a knob to twiddle to allow it to be coarsened, and then some tuning would be necessary.

Not sure I understand the comment but I believe what's being recommended is currently a feature of Oceananigans. Output / diagnostics can be "scheduled" in various ways, either on time intervals, iteration intervals, wall time intervals, or custom criteria designed by the user.

[tomchor]

I'm sure you can do it, but I expect the global sort is still necessary.

Every source that I've read so far that mentions how to compute it says sorting needs to be done. Everything else seems to be an approximation. The hard part is, of course, to compute the "unavailable potential energy", which is the minimum PE of the system. The canonical way to do it apparently is to flatten the array, sort it, and then put it back into the original shape. In Python this can be easily done with (using b instead of ρ)

PE_min = np.sort(b.data.flatten()).reshape(*b.shape)

(disclaimer: I haven't tested if the result is exactly correct, but it should be along those lines.)

I will have to program this diagnostic very soon and I'm thinking of doing it using KernelComputedFields (and then putting it up on a repo). My question is: what's the sorting/reshaping method that works for GPUs?

[francispoulin]

Great to know that python can do this in one line! Maybe it will turn out to be as simple in Julia for both CPUs and GPUs.

[ali-ramadhan]

Reshaping a CuArray should work just like reshaping a regular Array but sorting on GPUs is pretty non-trivial since it has to be made parallel across thousands of little threads.

Looks like CUDA.jl has a quicksort implementation as of 9 days ago though! JuliaGPU/CUDA.jl#431

[francispoulin]

Luckily we are not having this discussion two weeks ago!

[glwagner]

@tomchor if KernelComputedField doesn't work (I suspect it won't, since this problem may require more than one kernel), remember that you can also define your own field that subtypes AbstractField and defines an appropriate compute! method.

[hdrake]

Great to know that python can do this in one line! Maybe it will turn out to be as simple in Julia for both CPUs and GPUs.

Just wanted to point out that APE calculations are only this simple with simple for Cartesian grids with a flat bottom and a linear equation of state. Otherwise, the hard part is actually finding the depths/pressure associated with the sorted profile density profile (see Huang 2005). I've been working on implementing this in an MITgcm simulation with a linear equation of state but complicated topography, and it has been a mess (and very slow).

[francispoulin]

Thanks @hdrake for sharing both the reference and your insights. I guess Oceananigans does not have topography yet but it would be nice to dl thi sfor different equations of state and stratifications.

Please let us know how it goes an dmaybe when you have something you can share we can try and do something similar here?

[glwagner]

@hdrake can you clarify specifically what is difficult about the calculation? If we know the specific calculations we need to make we might be able to support / facilitate them with abstractions.

For example, it may be possible to sort abstract operations --- like the product of a field and a grid metric like cell area or volume, or another multiplicative factor that represents the effect of bathymetry on cell interface area / volume.

[hdrake]

@glwagner I think you're probably right that this can still be done, but the real problem has to do with the fact that density sorting is not unique due to non-linear pressure effects, and so you have to do some kind of iteration with different reference pressures. Maybe that just means an outer loop around the sorting of abstract operations though.

[glwagner]

Ah, I understand. No abstraction can help with that.

This sounds similar to the problems that motivate the formulation of neutral density?

[tomchor]

Thanks, @hdrake, for bringing up this important point. I guess we should also remember that the sorting procedure is only an approximation to the actual calculation of background potential energy. It's super useful and easy in most cases, so most people jump straight to that, but the more "proper" calculation of BPE involves calculating the quantity z_*:

from Winters et alia, (1995).

I think calculating BPE this way still has some complications, but at least I believe the calculation is unique. It might be the preferred for cases with topography, stretched grids, etc.

[tomchor]

Also I'll post here my implementation of the sorting just in case it's useful for some people.

    function flattenedsort(A, dim_order::Union{Tuple, AbstractVector})
        return reshape(sort(Array(permutedims(A, dim_order)[:])), (grid.Nx, grid.Ny, grid.Nz))
    end

    function sort_b(model; average=false)
        b = model.tracers.b
        sorted_B = flattenedsort(interior(b), [3,2,1])
        if !average
            return sorted_B
        else
            return dropdims(mean(sorted_B, dims=(1,2)), dims=(1,2))
        end
    end
    mean_sort_b = (mod)->sort_b(mod; average=true)

I bet this can be done more elegantly but here it is. Thanks to @glwagner for pointing me in the right direction for this.

Note: Right now I have to explicitly transform b into an Array before sorting because the version of CUDA that Oceananigans currently uses doesn't have sort implemented. After we upgrade to Julia 1.6 we'll be able to use the most recent CUDA version and that transformation won't be necessary anymore (I think).

[glwagner]

@hdrake's point references Huang (2005), who point out that a reference state (associated with z* in Winters et al 1995) may not be easily or uniquely definable for seawater with a nonlinear equation of state that depends on salinity, temperature, and pressure:

In practical terms I think the issue here is determine the "adiabatic rearrangement" mentioned by Winters et al. 2005. For a nonlinear equation of state like TEOS-10 where the buoyancy field b in the code depends on depth / hydrostatic pressure, the adiabatic rearrangement of the water column may not be obtained simply by sorting. However, @hdrake points out that an iterative procedure in which parcels are sorted, their buoyancy recalculated at the new depth, and then sorted again, may eventually converge to an adiabatic reference state.

[tomchor]

@glwagner thanks for clarifying. Coming from an atmospheric sciences background, I'm a bit unclear on the details of TEOS-10 and the nuances associated with it. In particular it seems counter-intuitive to me how an equation of state that depends on depth fits in a Boussinesq fluid (where adiabatic rearrangement of parcels shouldn't change their buoyancy). But this probably isn't the best place to have a long discussion about it.

Using the more fundamental definition of z_* does circumvent some of the issues mentioned by @hdrake, no? (Like topography.)

[glwagner]

I think the difficult part is finding z*; this is what Huang 2005 seems to be discussing.

I'm hazy on the details (maybe @hdrake can chime in) but I think the point is that while an adiabatic rearrangement into a stable density profile may exist, it may not be obtainable by a single sorting procedure due to the pressure dependence of the equation of state. Some of the subtleties of the Boussinesq approximation for seawater are discussed in http://pordlabs.ucsd.edu/wryoung/reprintPDFs/SeawaterBoussinesq.pdf

[hdrake]

@tomchor

Using the more fundamental definition of z_* does circumvent some of the issues mentioned by @hdrake, no? (Like topography.)

Somewhat, but even with the "fundamental" definition below there are some implicit assumptions. Winters et al. assume the area A is constant with depth. It's not clear to me that this formula still holds if A=A(z). Also not clear how to calculate a continuous A(z) from a discrete grid in a consistent way, especially with unstructured grids / or partial cells.

[glwagner]

I'm converting this to a discussion!