Tracer Flux Boundary Condition with HydrostaticFreeSurfaceModel

[darosaj]

Hi,

I'm trying to run some different variations on Oceananigans solvers, BCs, and other forcing to try and learn the package better.

I'm attempting to simulate the same setup (BC's, forcing, etc) I had done with the NonhydrostaticModel, but with the HydroStaticFreeSurfaceModel.

I'm trying to force my simulation with wind stress at the top boundary in the x-direction, with no-slip boundaries on the bottom immersed boundary field. I also am trying to impose a temperature and salinity flux at the surface.

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx), immersed=ValueBoundaryCondition(0))

T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵀ), immersed = GradientBoundaryCondition(dTdz))

evaporation_bc = FluxBoundaryCondition(Jˢ, field_dependencies=:S, parameters=evaporation_rate)
S_bcs = FieldBoundaryConditions(top=evaporation_bc)

However, when I use the HydrostaticFreeSurfaceModel I get this error:

[ Info: Initializing simulation...
ERROR: InvalidIRError: compiling MethodInstance for Oceananigans.BoundaryConditions.gpu__apply_z_bcs!(::KernelAbstractions.CompilerMetadata{…}, ::OffsetArrays.OffsetArray{…}, ::Tuple{…}, ::ImmersedBoundaryGrid{…}, ::BoundaryCondition{…}, ::BoundaryCondition{…}, ::Tuple{…}) resulted in invalid LLVM IR
Reason: unsupported call to an unknown function (call to ijl_get_nth_field_checked)

I've tried removing / changing the boundary conditions for my simulation, and it seems like I only get this error when trying to set the tracer flux boundary condition at the top surface.

Is it not possible to set the flux BC at the top surface for the tracer with the HydrostaticFreeSurfaceModel, or is there something different I need to do for the free surface solver when trying to impose this boundary condition.

Thank you for all of your help!

All the best,

Josh

[glwagner]

You can definitely set tracer boundary conditions. Can you show us the code for Jᵀ?

[darosaj]

Hi,

Thank you for the help and advice. I was attempting to just use the same parameters as I had been using for the NonHydrostaticModel, to get it up and running, but that was a mistake. I should have adjusted to correctly use the appropriate closure and BCs.

I just ran it without the S_bcs and it ran successfully.

Why does this BC cause more of an issue with the HydrostaticFreeSurfaceSolver?

Thank you again for all of your help. I'm really enjoying learning the package so far! (also I'm excited to read that new paper you had posted in the discussion board as well)

[glwagner]

Im surprised there is an issue with the BC in Nonhydrostatic vs hydrostatic...

If you can paste your whole script (ideally removing unnecessary parts) then I can try to reproduce the issue

[darosaj]

Hi, I definitely can. Cleaning up the script to paste it in now.

In the meantime:

I found if I changed the Salinity BC to:

@inline Jˢ(x, y, t) = - 1e-3 / hour # [salinity unit] m s⁻¹
evaporation_bc = FluxBoundaryCondition(Jˢ)
S_bcs = FieldBoundaryConditions(top=evaporation_bc)

or

evaporation_rate = 1e-3 / hour # m s⁻¹
@inline Jˢ(x, y, t, evaporation_rate) = - evaporation_rate * 35 # [salinity unit] m s⁻¹
evaporation_bc = FluxBoundaryCondition(Jˢ)
S_bcs = FieldBoundaryConditions(top=evaporation_bc)

the script is able to run.

However, when I pass in

field_dependencies=:S

as seen here:

@inline Jˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S # [salinity unit] m s⁻¹

evaporation_rate = 1e-3 / hour # m s⁻¹

evaporation_bc = FluxBoundaryCondition(Jˢ, field_dependencies=:S, parameters=evaporation_rate)
S_bcs = FieldBoundaryConditions(top=evaporation_bc)

I get that same error I had posted above.

[simone-silvestri]

I think this error has to do with the bug that should be corrected in #4008. The problem is I found some inconsistent behavior with other parts of the model I am still trying to solve.
A discrete boundary function written like

@inline Jˢ(i, j, grid, clock, fields, evaporation_rate) = @inbonds - evaporation_rate * fields.S[i, j, grid.Nz] # [salinity unit] m s⁻¹

evaporation_rate = 1e-3 / hour # m s⁻¹

evaporation_bc = FluxBoundaryCondition(Jˢ, discrete_form=true, parameters=evaporation_rate)
S_bcs = FieldBoundaryConditions(top=evaporation_bc)

