Skip to content

Commit 5be36f5

Browse files
simone-silvestrisimone-silvestri
committed
add a prescribed ocean
1 parent a89e485 commit 5be36f5

File tree

8 files changed

+126
-76
lines changed

8 files changed

+126
-76
lines changed

‎experiments/flux_climatology/flux_climatology.jl‎

Lines changed: 0 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -115,32 +115,6 @@ end
115115
##### A prescribed ocean...
116116
#####
117117

118-
struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing}
119-
architecture :: A
120-
grid :: G
121-
clock :: Clock{C}
122-
velocities :: U
123-
tracers :: T
124-
timeseries :: F
125-
end
126-
127-
function PrescribedOcean(timeseries;
128-
grid,
129-
clock=Clock{Float64}(time = 0))
130-
131-
τx = Field{Face, Center, Nothing}(grid)
132-
τy = Field{Center, Face, Nothing}(grid)
133-
Jᵀ = Field{Center, Center, Nothing}(grid)
134-
Jˢ = Field{Center, Center, Nothing}(grid)
135-
136-
u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx)))
137-
v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy)))
138-
T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
139-
S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))
140-
141-
PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
142-
end
143-
144118
# ...with prescribed velocity and tracer fields
145119
version = ECCO4Monthly()
146120
dates = all_ECCO_dates(version)[1:24]
@@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch)
157131
ocean_model = PrescribedOcean((; u, v, T, S); grid)
158132
ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days)
159133

160-
#####
161-
##### Need to extend a couple of methods
162-
#####
163-
164-
function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true)
165-
tick!(model.clock, Δt)
166-
time = Time(model.clock.time)
167-
168-
possible_fts = merge(model.velocities, model.tracers)
169-
time_series_tuple = extract_field_time_series(possible_fts)
170-
171-
for fts in time_series_tuple
172-
update_field_time_series!(fts, time)
173-
end
174-
175-
parent(model.velocities.u) .= parent(model.timeseries.u[time])
176-
parent(model.velocities.v) .= parent(model.timeseries.v[time])
177-
parent(model.tracers.T) .= parent(model.timeseries.T[time])
178-
parent(model.tracers.S) .= parent(model.timeseries.S[time])
179-
180-
return nothing
181-
end
182-
183-
update_state!(::PrescribedOcean) = nothing
184-
timestepper(::PrescribedOcean) = nothing
185-
186-
reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6
187-
heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6
188-
189134
#####
190135
##### A prescribed atmosphere...
191136
#####

‎src/ClimaOcean.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,16 @@ using DataDeps
4141
using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries
4242
using Oceananigans.Grids: node
4343

44+
# Does not really matter if there is no model
45+
reference_density(::Nothing) = 0
46+
heat_capacity(::Nothing) = 0
47+
48+
reference_density(unsupported) =
49+
throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))
50+
51+
heat_capacity(unsupported) =
52+
throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))
53+
4454
const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries,
4555
GPUAdaptedFieldTimeSeries}
4656

‎src/OceanSeaIceModels/OceanSeaIceModels.jl‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,6 @@ using ClimaSeaIce: SeaIceModel
2222
using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
2323

2424
using ClimaOcean: stateindex
25-
2625
using KernelAbstractions: @kernel, @index
2726
using KernelAbstractions.Extras.LoopInfo: @unroll
2827

‎src/OceanSeaIceModels/ocean_sea_ice_model.jl‎

Lines changed: 0 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -73,25 +73,6 @@ function initialize!(model::OSIM)
7373
return nothing
7474
end
7575

76-
reference_density(unsupported) =
77-
throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))
78-
79-
heat_capacity(unsupported) =
80-
throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))
81-
82-
reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation)
83-
reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state)
84-
reference_density(eos::TEOS10EquationOfState) = eos.reference_density
85-
reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density
86-
87-
heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation)
88-
heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state)
89-
heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity
90-
91-
# Does not really matter if there is no model
92-
reference_density(::Nothing) = 0
93-
heat_capacity(::Nothing) = 0
94-
9576
function heat_capacity(::TEOS10EquationOfState{FT}) where FT
9677
cₚ⁰ = SeawaterPolynomials.TEOS10.teos10_reference_heat_capacity
9778
return convert(FT, cₚ⁰)

‎src/OceanSimulations/OceanSimulations.jl‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
module OceanSimulations
22

3-
export ocean_simulation
3+
export ocean_simulation, PrescribedOceanModel
44

55
using Oceananigans
66
using Oceananigans.Units
@@ -21,6 +21,15 @@ using Oceananigans.BuoyancyFormulations: g_Earth
2121
using Oceananigans.Coriolis: Ω_Earth
2222
using Oceananigans.Operators
2323

