Calibrate IsopycnalSkewSymmetricDiffusivity in eddying channel
#159
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
The script Running the script gives: julia> using OceanTurbulenceParameterEstimation
julia> using Oceananigans
julia> using Oceananigans.Units
julia> using Oceananigans.Models.HydrostaticFreeSurfaceModels: SliceEnsembleSize
julia> using Oceananigans.TurbulenceClosures: FluxTapering
julia> using LinearAlgebra, CairoMakie, DataDeps, Distributions, JLD2
julia> using ElectronDisplay
julia> architecture = CPU()
CPU()
julia> # download from https://www.dropbox.com/s/91altratyy1g0fc/eddying_channel_catke_zonal_average.jld2?dl=0
# and change path below accordingly
filepath = "/Users/navid/Research/mesoscale-parametrization-OSM2022/eddying_channel/Ny200Nx100_Lx1000_Ly2000/eddying_channel_catke_zonal_average.jld2"
"/Users/navid/Research/mesoscale-parametrization-OSM2022/eddying_channel/Ny200Nx100_Lx1000_Ly2000/eddying_channel_catke_zonal_average.jld2"
julia> b_timeseries = FieldTimeSeries(filepath, "b")
1×200×60×1043 FieldTimeSeries{InMemory} located at (⋅, Center, Center)
├── architecture: Float64
└── grid: 100×200×60 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
julia> field_names = (:b, :u, :v, :w)
(:b, :u, :v, :w)
julia> # field_names = (:b, :c, :u, :v, :w)
normalization = (b = ZScore(),
# c = ZScore(),
u = ZScore(),
v = ZScore(),
w = RescaledZScore(1e-2))
(b = ZScore{Nothing}(nothing, nothing), u = ZScore{Nothing}(nothing, nothing), v = ZScore{Nothing}(nothing, nothing), w = RescaledZScore{Float64, Nothing}(0.01, nothing))
julia> times = b_timeseries.times[500:10:800]
31-element Vector{Float64}:
3.017952e8
3.078432e8
3.138912e8
3.199392e8
3.259872e8
3.320352e8
3.380832e8
3.441312e8
3.501792e8
3.562272e8
3.622752e8
3.683232e8
3.743712e8
3.804192e8
3.864672e8
3.925152e8
3.985632e8
4.046112e8
4.106592e8
4.167072e8
4.227552e8
4.288032e8
4.348512e8
4.408992e8
4.469472e8
4.529952e8
4.590432e8
4.650912e8
4.711392e8
4.771872e8
4.832352e8
julia> observations = SyntheticObservations(filepath; normalization, times, field_names)
SyntheticObservations with fields (:b, :u, :v, :w)
├── times: [3.017952e8, 3.078432e8, 3.138912e8, 3.199392e8, 3.259872e8, 3.320352e8, 3.380832e8, 3.441312e8, 3.501792e8, 3.562272e8, 3.622752e8, 3.683232e8, 3.743712e8, 3.804192e8, 3.864672e8, 3.925152e8, 3.985632e8, 4.046112e8, 4.106592e8, 4.167072e8, 4.227552e8, 4.288032e8, 4.348512e8, 4.408992e8, 4.469472e8, 4.529952e8, 4.590432e8, 4.650912e8, 4.711392e8, 4.771872e8, 4.832352e8]
├── grid: 100×200×60 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── path: "/Users/navid/Research/mesoscale-parametrization-OSM2022/eddying_channel/Ny200Nx100_Lx1000_Ly2000/eddying_channel_catke_zonal_average.jld2"
├── metadata: (:grid, :coriolis, :closure)
└── normalization: Dict{Symbol, Any} with 4 entries
julia> #####
##### Simulation
#####
file = jldopen(filepath)
JLDFile /Users/navid/Research/mesoscale-parametrization-OSM2022/eddying_channel/Ny200Nx100_Lx1000_Ly2000/eddying_channel_catke_zonal_average.jld2 (read-only)
├─📂 grid
│ ├─🔢 Nx
│ ├─🔢 Ny
│ ├─🔢 Nz
│ ├─🔢 Hx
│ ├─🔢 Hy
│ ├─🔢 Hz
│ ├─🔢 Lx
│ └─ ⋯ (14 more entries)
└─ ⋯ (4 more entries)
julia> coriolis = file["serialized/coriolis"]
BetaPlane{Float64}: f₀ = -1.00e-04, β = 1.00e-11
julia> # Domain
const Ly, Lz= file["grid/Ly"], file["grid/Lz"]
(2.0e6, 2000.0)
julia> # number of grid points
Ny, Nz= file["grid/Ny"], file["grid/Nz"]
(200, 60)
julia> Nensemble = 6
6
julia> slice_ensemble_size = SliceEnsembleSize(size=(Ny, Nz), ensemble=Nensemble)
SliceEnsembleSize(6, 200, 60, 1, 1)
julia> ensemble_grid = RectilinearGrid(architecture,
size = slice_ensemble_size,
topology = (Flat, Bounded, Bounded),
y = (0, Ly),
z = (-Lz, 0),
halo=(3, 3))
6×200×60 RectilinearGrid{Float64, Flat, Bounded, Bounded} on CPU with 0×3×3 halo
├── Flat x
├── Bounded y ∈ [0.0, 2.0e6] regularly spaced with Δy=10000.0
└── Bounded z ∈ [-2000.0, 0.0] regularly spaced with Δz=33.3333
julia> κ_skew = 1000.0 # [m² s⁻¹] skew diffusivity
1000.0
julia> κ_symmetric = 900.0 # [m² s⁻¹] symmetric diffusivity
900.0
julia> gerdes_koberle_willebrand_tapering = FluxTapering(1e-2)
FluxTapering{Float64}(0.01)
julia> gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew = κ_skew,
κ_symmetric = κ_symmetric,
slope_limiter = gerdes_koberle_willebrand_tapering)
IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01))
julia> closures = file["serialized/closure"]
(Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
julia> closure_ensemble = [(deepcopy(gent_mcwilliams_diffusivity), closures[1], closures[2]) for _ = 1:Nensemble]
6-element Vector{Tuple{IsopycnalSkewSymmetricDiffusivity{Float64, Float64, Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}, FluxTapering{Float64}}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}, AnisotropicDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Float64, Float64, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}}}}:
(IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01)), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
(IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01)), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
(IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01)), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
(IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01)), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
(IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01)), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
(IsopycnalSkewSymmetricDiffusivity: (κ_symmetric=900.0, κ_skew=1000.0, (isopycnal_tensor=Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}(0), slope_limiter=FluxTapering{Float64}(0.01)), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}(2.91, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}(1.16, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.15, 3.87, 0.4, 0.77, 0.13, 1.11, 0.72, 0.76), Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}(3.62, 1.31)), AnisotropicDiffusivity: (νx=30.0, νy=30.0, νz=0.0003), (κx=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κy=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6), κz=(b = 5.0e-6, e = 5.0e-6, c = 5.0e-6)))
julia> ensemble_model = HydrostaticFreeSurfaceModel(grid = ensemble_grid,
tracers = (:b, :e, :c),
buoyancy = BuoyancyTracer(),
coriolis = coriolis,
closure = closure_ensemble,
free_surface = ImplicitFreeSurface(),
)
ERROR: MethodError: no method matching with_tracers(::Tuple{Symbol, Symbol, Symbol}, ::Vector{Tuple{IsopycnalSkewSymmetricDiffusivity{Float64, Float64, Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}, FluxTapering{Float64}}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}, AnisotropicDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Float64, Float64, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}}}})
Closest candidates are:
with_tracers(::Any, ::IsopycnalSkewSymmetricDiffusivity) at /Users/navid/.julia/packages/Oceananigans/IHPoE/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl:35
with_tracers(::Any, ::ConvectiveAdjustmentVerticalDiffusivity{TD, CK, CN, BK, BN} where {CK, CN, BK, BN}) where TD at /Users/navid/.julia/packages/Oceananigans/IHPoE/src/TurbulenceClosures/turbulence_closure_implementations/convective_adjustment_vertical_diffusivity.jl:67
with_tracers(::Any, ::NamedTuple, ::Any; with_velocities) at /Users/navid/.julia/packages/Oceananigans/IHPoE/src/Utils/with_tracers.jl:10
...
Stacktrace:
[1] HydrostaticFreeSurfaceModel(; grid::RectilinearGrid{Float64, Flat, Bounded, Bounded, Float64, Float64, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, CPU}, clock::Clock{Float64}, momentum_advection::CenteredSecondOrder, tracer_advection::CenteredSecondOrder, buoyancy::BuoyancyTracer, coriolis::BetaPlane{Float64}, free_surface::ImplicitFreeSurface{Nothing, Float64, Nothing, Nothing, Symbol, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, forcing::NamedTuple{(), Tuple{}}, closure::Vector{Tuple{IsopycnalSkewSymmetricDiffusivity{Float64, Float64, Oceananigans.TurbulenceClosures.SmallSlopeIsopycnalTensor{Int64}, FluxTapering{Float64}}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.CATKEVerticalDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.MixingLength{Float64}, Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities.SurfaceTKEFlux{Float64}}, AnisotropicDiffusivity{Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization, Float64, Float64, Float64, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:b, :e, :c), Tuple{Float64, Float64, Float64}}}}}, boundary_conditions::NamedTuple{(), Tuple{}}, tracers::Tuple{Symbol, Symbol, Symbol}, particles::Nothing, velocities::Nothing, pressure::Nothing, diffusivity_fields::Nothing, auxiliary_fields::NamedTuple{(), Tuple{}})
@ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/.julia/packages/Oceananigans/IHPoE/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:152
[2] top-level scope
@ REPL[87]:1
julia> |
|
(Note, we don't need to implement the most general solution possible. I'd be happy to implement a workaround that, e.g., works for 3 closures but not for 4!) |
|
The problem is |
|
Should we support a vector of closure tuples or a tuple of vectors? |
|
Why do you think it's related with (But |
For single closures case we put them in a vector that is |
I suspect it's better to always have tuples of closures (and sometimes the closures are singletons, other times they are arrays). But we'll see... |
|
Sure. But example still not working. |
|
Maybe just open a new PR from the same branch? |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This is the first attempt to calibrate
IsopycnalSkewSymmetricDiffusivityusing output from the eddying channel solution.I am running into problems (briefly discussed also in #158) when constructing ensemble model with more than one closures. What I would like is to create an ensemble model in which each model should have 3 closures:
CATKE+AnisotropicDiffusivitythat are the same for all ensemble members andIsopycnalSkewSymmetricDiffusivitythat vary for each ensemble member.cc @glwagner