Skip to content
Open
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 0 additions & 55 deletions experiments/flux_climatology/flux_climatology.jl
Original file line number Diff line number Diff line change
Expand Up @@ -115,32 +115,6 @@ end
##### A prescribed ocean...
#####

struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing}
architecture :: A
grid :: G
clock :: Clock{C}
velocities :: U
tracers :: T
timeseries :: F
end

function PrescribedOcean(timeseries;
grid,
clock=Clock{Float64}(time = 0))

τx = Field{Face, Center, Nothing}(grid)
τy = Field{Center, Face, Nothing}(grid)
Jᵀ = Field{Center, Center, Nothing}(grid)
Jˢ = Field{Center, Center, Nothing}(grid)

u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx)))
v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy)))
T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))

PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
end

# ...with prescribed velocity and tracer fields
version = ECCO4Monthly()
dates = all_ECCO_dates(version)[1:24]
Expand All @@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch)
ocean_model = PrescribedOcean((; u, v, T, S); grid)
ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days)

#####
##### Need to extend a couple of methods
#####

function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true)
tick!(model.clock, Δt)
time = Time(model.clock.time)

possible_fts = merge(model.velocities, model.tracers)
time_series_tuple = extract_field_time_series(possible_fts)

for fts in time_series_tuple
update_field_time_series!(fts, time)
end

parent(model.velocities.u) .= parent(model.timeseries.u[time])
parent(model.velocities.v) .= parent(model.timeseries.v[time])
parent(model.tracers.T) .= parent(model.timeseries.T[time])
parent(model.tracers.S) .= parent(model.timeseries.S[time])

return nothing
end

update_state!(::PrescribedOcean) = nothing
timestepper(::PrescribedOcean) = nothing

reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6
heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6

#####
##### A prescribed atmosphere...
#####
Expand Down
10 changes: 10 additions & 0 deletions src/ClimaOcean.jl
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,16 @@ using DataDeps
using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries
using Oceananigans.Grids: node

# Does not really matter if there is no model
reference_density(::Nothing) = 0
heat_capacity(::Nothing) = 0

reference_density(unsupported) =
throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))

heat_capacity(unsupported) =
throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))

const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries,
GPUAdaptedFieldTimeSeries}

Expand Down
1 change: 0 additions & 1 deletion src/OceanSeaIceModels/OceanSeaIceModels.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ using ClimaSeaIce: SeaIceModel
using ClimaSeaIce.SeaIceThermodynamics: melting_temperature

using ClimaOcean: stateindex

using KernelAbstractions: @kernel, @index
using KernelAbstractions.Extras.LoopInfo: @unroll

Expand Down
19 changes: 0 additions & 19 deletions src/OceanSeaIceModels/ocean_sea_ice_model.jl
Original file line number Diff line number Diff line change
Expand Up @@ -73,25 +73,6 @@ function initialize!(model::OSIM)
return nothing
end

reference_density(unsupported) =
throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))

heat_capacity(unsupported) =
throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))