24+
import ClimaOcean: reference_density, heat_capacity
25+
26+
reference_density(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = reference_density(ocean.model.buoyancy.formulation)
27+
reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state)
28+
reference_density(eos::TEOS10EquationOfState) = eos.reference_density
29+
30+
heat_capacity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = heat_capacity(ocean.model.buoyancy.formulation)
31+
heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state)
32+
2433
struct Default{V}
2534
value :: V
2635
end
@@ -41,5 +50,6 @@ default_or_override(override, alternative_default=nothing) = override
4150

4251
include("barotropic_potential_forcing.jl")
4352
include("ocean_simulation.jl")
53+
include("prescribed_ocean_model.jl")
4454

4555
end # module

‎src/OceanSimulations/prescribed_ocean_model.jl‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series!
2+
3+
import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick!
4+
import Oceananigans.Models: timestepper, update_model_field_time_series!
5+
6+
import ClimaOcean: reference_density, heat_capacity
7+
import Oceananigans.Architectures: on_architecture
8+
9+
#####
10+
##### A prescribed ocean...
11+
#####
12+
13+
struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing}
14+
architecture :: A
15+
grid :: G
16+
clock :: Clock{C}
17+
velocities :: U
18+
tracers :: T
19+
timeseries :: F
20+
end
21+
22+
function PrescribedOcean(timeseries;
23+
grid,
24+
clock=Clock{Float64}(time = 0))
25+
26+
# Make sure all elements of the timeseries are on the same grid
27+
for (k, v) in timeseries
28+
if !isa(v, FieldTimeSeries)
29+
throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`"))
30+
end
31+
if v.grid != grid
32+
throw(ArgumentError("All variables in the timeseries reside on the provided grid"))
33+
end
34+
end
35+
36+
τx = Field{Face, Center, Nothing}(grid)
37+
τy = Field{Center, Face, Nothing}(grid)
38+
Jᵀ = Field{Center, Center, Nothing}(grid)
39+
Jˢ = Field{Center, Center, Nothing}(grid)
40+
41+
u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx)))
42+
v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy)))
43+
T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
44+
S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))
45+
46+
PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
47+
end
48+
49+
#####
50+
##### Need to extend a couple of methods
51+
#####
52+
53+
function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true)
54+
tick!(model.clock, Δt)
55+
time = Time(model.clock.time)
56+
57+
possible_fts = merge(model.velocities, model.tracers)
58+
time_series_tuple = extract_field_time_series(possible_fts)
59+
60+
for fts in time_series_tuple
61+
update_field_time_series!(fts, time)
62+
end
63+
64+
update_u_velocity = haskey(model.timeseries, :u)
65+
update_v_velocity = haskey(model.timeseries, :v)
66+
update_temperature = haskey(model.timeseries, :T)
67+
update_salinity = haskey(model.timeseries, :S)
68+
69+
update_u_velocity && iterpolate!(model.velocities.u, model.timeseries.u[time])
70+
update_v_velocity && iterpolate!(model.velocities.v, model.timeseries.v[time])
71+
update_temperature && iterpolate!(model.tracers.T, model.timeseries.T[time])
72+
update_salinity && iterpolate!(model.tracers.S, model.timeseries.S[time])
73+
74+
return nothing
75+
end
76+
77+
update_state!(::PrescribedOcean) = nothing
78+
timestepper(::PrescribedOcean) = nothing
79+
80+
reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6
81+
heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6

‎src/SeaIceSimulations.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,11 @@ using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium
2020

2121
using ClimaOcean.OceanSimulations: Default
2222

23+
import ClimaOcean: reference_density, heat_capacity
24+
25+
reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density
26+
heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity
27+
2328
function sea_ice_simulation(grid;
2429
Δt = 5minutes,
2530
ice_salinity = 0, # psu

‎test/test_ocean_sea_ice_model.jl‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -100,5 +100,24 @@ using ClimaSeaIce.Rheologies
100100
true
101101
end
102102
end
103+
104+
@testset "Test a prescribed ocean model" begin
105+
grid = LatitudeLongitudeGrid(size = (10, 10, 10), latitude = (20, 30), longitude = (40, 50), z = (-100, 0))
106+
T = ECCOFieldTimeSeries(:temperature, grid; time_indices_in_memory = 2)
107+
S = ECCOFieldTimeSeries(:salinity, grid; time_indices_in_memory = 2)
108+
ocean_model = PrescribedOcean((; T, S); grid)
109+
110+
time_step!(ocean_model, T.times[2])
111+
112+
@test ocean_model.clock.time == T.times[2]
113+
@test ocean_model.tracers.T.data == T[2].data
114+
@test ocean_model.tracers.S.data == S[2].data
115+
116+
time_step!(ocean_model, 10days)
117+
118+
@test ocean_model.clock.time == 10days
119+
@test T.backend.start == 2
120+
@test ocean_model.tracers.T.data != T[3].data
121+
end
103122
end
104123

0 commit comments

Comments
 (0)