Oceananigans.jl: Sponge layer not behaving as expected

[tomchor]

I set-up a sponge layer for my simulation following the examples here, but I'm getting results that are different from what I expected.

Here's how I'm setting my grid, model and sponge:

#-----
Lx = Ly = 16_000; Lz = 1_000
topology = (Periodic, Periodic, Bounded)
factor = 1
grid = RegularCartesianGrid(size=(2048÷factor, 1, 256÷factor),
                            x=(0, Lx),
                            y=(0, 1),
                            z=(-Lz/2, +Lz/2),
                            topology=topology)
#-----


#-----
ubc = UVelocityBoundaryConditions(grid,
                                  bottom = BoundaryCondition(Gradient, 0),
                                  top = BoundaryCondition(Gradient, 0),
                                 )
vbc = VVelocityBoundaryConditions(grid,
                                  bottom = BoundaryCondition(Gradient, 0),
                                  top = BoundaryCondition(Gradient, 0),
                                 )
wbc = WVelocityBoundaryConditions(grid,
                                  bottom = BoundaryCondition(Value, 0),
                                  top = BoundaryCondition(Value, 0),
                                 )
bbc = TracerBoundaryConditions(grid,
                               bottom = BoundaryCondition(Gradient, N2),
                               top = BoundaryCondition(Gradient, N2),
                               )
#-----

#-----
# Set-up sponge layer
bottom_mask = GaussianMask{:z}(center=minimum(znodes(Face, grid)), width=grid.Lz/10)
mom_sponge = Relaxation(rate=1/10, mask=bottom_mask, target=0)
#----

