KE budget with open boundary condition

[peitianyi23-sudo]

Hi Oceananiagns users,

I was playing with the OpenBoundaryCondition in a 2D box. I considered an inflow from a orifice at the bottom and a outflow at the top, both using a Gauss-like profile (To be honset I am not sure whether the outflow BC is correct or not. I was thinking that maybe I need to keep the mass conservation).

I am having trouble to close the KE budget using Oceanostic. I haven't studied Oceanostics in depth yet. But if i am getting this right, the KE could be written as: $u_j \partial_t u_j + u_j \partial_i(u_i u_j) = -u_j \partial_j P + u_j \partial_i \tau_{ij}$
, for my system with no body force. Thus in the steady state, I should have the intergal outputs satisfying (my code attached below):
$∫fnu - ∫fp - ∫fu \approx 0$. However, the result is wrong for now (pls see the fig).

I am guessing I make some mistake with the quantities computed from Oceanostics (sign error??). I would greatly appreciate it if someone could point out! :)

Btw, Oceanostic can also compute a KineticEnergyTendency $u\cdot G$. Is $\int u\cdot G = \int u_j \partial_i \tau_{ij} - \int u_j \partial_i(u_i u_j)$ ?

using Oceananigans
using CUDA
using NCDatasets
using Oceanostics
using CairoMakie

##
function getOrifice(Nx,orifice)
    x_edges = range(-0.5, 0.5; length = Nx + 1)
    x = 0.5 .* (x_edges[1:end-1] .+ x_edges[2:end])
    y = zeros(length(x))
    y = exp.(-((x ./(orifice/2)).^2))
    return x, y
end

Nx = 256
Nz = 256
Lx = Lz = 1.0
Sz = 1.5
orifice = 0.2
tanh_z_faces(k) = Lz * (1 + tanh(Sz * (2 * (k - 1) / Nz - 1))/tanh(Sz)) / 2 - 0.5

x0, jets = getOrifice(Nx,orifice)
const Re = 1000.0
const nu = 1.0 / Re
const deltaT = 1/Nx*0.5

filename = "JetSingle2D_Re$(round(Int, Re))_$(Nx)x$(Nz)" 

grid = RectilinearGrid(GPU();
    size = (Nx, Nz),
    x = (-Lx/2, Lx/2),
    z = tanh_z_faces,
    topology = (Bounded, Flat, Bounded)
)

no_slip_bc = ValueBoundaryCondition(0)
no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc)

jets_gpu      = CuArray(jets)
jets_bcs = FieldBoundaryConditions(no_slip_bc,
           top=OpenBoundaryCondition(jets_gpu),
           bottom=OpenBoundaryCondition(jets_gpu))

closure = ScalarDiffusivity(VerticallyImplicitTimeDiscretization(),ν=nu) 

model = NonhydrostaticModel(; grid,
                            boundary_conditions=(u=no_slip_field_bcs, v=no_slip_field_bcs, w=jets_bcs),
                            closure)

                            
simulation = Simulation(model, Δt=deltaT, stop_time=100)
wizard = TimeStepWizard(cfl=0.5, diffusive_cfl=0.5)   
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

progress = ProgressMessengers.BasicMessenger()
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

u, v, w = model.velocities
ω = @at (Center, Center, Center) ∂x(w) - ∂z(u)
∫KE = Integral(KineticEnergyEquation.KineticEnergy(model))
∫ε = Integral(KineticEnergyEquation.DissipationRate(model))
∫uG = Integral(KineticEnergyEquation.KineticEnergyTendency(model))
∫fp = Integral(KineticEnergyEquation.PressureRedistribution(model))
∫fu = Integral(KineticEnergyEquation.Advection(model))
∫εᴰ = Integral(KineticEnergyEquation.Stress(model))

output_fields = (; ω, ∫KE, ∫ε, ∫εᴰ, ∫fp, ∫uG, ∫fu)
checkpoint_dir = joinpath(@__DIR__, filename)
isdir(checkpoint_dir) || mkpath(checkpoint_dir)
checkpoint_prefix = joinpath(checkpoint_dir, filename)

simulation.output_writers[:checkpointer] =
    Checkpointer(model, schedule=TimeInterval(50), prefix=checkpoint_prefix)

simulation.output_writers[:nc] = NetCDFWriter(model, output_fields,
                                              filename = joinpath(checkpoint_dir, filename),
                                              schedule = TimeInterval(0.5),
                                              overwrite_existing = true)

run!(simulation)

[peitianyi23-sudo]

Oh, one additional thing is that here I use an non-uniform grid in z direction, with the VerticallyImplicitTimeDiscretization. When setting up the NonhydrostaticModel, the code shows that "Implicit-explicit time-stepping with RungeKutta3TimeStepper is not tested". I'm not sure whether this is a big issue? The code runs much slower without the VerticallyImplicitTimeDiscretization method.

[glwagner]

The code runs much slower without the VerticallyImplicitTimeDiscretization method

Does it execute each time step slower, or do you just need to use a shorter time step?

[peitianyi23-sudo]

hmm, thats weird. I remeber last week I used a 1536^2 grid with the same tanh_z_face spacing (similar problem but more inflows & outflows), and the code would be stucked after "... initial time step complete" if i use an explicit form. But just now I test the single inflow case presented above, fixing the same tanh spacing, using two grids 256^2 and 1024^2, with or without VerticallyImplicitTimeDiscretization, to run 1000 iterations on my 4060ti pc.