reference_density(ocean::Simulation) = reference_density(ocean.model.buoyancy.formulation)
reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state)
reference_density(eos::TEOS10EquationOfState) = eos.reference_density
reference_density(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density

heat_capacity(ocean::Simulation) = heat_capacity(ocean.model.buoyancy.formulation)
heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state)
heat_capacity(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity

# Does not really matter if there is no model
reference_density(::Nothing) = 0
heat_capacity(::Nothing) = 0

function heat_capacity(::TEOS10EquationOfState{FT}) where FT
cₚ⁰ = SeawaterPolynomials.TEOS10.teos10_reference_heat_capacity
return convert(FT, cₚ⁰)
Expand Down
12 changes: 11 additions & 1 deletion src/OceanSimulations/OceanSimulations.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module OceanSimulations

export ocean_simulation
export ocean_simulation, PrescribedOceanModel

using Oceananigans
using Oceananigans.Units
Expand All @@ -21,6 +21,15 @@ using Oceananigans.BuoyancyFormulations: g_Earth
using Oceananigans.Coriolis: Ω_Earth
using Oceananigans.Operators

import ClimaOcean: reference_density, heat_capacity

reference_density(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = reference_density(ocean.model.buoyancy.formulation)
reference_density(buoyancy_formulation::SeawaterBuoyancy) = reference_density(buoyancy_formulation.equation_of_state)
reference_density(eos::TEOS10EquationOfState) = eos.reference_density

heat_capacity(ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = heat_capacity(ocean.model.buoyancy.formulation)
heat_capacity(buoyancy_formulation::SeawaterBuoyancy) = heat_capacity(buoyancy_formulation.equation_of_state)

struct Default{V}
value :: V
end
Expand All @@ -41,5 +50,6 @@ default_or_override(override, alternative_default=nothing) = override

include("barotropic_potential_forcing.jl")
include("ocean_simulation.jl")
include("prescribed_ocean_model.jl")

end # module
96 changes: 96 additions & 0 deletions src/OceanSimulations/prescribed_ocean_model.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
using Oceananigans.Models: AbstractModel
using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series!

import Oceananigans.TimeSteppers: time_step!, update_state!, reset!, tick!
import Oceananigans.Models: timestepper, update_model_field_time_series!

import ClimaOcean: reference_density, heat_capacity
import Oceananigans.Architectures: on_architecture

#####
##### A prescribed ocean...
#####

struct PrescribedOceanModel{G, C, U, T, F, Arch} <: AbstractModel{Nothing, Arch}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We call it PrescribedAtmosphere. Is this PrescribedOcean or PrescribedOceanModel?

I'm not sure really sure what a "prescribed model" is.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, this is because we add a simulation on top, while the atmosphere is a "simulation" rather than a model. I can switch it to a simulation rather than a model

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is because the infrastructure we have now expects an ocean with a pattern of ocean.model.stuff..., while there is no assumption on the atmosphere. I did not want to make a huge refactor here, I can call it PrescribedOcean though and just have one field called model inside of it

architecture :: Arch
grid :: G
clock :: Clock{C}
velocities :: U
tracers :: T
timeseries :: F
end

"""
PrescribedOceanModel(timeseries=NamedTuple(); grid, clock=Clock{Float64}(time = 0))

Create a prescribed ocean model to be used in combination with ClimaOcean's `OceanSeaIceModel`
on a `grid` with a `clock`.

Arguments
=========
- `timeseries`: A `NamedTuple` containing time series data for various fields. The named tuple can be empty
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the name is confusing because when we say "timeseries" it means just one series. But this is many timeseries?

(meaning that the ocean does not evolve in time) or include any combination of the
following fields: `u`, `v`, `T`, `S`. All elements provided must be of type `FieldTimeSeries`
and reside on the provided `grid`.
"""
function PrescribedOceanModel(timeseries=NamedTuple();
grid,
clock=Clock{Float64}(time = 0))

# Make sure all elements of the timeseries are on the same grid
# If we decide to pass a timeseries
if !isempty(timeseries)
for (k, v) in timeseries
isa(v, FieldTimeSeries) ||
throw(ArgumentError("All variables in the `timeseries` argument must be `FieldTimeSeries`"))
v.grid == grid ||
throw(ArgumentError("All variables in the timeseries reside on the provided grid"))
end
end

τx = Field{Face, Center, Nothing}(grid)
τy = Field{Center, Face, Nothing}(grid)
Jᵀ = Field{Center, Center, Nothing}(grid)
Jˢ = Field{Center, Center, Nothing}(grid)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why this? I think this is wrong


u = XFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Face, Center, Center), top = FluxBoundaryCondition(τx)))
v = YFaceField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Face, Center), top = FluxBoundaryCondition(τy)))
T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))

return PrescribedOceanModel(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
end

#####
##### Need to extend a couple of methods
#####

function time_step!(model::PrescribedOceanModel, Δt; callbacks=[], euler=true)
tick!(model.clock, Δt)
time = Time(model.clock.time)

possible_fts = merge(model.velocities, model.tracers)
time_series_tuple = extract_field_time_series(possible_fts)

for fts in time_series_tuple
update_field_time_series!(fts, time)
end

update_u_velocity = haskey(model.timeseries, :u)
update_v_velocity = haskey(model.timeseries, :v)
update_temperature = haskey(model.timeseries, :T)
update_salinity = haskey(model.timeseries, :S)

update_u_velocity && parent(model.velocities.u) .= parent(model.timeseries.u[time])
update_v_velocity && parent(model.velocities.v) .= parent(model.timeseries.v[time])
update_temperature && parent(model.tracers.T) .= parent(model.timeseries.T[time])
update_salinity && parent(model.tracers.S) .= parent(model.timeseries.S[time])

return nothing
end

update_state!(::PrescribedOceanModel) = nothing
timestepper(::PrescribedOceanModel) = nothing

reference_density(ocean::Simulation{<:PrescribedOceanModel}) = 1025.6
heat_capacity(ocean::Simulation{<:PrescribedOceanModel}) = 3995.6
5 changes: 5 additions & 0 deletions src/SeaIceSimulations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,11 @@ using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium

using ClimaOcean.OceanSimulations: Default

import ClimaOcean: reference_density, heat_capacity

reference_density(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_density
heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.ice_heat_capacity

function sea_ice_simulation(grid;
Δt = 5minutes,
ice_salinity = 0, # psu
Expand Down
19 changes: 19 additions & 0 deletions test/test_ocean_sea_ice_model.jl
Original file line number Diff line number Diff line change
Expand Up @@ -100,5 +100,24 @@ using ClimaSeaIce.Rheologies
true
end
end

@testset "Test a prescribed ocean model" begin
grid = LatitudeLongitudeGrid(size = (10, 10, 10), latitude = (20, 30), longitude = (40, 50), z = (-100, 0))
T = ECCOFieldTimeSeries(:temperature, grid; time_indices_in_memory = 2)
S = ECCOFieldTimeSeries(:salinity, grid; time_indices_in_memory = 2)
ocean_model = PrescribedOcean((; T, S); grid)

time_step!(ocean_model, T.times[2])

@test ocean_model.clock.time == T.times[2]
@test ocean_model.tracers.T.data == T[2].data
@test ocean_model.tracers.S.data == S[2].data

time_step!(ocean_model, 10days)

@test ocean_model.clock.time == 10days
@test T.backend.start == 2
@test ocean_model.tracers.T.data != T[3].data
end
end