model = IncompressibleModel(
                            architecture = GPU(),
                            grid = grid,
                            advection = UpwindBiasedThirdOrder(),
                            timestepper = :RungeKutta3,
                            closure = IsotropicDiffusivity(ν=1e-4, κ=1e-4),
                            coriolis = FPlane(f=f0),
                            tracers = (:b,),
                            buoyancy = BuoyancyTracer(),
                            boundary_conditions = (b=bbc, u=ubc, v=vbc, w=wbc),
                            forcing = (u=mom_sponge, v=mom_sponge, w=mom_sponge),

From the docs, I expected this to

relax velocity fields to zero [...] in the bottom 1/10th of the domain

which leads me to believe that there's some a sharp cut-off in the function mask so that velocities above the range [-500, -400]aren't affected.

However, when I plot the vertical velocities (and other variables with zero mean) I get the results below, which indicate that the sponge layer is active in the approximate range [-500, -200], which is around 3 times larger than what I (thought I) set up.

Am I missing something here?

Thanks!

[francispoulin]

Can you plot bottom_mask to make sure that looks like what you expect it to be?

[tomchor]

I honestly don't know how to do that. It's type is GaussianMask and all of the things I tried to make it into something plottable have failed. How can I plot it?

[ali-ramadhan]

Ah not sure we can plot GaussianMask since it's a struct. But yeah,

relax velocity fields to zero [...] in the bottom 1/10th of the domain

might not be the best description of GaussianMask since width is the standard deviation of the Gaussian sponge (so the sponge decays over a lengthscale width).

I guess if you want to reduce the spatial effect of the sponge you could reduce the width.

[francispoulin]

Sorry @tomchor for my suggestion. It might be nice to be able to plot the mask so the user knows what they are actually doing. Do you people agree? If yes I can create an issue asking whether we want to do this or not.

[francispoulin]

From relaxation.jl, I see that GaussianMask in the z direction is defined as follows:

@inline (g::GaussianMask{:z})(x, y, z) = exp(-(z - g.center)^2 / (2 * g.width^2))

I see that your center is center=minimum(znodes(Face, grid)).

Maybe you can define this function with replacing z with the znodes and then plot it to make sure everything looks right?

[ali-ramadhan]

Sorry @tomchor for my suggestion. It might be nice to be able to plot the mask so the user knows what they are actually doing. Do you people agree? If yes I can create an issue asking whether we want to do this or not.

That could be a pretty neat feature (and helpful for debugging). Might be easiest to do by adding a plotting recipe (#1115)?

[glwagner]

This behavior might be expected. Within 100 m of the bottom boundary, the Gaussian has a value

julia> exp(-0.5 * 1^2)
0.6065306597126334

so I think we expect the solution to be strongly damped.

within 200 m, it falls to 14%

julia> exp(-0.5 * 2^2)
0.1353352832366127

Which is still significant. At 300 m its about 1%, so it probably has a negligible impact above -200 m.

[navidcy]

@francispoulin

It might be nice to be able to plot the mask so the user knows what they are actually doing. Do you people agree? If yes I can create an issue asking whether we want to do this or not.

I couldn't agree more. It's amazing how much you can see when you look.

[tomchor]

Thanks for the help guys. I also agree with @francispoulin that plotting the mask would make it very helpful!

Here are some thought that I had that might be worth discussing:

1 - Sponge layers (afaik) should only alter the NS equations in a finite part of the domain, right? Otherwise, even far from the relaxed edge, you're not really solving the NS equations --- it's being modified by an extra unphysical forcing. So one question is: is this really the behavior you want your sponge layer to have? (And here I'm assuming you're using relaxation as a synonym for sponge, because that's what the docs imply.)
2 - The nomenclature is kinda misleading, because width, coupled with the wording used in the docs really makes one believe that the effects of the relaxation function are confined to a layer of width width. For a reference, here's what's written on the docs:

We illustrate usage of mask and target by implementing a sponge layer that relaxes velocity fields to zero and restores temperature to a linear gradient in the bottom 1/10th of the domain

Hopefully that makes sense.

[navidcy]

@tomchor, indeed you are correct on 1. above, but let me remark that sponge layers could be a bit tricky. I explain:

What you want sponge layers to do is to allow flow propagation into the sponge region but where drag would dissipate any fluid motion there. The "tricky" part is to make the sponge layer transition smooth enough so that the fluid does not experience it as a "wall". If the transition region is very short then the fluid "sees" the sponge as a wall and waves are reflected back into the fluid. But how "short" is short-enough depends on the scales of the fluid motion... In my experience always there is some fiddling to be done for every problem or even for same problem at different parameter regimes.

[glwagner]

Relaxation applies to more than just sponge layers --- it's a convenience wrapper for a forcing function that "relaxes" a field to some target at a specified rate, limited to a region outside a mask function. Relaxation is actually a special type of callable ContinuousForcing:

https://github.com/CliMA/Oceananigans.jl/blob/35f749ef09fe021d13e8bad6433d889e8c39d6ca/src/Forcings/relaxation.jl#L81-L82

The snippet shows that mask is callable with the signature x, y, z and target is callable with the signature x, y, z, t. Any appropriate mask and target will work. Also, we've provided a parameterized convenience type for a Gaussian mask.

So, sponge layers are simply one application of Relaxation.

@tomchor, perhaps what you're saying is that a Gaussian is not appropriate for a sponge layer mask, since one might want a function with a sharper cutoff (subject to @navidcy's caveats). I think this is a valid criticism --- it might make sense to add more appropriate parameterized types (using tanh's, etc?) This could be added in a PR that also updates the documentation and provides a better example for implementing a sponge layers.

[navidcy]

I admit I haven't read @tomchor script above in detail before making my comment...

From my experience, sponge layers are usually applied as linear drag with a variable linear coefficient. E.g.,

∂u/∂t = ... + μ(z) u
∂v/∂t = ... + μ(z) v

where μ(z) is a function that is 0 everywhere and becomes 1 in the sponge with some transitional region of certain width. As, @glwagner points out, tanh((z-z_sponge) / width_sponge) is a good starting point.

[tomchor]

@glwagner yes, that was more along the lines of what I was saying. From the docs it seemed to me using a Gaussian mask in the relaxation should be equivalent to a sponge layer, which is why I expected the gaussian function to be set-up with a sharp cut-off. But I think we're on the same page about what a sponge layer should be.

So I think it may be a matter of maybe making the docs clearer? Maybe also, like you said, constructing an example that uses sponge layers (I think there's an issue about there somewhere, right?). The Langmuir example might a good one. In the original paper by Jim they use a radiation boundary condition at the bottom which doesn't have an analog in the Langmuir example at the moment.

[tomchor]

Based on the answers here I tried to set my own mask function. It appears to run fine on the CPU but when I try to run it on the GPU I get this error when running the simulation (with a few subsequent errors after these lines, but this is the first one):

ERROR: LoadError: InvalidIRError: compiling kernel gpu_calculate_Gw!(Cassette.Context{nametype(CUDACtx),KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIteration.StaticSize{(2048, 1, 256)},KernelAbstractions.NDIteration.DynamicCheck,Nothing,Nothing,KernelAbstractions.NDIteration.NDRange{3,KernelAbstractions.NDIteration.StaticSize{(2048, 1, 256)},KernelAbstractions.NDIteration.StaticSize{(1, 256, 1)},Nothing,Nothing}},Nothing,KernelAbstractions.var"##PassType#253",Nothing,Cassette.DisableHooks}, typeof(Oceananigans.TimeSteppers.gpu_calculate_Gw!), OffsetArrays.OffsetArray{Float64,3,CUDA.CuDeviceArray{Float64,3,1}}, RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, UpwindBiasedThirdOrder, FPlane{Float64}, Nothing, AnisotropicDiffusivity{Float64,Float64,Float64,NamedTuple{(:b,),Tuple{Float64}},NamedTuple{(:b,),Tuple{Float64}},NamedTuple{(:b,),Tuple{Float64}}}, NamedTuple{(:velocities, :tracers),Tuple{NamedTuple{(:u, :v, :w),Tuple{Oceananigans.Fields.ZeroField,Oceananigans.Fields.ZeroField,Oceananigans.Fields.ZeroField}},NamedTuple{(:b,),Tuple{Oceananigans.Fields.ZeroField}}}}, NamedTuple{(:u, :v, :w),Tuple{OffsetArrays.OffsetArray{Float64,3,CUDA.CuDeviceArray{Float64,3,1}},OffsetArrays.OffsetArray{Float64,3,CUDA.CuDeviceArray{Float64,3,1}},OffsetArrays.OffsetArray{Float64,3,CUDA.CuDeviceArray{Float64,3,1}}}}, NamedTuple{(:b,),Tuple{OffsetArrays.OffsetArray{Float64,3,CUDA.CuDeviceArray{Float64,3,1}}}}, Nothing, NamedTuple{(:u, :v, :w, :b),Tuple{typeof(Oceananigans.Forcings.zeroforcing),typeof(Oceananigans.Forcings.zeroforcing),Oceananigans.Forcings.ContinuousForcing{Cell,Cell,Face,Nothing,Relaxation{Float64,typeof(bottom_mask_func),Int64},Nothing,Tuple{Int64},Tuple{typeof(identity)}},typeof(Oceananigans.Forcings.zeroforcing)}}, NamedTuple{(:time, :iteration, :stage),Tuple{Float64,Int64,Int64}}) resulted in invalid LLVM IR
Reason: unsupported dynamic function invocation (call to overdub)
Stacktrace:
 [1] bottom_mask_func at /glade/work/tomasc/ocnn/ocnn_JD15/jd15_00.jl:56
 [2] Relaxation at /glade/u/home/tomasc/.julia/packages/Oceananigans/RxUYW/src/Forcings/relaxation.jl:84
 [3] ContinuousForcing at /glade/u/home/tomasc/.julia/packages/Oceananigans/RxUYW/src/Forcings/continuous_forcing.jl:113
 [4] w_velocity_tendency at /glade/u/home/tomasc/.julia/packages/Oceananigans/RxUYW/src/TimeSteppers/velocity_and_tracer_tendencies.jl:167
 [5] macro expansion at /glade/u/home/tomasc/.julia/packages/Oceananigans/RxUYW/src/TimeSteppers/calculate_tendencies.jl:167
 [6] gpu_calculate_Gw! at /glade/u/home/tomasc/.julia/packages/KernelAbstractions/jAutM/src/macros.jl:80
 [7] overdub at /glade/u/home/tomasc/.julia/packages/Cassette/158rp/src/overdub.jl:0
Reason: unsupported dynamic function invocation (call to overdub)

[1] bottom_mask_func at /glade/work/tomasc/ocnn/ocnn_JD15/jd15_00.jl:56 refers to the return line in the following definition that I'm using for the mask:

sponge_one = minimum(oc.Grids.znodes(Face, grid)); sponge_zero = sponge_one + grid.Lz/10
heaviside(x::AbstractFloat) = ifelse(x < 0, zero(x), one(x))
function bottom_mask_func(x, y, z)
    sponge_one = -grid.Lz/2; sponge_zero = sponge_one + grid.Lz/10
    return heaviside(-(z-sponge_zero)) * (z - sponge_zero)^2 / (sponge_one - sponge_zero)^2
end
mom_sponge = Relaxation(rate=1/20, mask=bottom_mask_func, target=0)

Am I missing something here? Like I said, it runs fine on the CPU.

Thanks!

[ali-ramadhan]

When running GPU kernels, KernelAbstractions.jl uses Casette.jl (which defines overdub) to substitute GPU/CUDA functions for CPU functions wherever necessary, e.g. CUDA.sin(1) instead of sin(1). So usually when I see an error like

Reason: unsupported dynamic function invocation (call to overdub)

I assume there was a a function it failed to convert to be GPU-friendly. But the only function I see in there is your heaviside which should be fine.

@tomchor Can you post or share your full script? Might help with debugging.

[tomchor]

@ali-ramadhan Thanks for the info. Here's the my code: https://pastebin.com/ACKJhEj7
Still haven't been able to figure it out...

Thanks!

[ali-ramadhan]

Thanks @tomchor!

I think the issue was that the GPU was not able to access the grid inside bottom_mask_func, but obviously the error did not give any hints at all.

In general, for functions that are to be embedded within GPU kernel (like your bottom_mask_func) I would not use anything but the input arguments and const values in the function.

So with this change I was able to run your script on the GPU:

# Set-up sponge layer
heaviside(x) = ifelse(x < 0, zero(x), one(x))

const H = grid.Lz

sponge_one = minimum(oc.Grids.znodes(Face, grid))
sponge_zero = sponge_one + grid.Lz/10

function bottom_mask_func(x, y, z)
    sponge_one = -H/2
    sponge_zero = sponge_one + H/10
    return heaviside(-(z-sponge_zero)) * (z - sponge_zero)^2 / (sponge_one - sponge_zero)^2
end
#bottom_mask = GaussianMask{:z}(center=-grid.Lz/2, width=grid.Lz/40)
mom_sponge = Relaxation(rate=1/20, mask=bottom_mask_func, target=0)

You could also declare sponge_one and sponge_zero as const values outside bottom_mask_func if they won't change, and to save a few measly FLOPS.

[tomchor]

Thanks! Indeed that seems to work :+1:

Apparently running on GPUs has a few caveats even after you guys did most of the hard work in Oceananigans

[ali-ramadhan]

Awesome! Yeah there are some GPU idiosyncrasies (biggest one being that GPU errors can be very obscure) with user-defined forcing functions and boundary conditions.

Should this issue be closed or do you want to modify the title to reflect that it's about improving the language in the documentation?

[glwagner]

Might be good to clarify some of this discussion: a "sponge layer" refers to any region in a simulation domain that dissipates energy. A wide variety of forcing functions are "valid" sponge layers, including the example in the documentation.

The Gaussian mask appears to be behaving correctly, and as one expects given the documentation, source code, and known/measurable/plottable properties of Gaussians.

Using a Gaussian to define the masked region is valid strategy used successfully in research. Of course, this does not mean that other masking functions may be more appropriate, or that a "smooth step" function might be more useful to users, since the width of the sponge, and width of the sponge layer transition region may be modulated independently. This is not the case for the Gaussian, which has only one width parameter.

A sharp heaviside is not a good choice for a sponge layer when the sponge is intended to absorb radiating internal waves. The reason is that waves can reflect off a sponge layer if the transition between undamped and damping regions is too abrupt (as pointed out by @navidcy). was under the impression this might even be an issue for problems without waves... ?

The requirement that sponge layers involve smooth transitions does mean that sponge layers can take up a substantial portion of the computational domain. This is often cited as a downside of sponge layers (eg Klemp and Durran 1993).

[glwagner]

From the docs it seemed to me using a Gaussian mask in the relaxation should be equivalent to a sponge layer, which is why I expected the gaussian function to be set-up with a sharp cut-off.

@tomchor can you clarify what you mean when you say you expected the "gaussian function to be set-up with a sharp cut-off"? Do you mean you expected the width of the Gaussian to be narrower?

We can change the example in the docs to have a narrower width; perhaps Lz/20 would be appropriate. If the width is too narrow, the sponge will not adequately damp radiating waves, so there is a balance to be achieved there.

From the docs it seemed to me using a Gaussian mask in the relaxation should be equivalent to a sponge layer

To clarify: the implementation in the documentation is a valid implementation of a sponge layer, correct?

[glwagner]

Apparently running on GPUs has a few caveats even after you guys did most of the hard work in Oceananigans

@tomchor the reason your code did not compile is because this function

function bottom_mask_func(x, y, z)
    sponge_one = -grid.Lz/2; sponge_zero = sponge_one + grid.Lz/10
    return heaviside(-(z-sponge_zero)) * (z - sponge_zero)^2 / (sponge_one - sponge_zero)^2
end

references grid.Lz as a global variable. Global variables can only be used inside functions on the GPU if they are declared const as @ali-ramadhan has done.

One strategy to avoid having to use const is to build callable objects like GaussianMask that store their parameters. This is why a parameterized GaussianMask callable object is provided for users to use as mask.

It might be nice to add a SigmoidMask masking function as well that uses tanh and has a width, center, and direction.

[tomchor]

@glwagner Thanks for the clarifying the global-var-limitation for GPUs! This mat be something to add to the docs at some point (I didn't see it there!)

(I'll answer the sponge layer questions on a separate answer in a minute!)

[tomchor]

@tomchor can you clarify what you mean when you say you expected the "gaussian function to be set-up with a sharp cut-off"? Do you mean you expected the width of the Gaussian to be narrower?

No, I mean I expected the sponge layer to be confined to a finite part of the domain. That is, I expected the sponge layer to be _identically_ zero for most of the domain.

Us not "being on the same page" here (for the lack of a better description! haha) may be due to different philosophies of what a sponge layer should be. So, due to my background, a sponge layer should be a function that _only_ acts in a _finite_ part of the domain. So, for example, in my code I set-up the sponge layer mask as

mask = H(-(z-z₀)) ⋅ (z - z₀)² / (z₁ - z₀)²

with H(z) being the Heaviside function. Comparing this to the Gaussian mask, this has a disadvantage of not being as smooth (although it should be smooth enough because of the (z - z₀)² / (z₁ - z₀)² term, and we can always increase the order of the exponents). However it has a very important advantage: I know _exactly_ when it stops influencing my solution. So I know that in the range (z₀, max(z)] my solver is solving the NS equations (plus or minus some approximations) without any added nonphysical term.

I was always taught (and I still agree) that setting up a sponge layer mask as something that only reaches zero at infinity (like a Gaussian) is not recommended since you don't know when do start believing your solution. In practice if you're 4σs away from the center of the Gaussian you're pretty sure that the sponge layer influences are minimal (but they're not exactly zero!). But a harder question is: how far away is far enough?

A note here should be that if my mask was simply mask = H(-(z-z₀)) it would probably fail as a sponge layer because, as @navidcy and @glwagner have pointed out, it would reflect a ton of waves. So the transition should be done with care.

So, again, I think it boils down to your background and modeling philosophy! But it would be nice to have options for finite-range masks, such as the one I imposed, as well. (Although there are many other functions that are used and mine is probably the easiest one you can come up with...)

I hope I made the issue clearer! And sorry about any confusion.