256^2 explicit: 13.028 seconds; 256^2 implicit: 23.057 seconds.
1024^2 explicit: 22.220 secones; 1024^2 implicit: 58.955 seconds.

So as for speed per iteration, actually explicit form runs faster. I will check the old code to fix it or repeat the stucking issue.

[glwagner]

I would expect explicit to run faster. The reason is that the implicit discretization does a tridiagonal solve 3x per time-step, requiring a 2D kernel with a reduction-like algorithm that does 3 loops over the grid in z. The explicit discretization doesn't require a solver for taking a time-step (just for pressure).

[glwagner]

I think this is really a question for Oceanostics. If you are interested in the entire budget though, you may be doing more work than needed? You only need to compute the dot product of $u$ and the tendency $G_u$ to obtain the residual of the kinetic energy.

I'm confused that the pressure and advection terms don't integrate to zero. Are they formulated in divergence form and at the current staggered grid location?

I've commented on this before and have been reassured that it's not a big deal, but I'm still confused here about whether or not this method of computing the KE budget respects the time discretization. For example, you're using vertically-implicit time discretization for the diffusion term. This has implications for what the term looks like when you formulate the discrete-in-time KE budget. Also does the Dissipation diagnostic do the right thing for implicit time discretization? For this to make sense, I think you would have to re-form the closure with an explicit time discretization and compute a "faux explicit" contribution to the momentum tendency, and use that as the basis for the dissipation term.

Just to illustrate what I mean, consider a discrete time equation for a forward Euler time-stepping method:

$$ \frac{u^{n+1} - u^n}{\Delta t} = \partial_z \left ( \nu \ \partial_z u^{n+1} \right ) + \ G^n_u $$

To find an equation for the variance $\left ( u^n \right )^2$ we need to multiply by $\left ( u^{n+1} + u^n \right ) / 2$ to obtain

$$ \frac{\left ( u^{n+1} \right )^2 - \left ( u^n \right )^2}{2 \Delta t} = \frac{1}{2} \left ( u^{n+1} + u^n \right ) \partial_z \left ( \nu \ \partial_z u^{n+1} \right ) + \frac{1}{2} \left ( u^{n+1} + u^n \right ) \ G_u $$

I think Oceanostics will make an approximation $u^{n+1} = u^n + O(\Delta t)$ --- you have to make this both for the time-interpolated velocity, $u^n + u^{n+1} \approx 2 u^n$, as well as for the velocity field that appears in the diffusion term, to get

$$ \frac{\left ( u^{n+1} \right )^2 - \left ( u^n \right )^2}{2 \Delta t} \approx u^n \partial_z \left ( \nu \ \partial_z u^n \right ) + u^n \ G_u $$

perhaps?

[peitianyi23-sudo]

The figure below shows the pressure field at the final time for the 256^2 explicit case. There is large pressure variation around the outflow as i plot $p|_{z=zmax}$ against x node center. I am not so sure what a domain scale gradient looks like?

I also add a local dissipation rate output to get the local kolmogorov scale. Still, large dissipation exists around the outflow, but the local kolmogorov scale is indeed larger than the grid spacing everywhere.

I will dive into the Oceanostics source code then. Thanks for your note about the boundary conditions. I will modify it using Field.

PerturbationAdvection scheme is really a new thing to me, but it looks cool! I will study the related documentation, but can you suggest a outflow&inflow_timescale of the top and bottom BC for me to play with?

[peitianyi23-sudo]

Still kind of stuck on this budget problem. I turn to the flow past cylinder case provided in https://github.com/CliMA/Oceananigans.jl/blob/main/validation/open_boundaries/cylinder.jl. I run this test with Re=100 and Nx=Ny=128 (everything remains the same).

Similarly I add:

    ∫uG = Integral(KineticEnergyEquation.KineticEnergyTendency(model))
    ∫fp = Integral(KineticEnergyEquation.PressureRedistribution(model))
    ∫fu = Integral(KineticEnergyEquation.Advection(model))
    ∫εᴰ = Integral(KineticEnergyEquation.Stress(model))
    ∫KE = Integral(KineticEnergyEquation.KineticEnergy(model))
    output_intergal = (; ∫KE, ∫εᴰ, ∫fp, ∫uG, ∫fu) 

and seems i get the same issue. I still cant figure out what i miss but I am really curious.

I study the PerturbationAdection 's document and how you guys develop it https://github.com/CliMA/Oceananigans.jl/pull/4687. I think indeed i need to use this scheme, maybe with inflow_timescale =0 and outflow_timescale = inf for both top and bottom bc. It seems that by setting outflow_timescale = inf, the sharp pressure variation and local dissipation around the top outflow orifice becomes better, and the pressure transport $\int u_i\partial_i p$ decreases to around zero. I am still testing, and thanks for telling me that :)

[glwagner]

Can you reduce your example until the budget is closed? This way you can figure out what aspect of the setup is causing the budget to be unclosed, which will in turn help you pinpoint where the issue with the diagnostics lies.

[peitianyi23-sudo]

By "reduce", you mean reduce the grid sizes and Re for my inflow&outflow orifice case?

[glwagner]

No, I mean reduce the complexity such as replacing boundary conditions with defaults and removing the closure.