should do the work.
In the meantime I will try to solve #4008

[darosaj]

Hi Simone, thanks so much for the help.

I'll implement that fix and look to update when it becomes available.

[darosaj]

Hi, Here is the full script that I have errors with when trying to included the Salinity BC as formulated:

using Oceananigans
using Oceananigans.Units
using Oceananigans.Units: minute, minutes, hour
using CairoMakie
using Random
using Printf
using NetCDF

"""
Create the grid / grid fitted bottom boundary condition
"""

bot = 10meters
architecture = GPU() 

underlying_grid = RectilinearGrid(architecture, topology = (Periodic, Periodic, Bounded),
                       size = (100, 100, 50),
                       halo = (4, 4, 4),
                       x = (-50meters, 50meters),
                       y = (-50meters, 50meters),
                       z = (-bot, 0)) 


# Height and width
H = 4meters
W = 5meters

mountain(x, y) = H * exp(-(x^2 + y^2) / 2W^2)
bottom(x,y) = -bot + 1meter + mountain(x,y)
grid_fitted_bot = GridFittedBottom(bottom)
grid = ImmersedBoundaryGrid(underlying_grid, grid_fitted_bot)



"""
Setup the bouyancy equation of state and other boundary conditions:
"""

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = 2e-4,
                                                                    haline_contraction = 8e-4))


Q = 200.0  # W m⁻², surface _heat_ flux
ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean
cᴾ = 3991.0 # J K⁻¹ kg⁻¹, typical heat capacity for seawater

Jᵀ = Q / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux

dTdz = 0.01 # K m⁻¹

T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵀ), immersed = GradientBoundaryCondition(dTdz))


u₁₀ = 10    # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3 # dimensionless drag coefficient
ρₐ = 1.225  # kg m⁻³, average density of air at sea-level

τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

#Immersed boundary condition to create nolip on the immersed bc faces
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx), immersed=ValueBoundaryCondition(0))
v_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0))


"""
Tidal Forcing / Internal forcing input
"""
coriolis = FPlane(f=1e-4)

T₂ = 12.421hours

ω₂ = 2π / T₂ # radians/sec

ϵ = 0.1 # excursion parameter

U_tidal = ϵ * ω₂ * 20kilometers

tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂

@inline tidal_forcing(x, y, z, t, p) = p.tidal_forcing_amplitude * sin(p.ω₂ * t)

u_forcing = Forcing(tidal_forcing, parameters=(; tidal_forcing_amplitude, ω₂))


"""
Salinity Boundary conditions:
"""

@inline Jˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S # [salinity unit] m s⁻¹

evaporation_rate = 1e-3 / hour # m s⁻¹

evaporation_bc = FluxBoundaryCondition(Jˢ, field_dependencies=:S, parameters=evaporation_rate)

S_bcs = FieldBoundaryConditions(top=evaporation_bc)


# Define the model with the immersed boundary grid
model = HydrostaticFreeSurfaceModel(; grid, buoyancy,
                            #free_surface = ExplicitFreeSurface(),
                            tracer_advection = WENO(),
                            momentum_advection = WENO(),
                            tracers = (:T, :S),
                            coriolis = FPlane(f=1e-4),
                            boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs, S=S_bcs),
                            closure = ScalarDiffusivity(ν=1e-6, κ=1e-6),
                            forcing = (; u = u_forcing))

"""
Intitial Conditions:
"""

## Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

## Temperature initial condition: a stable density gradient with random noise superposed.
Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

## Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = 0.1*U_tidal

## `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, T=Tᵢ, S=35)



"""
Setup callback and simulation outputs
"""

simulation = Simulation(model, Δt=10.0, stop_time=40minutes)

# The `TimeStepWizard` helps ensure stable time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 1.0.

wizard = TimeStepWizard(cfl=0.5, max_change=1.1, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# Nice progress messaging is helpful:

## Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n",
                                iteration(sim), prettytime(sim), prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))

add_callback!(simulation, progress_message, IterationInterval(20))

## Create a NamedTuple with free surface
free_surface = (; eta=model.free_surface.η)

filename = "ToyProblemTestFreeSurface2"

outputs = merge(model.velocities, model.tracers)

# Set up the NetCDFOutputWriter
simulation.output_writers[:slices] = NetCDFOutputWriter(model, outputs;
                                                        filename = filename * ".nc",
                                                        schedule = TimeInterval(1minute),
                                                        overwrite_existing = true)

"""
Run simulation:
"""

run!(simulation)

Thanks again for all of the help!