Reformulating continuous forcing in discrete form giving unexpected results #5369
Description
I am attempting to implement a forcing over a region of the temperature field of a hydrostatic free surface model on a rectilinear grid that restores toward a temperature profile with uniform vertical stratification. I initially implemented this using the Relaxation forcing helper, for example:
using Oceananigans
using Oceananigans.Units
const relaxation_time = 5day
@inline mask(x, y, z) = y < 200kilometers
@inline target(x, y, z, t) = 10.0 + z * 1e-3
forcing = (; T=Relaxation(1 / relaxation_time, mask, target))However, I found when trying to shift from running the simulation on CPU to GPU this caused a InvalidIRError (with Reason: unsupported call to an unknown function (call to ijl_get_nth_field_checked) equivalent to that described in #4165). Based on a suggestion in #4659 for an analogous error when using a continuous form boundary condition on the GPU, I have attempted to reformulate the forcing in discrete form:
@inline function continuous_forcing(x, y, z, t, T)
mask(x, y, z) * (target(x, y, z, t) - T) / relaxation_time
end
@inline function discrete_forcing(i, j, k, grid, clock, model_fields)
x, y, z = node(i, j, k, grid, Center(), Center(), Center())
t = clock.time
T = @inbounds model_fields.T[i, j, k]
continuous_forcing(x, y, z, t, T)
end
forcing = (; T=Forcing(discrete_forcing; discrete_form=true))However using this formulation gives different simulation behaviour from both using the Relaxation implementation and a manually implemented continuous forcing Forcing(continuous_forcing; field_dependencies=:T)). My attempted reformulation in to discrete form is based on reading the forcing documentation page and looking through the source code for the implementations of the forcings, but I am not sure if I am doing something wrong here?
Complete minimal reproducing example
using Oceananigans
using Oceananigans.Units
using Oceananigans.Grids
using Printf
const relaxation_time = 5day
@inline mask(x, y, z) = y < 200kilometers
@inline target(x, y, z, t) = 10.0 + z * 1e-3
@inline function continuous_forcing(x, y, z, t, T)
mask(x, y, z) * (target(x, y, z, t) - T) / relaxation_time
end
@inline function discrete_forcing(i, j, k, grid, clock, model_fields)
x, y, z = node(i, j, k, grid, Center(), Center(), Center())
t = clock.time
T = @inbounds model_fields.T[i, j, k]
continuous_forcing(x, y, z, t, T)
end
function message(sim)
@info @sprintf(
"Time: %s, max(|T|): %.2e", prettytime(sim), maximum(abs, sim.model.tracers.T),
)
end
grid = RectilinearGrid(
CPU();
size=(50, 100, 20),
x=(0, 1000kilometer),
y=(0, 2000kilometer),
z=(-2kilometer, 0),
halo=(6, 6, 3),
topology=(Bounded, Bounded, Bounded),
)
T_forcings = [
"Relaxation" => Relaxation(1 / relaxation_time, mask, target),
"Continuous form" => Forcing(continuous_forcing; field_dependencies=:T),
"Discrete form" => Forcing(discrete_forcing; discrete_form=true),
]
for (label, T_forcing) in T_forcings
@info label
model = HydrostaticFreeSurfaceModel(
grid; forcing=(; T=T_forcing), tracers=(:T, :S), buoyancy=SeawaterBuoyancy()
)
set!(model; u=0.0, v=0.0, T=(x, y, z) -> target(x, y, z, nothing), S=0.0)
simulation = Simulation(model; Δt=10minute, stop_time=1day)
add_callback!(simulation, message, IterationInterval(12))
run!(simulation)
endThis outputs (using Oceananigans v0.105.2 on Julia v1.12.4)
[ Info: Relaxation
[ Info: Initializing simulation...
[ Info: Time: 0 seconds, max(|T|): 9.95e+00
[ Info: ... simulation initialization complete (8.365 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.820 seconds).
[ Info: Time: 2 hours, max(|T|): 1.01e+01
[ Info: Time: 4 hours, max(|T|): 1.03e+01
[ Info: Time: 6 hours, max(|T|): 1.05e+01
[ Info: Time: 8 hours, max(|T|): 1.07e+01
[ Info: Time: 10 hours, max(|T|): 1.09e+01
[ Info: Time: 12 hours, max(|T|): 1.12e+01
[ Info: Time: 14 hours, max(|T|): 1.14e+01
[ Info: Time: 16 hours, max(|T|): 1.17e+01
[ Info: Time: 18 hours, max(|T|): 1.19e+01
[ Info: Time: 20 hours, max(|T|): 1.20e+01
[ Info: Time: 22 hours, max(|T|): 1.21e+01
[ Info: Simulation is stopping after running for 29.657 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.
[ Info: Time: 1 day, max(|T|): 1.22e+01
[ Info: Continuous form
[ Info: Initializing simulation...
[ Info: Time: 0 seconds, max(|T|): 9.95e+00
[ Info: ... simulation initialization complete (142.936 ms)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.055 seconds).
[ Info: Time: 2 hours, max(|T|): 1.01e+01
[ Info: Time: 4 hours, max(|T|): 1.03e+01
[ Info: Time: 6 hours, max(|T|): 1.05e+01
[ Info: Time: 8 hours, max(|T|): 1.07e+01
[ Info: Time: 10 hours, max(|T|): 1.09e+01
[ Info: Time: 12 hours, max(|T|): 1.12e+01
[ Info: Time: 14 hours, max(|T|): 1.14e+01
[ Info: Time: 16 hours, max(|T|): 1.17e+01
[ Info: Time: 18 hours, max(|T|): 1.19e+01
[ Info: Time: 20 hours, max(|T|): 1.20e+01
[ Info: Time: 22 hours, max(|T|): 1.21e+01
[ Info: Simulation is stopping after running for 17.647 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.
[ Info: Time: 1 day, max(|T|): 1.22e+01
[ Info: Discrete form
[ Info: Initializing simulation...
[ Info: Time: 0 seconds, max(|T|): 9.95e+00
[ Info: ... simulation initialization complete (160.312 ms)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.033 seconds).
[ Info: Time: 2 hours, max(|T|): 9.95e+00
[ Info: Time: 4 hours, max(|T|): 9.95e+00
[ Info: Time: 6 hours, max(|T|): 9.95e+00
[ Info: Time: 8 hours, max(|T|): 9.95e+00
[ Info: Time: 10 hours, max(|T|): 9.95e+00
[ Info: Time: 12 hours, max(|T|): 9.95e+00
[ Info: Time: 14 hours, max(|T|): 9.95e+00
[ Info: Time: 16 hours, max(|T|): 9.95e+00
[ Info: Time: 18 hours, max(|T|): 9.95e+00
[ Info: Time: 20 hours, max(|T|): 9.95e+00
[ Info: Time: 22 hours, max(|T|): 9.95e+00
[ Info: Simulation is stopping after running for 4.331 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.
[ Info: Time: 1 day, max(|T|): 9.95e+00
While simulation evolves equivalently when using the Relaxation and continuous forcing, using the discrete form forcing gives different behaviour, with the forcing not appearing to have any effect on the temperature field.