From 02b697c39efccbf5b8ea375c3c63c9bc5592f302 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 5 Jul 2023 12:41:29 -0500 Subject: [PATCH 01/49] Set up grid with immersed boundary condition for ASC model, starting to setup model itself --- examples/antarctic_slope_current_model.jl | 89 +++++++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 examples/antarctic_slope_current_model.jl diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl new file mode 100644 index 00000000..538e51fb --- /dev/null +++ b/examples/antarctic_slope_current_model.jl @@ -0,0 +1,89 @@ +using Oceananigans + +# This files sets up a model that resembles the Antarctic Slope Current (ASC) model in the +# 2022 paper by Ian Eisenman + +arch = CPU() + +g_Earth = 9.80665 + +Lx = 400000 +Ly = 450000 +Lz = 4000 + +Nx = 400 +Ny = 450 +Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor + +# +# Setting up the grid and bathymetry: +# +# Using a linear slope that approximates the min and max spacings in the paper: +init_adjustment = (0.5*(100-10)/(Nz-1)) + (10 - (100-10)/(Nz-1)) +linear_slope(k) = (0.5*(100-10)/(Nz-1))*(k)^2 + (10 - (100-10)/(Nz-1))*(k) - init_adjustment + +spacing_adjustment = Lz / linear_slope(Nz+1) +linear_slope_z_faces(k) = -spacing_adjustment * linear_slope(k) +# Can reverse this to get grid from -4000 to 0 later + +underlying_grid = RectilinearGrid(arch, + size = (Nx, Ny, Nz), + topology = (Periodic, Bounded, Bounded), + x = (-Lx/2, Lx/2), + y = (0, Ly), + z = linear_slope_z_faces) + +println(underlying_grid) + +# We want the underwater slope to provide a depth of 500 m at y = 0 and the full 4 km by y =200. It follows +# a hyperbolic tangent curve with its center at y ~= 150 at a depth of ~ -4000 + (4000 - 500)/2 = -2250 m +underwater_slope(x, y) = -1750tanh(0.00004*(y - 150000)) - 2250 +# TODO: add 50km troughs to the slope, dependent on x and y. Use appendix C to add detail here + +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) + +println(grid) + +# TODO: add underwater slope at y = 0 with troughs + +# +# Running the Hydrostatic model: +# + +# Forcings: +u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) + +# Boundary Conditions: +no_slip_bc = ValueBoundaryCondition(0.0) +free_slip_bc = FluxBoundaryCondition(nothing) + +free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)) +no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc) + + +# Assuming no particles or biogeochemistry +model = HydrostaticFreeSurfaceModel(; grid, + clock = Clock{eltype(grid)}(0, 0, 1), + momentum_advection = CenteredSecondOrder(), + tracer_advection = CenteredSecondOrder(), + buoyancy = SeawaterBuoyancy(eltype(grid)), + coriolis = nothing, + free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), + #forcing = (,), # NamedTuple(), + closure = nothing, + boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs), # NamedTuple(), + tracers = (:T, :S), + velocities = nothing, + pressure = nothing, + diffusivity_fields = nothing, + #auxiliary_fields = nothing, # NamedTuple(), +) + +println(model) + +println("u boundary conditions:") +println(model.velocities.u.boundary_conditions) +println("v boundary conditions:") +println(model.velocities.v.boundary_conditions) +println("w boundary conditions:") +println(model.velocities.w.boundary_conditions) From 3f7b23aef6f8bcf73735fada10c4d83fddb5a3a5 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 5 Jul 2023 15:29:07 -0500 Subject: [PATCH 02/49] Adding sponge layers and coriolis, updated equations of state to be high order polynomials --- examples/antarctic_slope_current_model.jl | 33 +++++++++++++++++++---- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 538e51fb..35deeb02 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,6 +1,8 @@ + using Oceananigans +using SeawaterPolynomials.TEOS10 -# This files sets up a model that resembles the Antarctic Slope Current (ASC) model in the +# This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Ian Eisenman arch = CPU() @@ -15,6 +17,8 @@ Nx = 400 Ny = 450 Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +sponge_width = 20000 + # # Setting up the grid and bathymetry: # @@ -51,7 +55,21 @@ println(grid) # # Forcings: -u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) +u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) # Not actual forcing, just example + +# TODO: need to use forcings to enact the sponge layers. Their support is within 20000 m +# of the N/S boundaries, they impose a cross-slope buoyancy gradient, and the relaxation +# tiem scales decrease linearly with distance from the interior termination of the sponge layers + +# We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity +damping_rate = 1 / 100 # Relaxation time scale of 100s, need to make this decrease linearly toward outermost boundary +south_mask = GaussianMask{:y}(center=0, width=sponge_width) +north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) +south_sponge_layer = Relaxation(; rate=damping_rate, mask=south_mask) +north_sponge_layer = Relaxation(; rate=damping_rate, mask=north_mask) +sponge_layers = south_sponge_layer +# TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers + # Boundary Conditions: no_slip_bc = ValueBoundaryCondition(0.0) @@ -60,16 +78,21 @@ free_slip_bc = FluxBoundaryCondition(nothing) free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)) no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc) +# Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 +eos = TEOS10EquationOfState() + +# Coriolis Effect, using basic f-plane +coriolis = FPlane(f=-1.3e14) # Assuming no particles or biogeochemistry model = HydrostaticFreeSurfaceModel(; grid, clock = Clock{eltype(grid)}(0, 0, 1), momentum_advection = CenteredSecondOrder(), tracer_advection = CenteredSecondOrder(), - buoyancy = SeawaterBuoyancy(eltype(grid)), - coriolis = nothing, + buoyancy = SeawaterBuoyancy(equation_of_state=eos), + coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), - #forcing = (,), # NamedTuple(), + forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = nothing, boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs), # NamedTuple(), tracers = (:T, :S), From 375b630565098bafd3130b86dd428c761c4d0163 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 5 Jul 2023 15:48:17 -0500 Subject: [PATCH 03/49] Added both sponge layers and relaxation forcings, and diffusivities as closures --- examples/antarctic_slope_current_model.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 35deeb02..53146c8e 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -67,7 +67,7 @@ south_mask = GaussianMask{:y}(center=0, width=sponge_width) north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) south_sponge_layer = Relaxation(; rate=damping_rate, mask=south_mask) north_sponge_layer = Relaxation(; rate=damping_rate, mask=north_mask) -sponge_layers = south_sponge_layer +sponge_layers = (south_sponge_layer, north_sponge_layer) # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers @@ -81,9 +81,14 @@ no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc) # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 eos = TEOS10EquationOfState() -# Coriolis Effect, using basic f-plane +# Coriolis Effect, using basic f-plane with precribed reference Coriolis parameter coriolis = FPlane(f=-1.3e14) +# Diffusivities as part of closure +# TODO: make sure this works for biharmonic diffusivities as the horizontal, +horizontal_closure = HorizontalScalarDiffusivity(ν=0.1, κ=0.1) +vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) + # Assuming no particles or biogeochemistry model = HydrostaticFreeSurfaceModel(; grid, clock = Clock{eltype(grid)}(0, 0, 1), @@ -93,7 +98,7 @@ model = HydrostaticFreeSurfaceModel(; grid, coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() - closure = nothing, + closure = (horizontal_closure, vertical_closure), boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs), # NamedTuple(), tracers = (:T, :S), velocities = nothing, From 2e28ce6eb34382ea4d5ab5520d761f4b9543704a Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 6 Jul 2023 18:29:12 -0500 Subject: [PATCH 04/49] Created simulation run --- examples/antarctic_slope_current_model.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 53146c8e..fd56d52e 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -51,7 +51,7 @@ println(grid) # TODO: add underwater slope at y = 0 with troughs # -# Running the Hydrostatic model: +# Creating the Hydrostatic model: # # Forcings: @@ -108,10 +108,18 @@ model = HydrostaticFreeSurfaceModel(; grid, ) println(model) - +#= println("u boundary conditions:") println(model.velocities.u.boundary_conditions) println("v boundary conditions:") println(model.velocities.v.boundary_conditions) println("w boundary conditions:") println(model.velocities.w.boundary_conditions) +=# + +# +# Now create a simulation and run the model +# + +simulation = Simulation(model; Δt=0.01, stop_iteration=100) +run!(simulation) \ No newline at end of file From 952b7965ccb57349e9981fbe95a5b9a57bc076ea Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 6 Jul 2023 18:31:07 -0500 Subject: [PATCH 05/49] Fixed simulation time scale, set to run on GPU --- examples/antarctic_slope_current_model.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index fd56d52e..be682fae 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -5,7 +5,7 @@ using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Ian Eisenman -arch = CPU() +arch = GPU() g_Earth = 9.80665 @@ -62,7 +62,7 @@ u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) # Not actual forcing, just exam # tiem scales decrease linearly with distance from the interior termination of the sponge layers # We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity -damping_rate = 1 / 100 # Relaxation time scale of 100s, need to make this decrease linearly toward outermost boundary +damping_rate = 1 / 43200 # Relaxation time scale in seconds, need to make this decrease linearly toward outermost boundary south_mask = GaussianMask{:y}(center=0, width=sponge_width) north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) south_sponge_layer = Relaxation(; rate=damping_rate, mask=south_mask) @@ -121,5 +121,5 @@ println(model.velocities.w.boundary_conditions) # Now create a simulation and run the model # -simulation = Simulation(model; Δt=0.01, stop_iteration=100) +simulation = Simulation(model; Δt=100.0, stop_iteration=100) run!(simulation) \ No newline at end of file From 17da147b2a14e84c34d10bb4032966308d6d38ad Mon Sep 17 00:00:00 2001 From: Joseph Kump Date: Fri, 7 Jul 2023 11:34:12 -0500 Subject: [PATCH 06/49] Printing out simulation --- examples/antarctic_slope_current_model.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index be682fae..fac5a98b 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -122,4 +122,6 @@ println(model.velocities.w.boundary_conditions) # simulation = Simulation(model; Δt=100.0, stop_iteration=100) -run!(simulation) \ No newline at end of file +run!(simulation) + +println(simulation) \ No newline at end of file From 2c2cfea85fb1ded76d2de03210b7a5e512598abd Mon Sep 17 00:00:00 2001 From: jlk9 Date: Fri, 7 Jul 2023 12:35:16 -0500 Subject: [PATCH 07/49] Adding storage of output --- examples/antarctic_slope_current_model.jl | 42 +++++++++++++++++++---- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index fac5a98b..0ff00b28 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,11 +1,13 @@ +using CairoMakie using Oceananigans +using Oceananigans.Units: minute, minutes, hour using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Ian Eisenman -arch = GPU() +arch = CPU() g_Earth = 9.80665 @@ -18,7 +20,7 @@ Ny = 450 Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20000 - +#= # # Setting up the grid and bathymetry: # @@ -29,13 +31,28 @@ linear_slope(k) = (0.5*(100-10)/(Nz-1))*(k)^2 + (10 - (100-10)/(Nz-1))*(k) - ini spacing_adjustment = Lz / linear_slope(Nz+1) linear_slope_z_faces(k) = -spacing_adjustment * linear_slope(k) # Can reverse this to get grid from -4000 to 0 later +=# +refinement = 20.0 # controls spacing near surface (higher means finer spaced), 12.0 +stretching = 3 # controls rate of stretching at bottom, 3 + +# Normalized height ranging from 0 to 1 +h(k) = (k - 1) / Nz + +# Linear near-surface generator +ζ₀(k) = 1 + (h(k) - 1) / refinement + +# Bottom-intensified stretching function +Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching)) + +# Generating function +z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1) underlying_grid = RectilinearGrid(arch, size = (Nx, Ny, Nz), topology = (Periodic, Bounded, Bounded), x = (-Lx/2, Lx/2), y = (0, Ly), - z = linear_slope_z_faces) + z = z_faces) println(underlying_grid) @@ -108,20 +125,31 @@ model = HydrostaticFreeSurfaceModel(; grid, ) println(model) -#= + println("u boundary conditions:") println(model.velocities.u.boundary_conditions) println("v boundary conditions:") println(model.velocities.v.boundary_conditions) println("w boundary conditions:") println(model.velocities.w.boundary_conditions) -=# # # Now create a simulation and run the model # +simulation = Simulation(model; Δt=100.0, stop_time=10hours) + +# Create a NamedTuple with eddy viscosity +eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ) + +filename = "asc_model_run" + +simulation.output_writers[:slices] = + JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity), + filename = filename * ".jld2", + indices = (:, grid.Ny/2, :), + schedule = TimeInterval(1minute), + overwrite_existing = true) -simulation = Simulation(model; Δt=100.0, stop_iteration=100) run!(simulation) -println(simulation) \ No newline at end of file +println(simulation) From f7ae29f77e03f8304df209e9f5860b250407e678 Mon Sep 17 00:00:00 2001 From: Joseph Kump Date: Fri, 7 Jul 2023 16:59:19 -0500 Subject: [PATCH 08/49] Minor edits --- examples/antarctic_slope_current_model.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 0ff00b28..99dbf5c3 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,5 +1,5 @@ -using CairoMakie +#using CairoMakie using Oceananigans using Oceananigans.Units: minute, minutes, hour using SeawaterPolynomials.TEOS10 @@ -7,7 +7,7 @@ using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Ian Eisenman -arch = CPU() +arch = GPU() g_Earth = 9.80665 @@ -136,15 +136,14 @@ println(model.velocities.w.boundary_conditions) # # Now create a simulation and run the model # -simulation = Simulation(model; Δt=100.0, stop_time=10hours) +simulation = Simulation(model; Δt=100.0, stop_time=600minutes) -# Create a NamedTuple with eddy viscosity -eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ) +# Create a NamedTuple filename = "asc_model_run" simulation.output_writers[:slices] = - JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity), + JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * ".jld2", indices = (:, grid.Ny/2, :), schedule = TimeInterval(1minute), From 8e27b2a26200e650203f44e6fdd434e797068ada Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:04:10 -0500 Subject: [PATCH 09/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 99dbf5c3..30613e3f 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -111,7 +111,7 @@ model = HydrostaticFreeSurfaceModel(; grid, clock = Clock{eltype(grid)}(0, 0, 1), momentum_advection = CenteredSecondOrder(), tracer_advection = CenteredSecondOrder(), - buoyancy = SeawaterBuoyancy(equation_of_state=eos), + buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() From 528c1942818bfc229a4b7f11d4c5857f550a02d3 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:04:42 -0500 Subject: [PATCH 10/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 30613e3f..750d640b 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -58,7 +58,16 @@ println(underlying_grid) # We want the underwater slope to provide a depth of 500 m at y = 0 and the full 4 km by y =200. It follows # a hyperbolic tangent curve with its center at y ~= 150 at a depth of ~ -4000 + (4000 - 500)/2 = -2250 m -underwater_slope(x, y) = -1750tanh(0.00004*(y - 150000)) - 2250 +y₀ = 150kilometers +Δy = 25kilometers + +slope_depth = 500 +basin_depth = 4000 + +""" Varies smoothly between 0 when y << y₀ to 1 when y >> y₀ + Δy, over a region of width Δy. """ +step(y, y₀, Δy) = (1 + tanh((y - y₀) / Δy)) / 2 + +underwater_slope(x, y) = -slope_depth + (slope_depth - basin_depth) * tanh((y - y₀) / Δy) # TODO: add 50km troughs to the slope, dependent on x and y. Use appendix C to add detail here grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) From 9656960ce4c5300d191e94849f2a4b4985c8824d Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:05:04 -0500 Subject: [PATCH 11/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 750d640b..77fa63a6 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -72,7 +72,7 @@ underwater_slope(x, y) = -slope_depth + (slope_depth - basin_depth) * tanh((y - grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) -println(grid) +@show grid # TODO: add underwater slope at y = 0 with troughs From ba3f10a61f06335f8c36e964ec6a853504635a6e Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:06:08 -0500 Subject: [PATCH 12/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 77fa63a6..6eaa8bfe 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -118,8 +118,8 @@ vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) # Assuming no particles or biogeochemistry model = HydrostaticFreeSurfaceModel(; grid, clock = Clock{eltype(grid)}(0, 0, 1), - momentum_advection = CenteredSecondOrder(), - tracer_advection = CenteredSecondOrder(), + momentum_advection = WENO(), + tracer_advection = WENO(), buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), From 1518f01d358c951fa36a6b716b60733950073282 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:06:54 -0500 Subject: [PATCH 13/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 6eaa8bfe..bf429f9f 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -124,7 +124,7 @@ model = HydrostaticFreeSurfaceModel(; grid, coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() - closure = (horizontal_closure, vertical_closure), + closure = CATKEVerticalDiffusivity(), boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs), # NamedTuple(), tracers = (:T, :S), velocities = nothing, From 9eeebf295abf7331973efb06697c10c362e8d6e1 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:08:00 -0500 Subject: [PATCH 14/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index bf429f9f..569e5628 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -160,4 +160,4 @@ simulation.output_writers[:slices] = run!(simulation) -println(simulation) +@show simulation From fed41678d63b8c1d74a541b2362d037218935097 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:08:17 -0500 Subject: [PATCH 15/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 569e5628..78703b4a 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -54,7 +54,7 @@ underlying_grid = RectilinearGrid(arch, y = (0, Ly), z = z_faces) -println(underlying_grid) +@show underlying_grid # We want the underwater slope to provide a depth of 500 m at y = 0 and the full 4 km by y =200. It follows # a hyperbolic tangent curve with its center at y ~= 150 at a depth of ~ -4000 + (4000 - 500)/2 = -2250 m From 9fb7d9f153ba643c24119a1c17cb5376482baf9b Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:09:08 -0500 Subject: [PATCH 16/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 78703b4a..5f075ffd 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -11,7 +11,7 @@ arch = GPU() g_Earth = 9.80665 -Lx = 400000 +Lx = 400kilometers Ly = 450000 Lz = 4000 From 432b51f4ee533ba9cb43d013a26627ecf5faa161 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:09:27 -0500 Subject: [PATCH 17/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 5f075ffd..4b57156f 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -12,7 +12,7 @@ arch = GPU() g_Earth = 9.80665 Lx = 400kilometers -Ly = 450000 +Ly = 450kilometers Lz = 4000 Nx = 400 From 5ae2bb0ba8704a3f6ad9cba373f40c49bfb367e3 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:09:39 -0500 Subject: [PATCH 18/49] Update examples/antarctic_slope_current_model.jl Co-authored-by: Gregory L. Wagner --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 4b57156f..c7c6c542 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -19,7 +19,7 @@ Nx = 400 Ny = 450 Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor -sponge_width = 20000 +sponge_width = 20kilometers #= # # Setting up the grid and bathymetry: From c7ee0f713d8fd4cfe7ddec97931ca7a48eb160c1 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 10 Jul 2023 10:44:30 -0500 Subject: [PATCH 19/49] Adding some modules for new changes --- examples/antarctic_slope_current_model.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index c7c6c542..e8695815 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,13 +1,14 @@ #using CairoMakie using Oceananigans -using Oceananigans.Units: minute, minutes, hour +using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities # FIND submodule with CATKE +using Oceananigans.Units: minute, minutes, hour, kilometers using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Ian Eisenman -arch = GPU() +arch = CPU() g_Earth = 9.80665 @@ -126,15 +127,11 @@ model = HydrostaticFreeSurfaceModel(; grid, forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = CATKEVerticalDiffusivity(), boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs), # NamedTuple(), - tracers = (:T, :S), - velocities = nothing, - pressure = nothing, - diffusivity_fields = nothing, - #auxiliary_fields = nothing, # NamedTuple(), + tracers = (:T, :S) ) -println(model) - +@show model +#= println("u boundary conditions:") println(model.velocities.u.boundary_conditions) println("v boundary conditions:") @@ -154,10 +151,11 @@ filename = "asc_model_run" simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * ".jld2", - indices = (:, grid.Ny/2, :), + indices = (:, :, grid.Nz), schedule = TimeInterval(1minute), overwrite_existing = true) run!(simulation) @show simulation +=# \ No newline at end of file From 9035fcdd3040106b8ce66a0d5c83166dec425025 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 13 Jul 2023 11:54:47 -0500 Subject: [PATCH 20/49] Created low-resolution model run that generates movie of model variables at surface --- examples/antarctic_slope_current_model.jl | 187 ++++++++++++++++++---- 1 file changed, 158 insertions(+), 29 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index e8695815..8e8732ae 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,12 +1,15 @@ #using CairoMakie using Oceananigans -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities # FIND submodule with CATKE +using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule with CATKE using Oceananigans.Units: minute, minutes, hour, kilometers using SeawaterPolynomials.TEOS10 +using GLMakie +using Printf + # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the -# 2022 paper by Ian Eisenman +# 2022 paper by Si, Stewart, and Eisenman arch = CPU() @@ -16,9 +19,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 400 -Ny = 450 -Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 64 #400 +Ny = 64 #450 +Nz = 32 #70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20kilometers #= @@ -49,11 +52,12 @@ h(k) = (k - 1) / Nz z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1) underlying_grid = RectilinearGrid(arch, - size = (Nx, Ny, Nz), + size = (Nx, Ny, Nz), topology = (Periodic, Bounded, Bounded), - x = (-Lx/2, Lx/2), - y = (0, Ly), - z = z_faces) + x = (-Lx/2, Lx/2), + y = (0, Ly), + z = z_faces, + halo = (4, 4, 4)) @show underlying_grid @@ -68,7 +72,8 @@ basin_depth = 4000 """ Varies smoothly between 0 when y << y₀ to 1 when y >> y₀ + Δy, over a region of width Δy. """ step(y, y₀, Δy) = (1 + tanh((y - y₀) / Δy)) / 2 -underwater_slope(x, y) = -slope_depth + (slope_depth - basin_depth) * tanh((y - y₀) / Δy) +#underwater_slope(x, y) = -(slope_depth) + (slope_depth - basin_depth) * tanh((y - y₀) / Δy) +underwater_slope(x, y) = (-basin_depth - slope_depth + (slope_depth - basin_depth) * tanh((y - y₀) / Δy)) / 2 # TODO: add 50km troughs to the slope, dependent on x and y. Use appendix C to add detail here grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) @@ -97,19 +102,20 @@ north_sponge_layer = Relaxation(; rate=damping_rate, mask=north_mask) sponge_layers = (south_sponge_layer, north_sponge_layer) # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers - +#= Currently not used # Boundary Conditions: +# Eventually want open boundary conditions with tides, but can add that in later no_slip_bc = ValueBoundaryCondition(0.0) free_slip_bc = FluxBoundaryCondition(nothing) free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)) no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc) - +=# # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 -eos = TEOS10EquationOfState() +eos = TEOS10EquationOfState() # can compare to linear EOS later (linear not recommended for polar regions) # Coriolis Effect, using basic f-plane with precribed reference Coriolis parameter -coriolis = FPlane(f=-1.3e14) +coriolis = FPlane(latitude=-60) # Diffusivities as part of closure # TODO: make sure this works for biharmonic diffusivities as the horizontal, @@ -117,32 +123,75 @@ horizontal_closure = HorizontalScalarDiffusivity(ν=0.1, κ=0.1) vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) # Assuming no particles or biogeochemistry -model = HydrostaticFreeSurfaceModel(; grid, - clock = Clock{eltype(grid)}(0, 0, 1), +model = HydrostaticFreeSurfaceModel(; grid = grid, momentum_advection = WENO(), tracer_advection = WENO(), buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), - forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() + #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = CATKEVerticalDiffusivity(), - boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs), # NamedTuple(), - tracers = (:T, :S) + #boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs), # NamedTuple(), + tracers = (:T, :S, :e) ) -@show model +# INITIAL CONDITION FOR TEMPERATURE: a baroclinically unstable temperature distribution +""" + ramp(y, Δy) + +Linear ramp from 0 to 1 between Ly/2 - Δy/2 and Ly/2 + Δy/2. + +For example: +``` + y < Ly/2 - Δy/2 => ramp = 0 + Ly/2 - Δy/2 < y < Ly/2 + Δy/2 => ramp = y / Δy + y > Ly/2 + Δy/2 => ramp = 1 +``` +""" +ramp(y, Δy) = min(max(0, (y - Ly/2)/Δy + 1/2), 1) + +N² = 1e-5 # [s⁻²] vertical stratification +M² = 1e-7 # [s⁻²] meridional temperature gradient + +Δy = 100kilometers # width of the region of the front +ΔT = Δy * M² # temperature jump associated with the front +ϵT = 1e-2 * ΔT # noise amplitude + +Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) + ϵT * randn() + +set!(model, T=Tᵢ) + #= -println("u boundary conditions:") -println(model.velocities.u.boundary_conditions) -println("v boundary conditions:") -println(model.velocities.v.boundary_conditions) -println("w boundary conditions:") -println(model.velocities.w.boundary_conditions) +using GLMakie + +# Build coordinates with units of kilometers +x, y, z = 1e-3 .* nodes(grid, (Center(), Center(), Center())) + +T = model.tracers.T + +fig, ax, hm = heatmap(y, z, interior(T)[1, :, :], + colormap=:deep, + axis = (xlabel = "y [km]", + ylabel = "z [km]", + title = "T(x=0, y, z, t=0)", + titlesize = 24)) + +Colorbar(fig[1, 2], hm, label = "[m s⁻²]") + +current_figure() # hide +fig +=# + +# TODO: impose a temperature restoring at the surface that corresponds to Si/Stewart to help get baroclinic eddies +# Add initial condition to get turbulence running (should be doable on a laptop), then: +# Buoyancy gradient imposed at North/South boundary, surface temp as freezing temperature imposed + +@show model # # Now create a simulation and run the model # -simulation = Simulation(model; Δt=100.0, stop_time=600minutes) +simulation = Simulation(model; Δt=100.0, stop_iteration=100) # Create a NamedTuple @@ -152,10 +201,90 @@ simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * ".jld2", indices = (:, :, grid.Nz), - schedule = TimeInterval(1minute), + schedule = IterationInterval(1), overwrite_existing = true) run!(simulation) @show simulation -=# \ No newline at end of file + +# +# Make a figure and plot it +# + +filepath = filename * ".jld2" + +time_series = (u = FieldTimeSeries(filepath, "u"), + v = FieldTimeSeries(filepath, "v"), + w = FieldTimeSeries(filepath, "w"), + T = FieldTimeSeries(filepath, "T"), + S = FieldTimeSeries(filepath, "S")) + +@show time_series.u +@show time_series.w +@show time_series.T +@show time_series.S + +# Coordinate arrays +xw, yw, zw = nodes(time_series.w) +xT, yT, zT = nodes(time_series.T) + +times = time_series.w.times +intro = 1 + +n = Observable(intro) + +wₙ = @lift interior(time_series.w[$n], :, :, 1) +Tₙ = @lift interior(time_series.T[$n], :, :, 1) +Sₙ = @lift interior(time_series.S[$n], :, :, 1) +uₙ = @lift interior(time_series.u[$n], :, :, 1) +vₙ = @lift interior(time_series.v[$n], :, :, 1) + +fig = Figure(resolution = (1000, 500)) + +axis_kwargs = (xlabel="x (m)", + ylabel="y (m)", + aspect = AxisAspect(grid.Lx/grid.Ly), + limits = ((-grid.Lx/2, grid.Lx/2), (0, grid.Ly))) + +ax_w = Axis(fig[2, 1]; title = "Vertical velocity", axis_kwargs...) +ax_T = Axis(fig[2, 3]; title = "Temperature", axis_kwargs...) +#ax_S = Axis(fig[3, 1]; title = "Salinity", axis_kwargs...) +ax_v = Axis(fig[3, 1]; title = "Meridional velocity", axis_kwargs...) +ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) + +title = @lift @sprintf("t = %s", prettytime(times[$n])) + +wlims = (-3e-6, 3e-6) +Tlims = (-1.0, 1.0) +Slims = (35, 35.005) +ulims = (-0.0005, 0.0005) + + +hm_w = heatmap!(ax_w, xw, yw, wₙ; colormap = :balance, colorrange = wlims) +Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") + +hm_T = heatmap!(ax_T, xT, yT, Tₙ; colormap = :thermal)#, colorrange = Tlims) +Colorbar(fig[2, 4], hm_T; label = "ᵒC") + +#hm_S = heatmap!(ax_S, xT, yT, Sₙ; colormap = :haline)#, colorrange = Slims) +#Colorbar(fig[3, 2], hm_S; label = "g / kg") + +hm_v = heatmap!(ax_v, xw, yw, vₙ; colormap = :balance, colorrange = ulims) +Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") + +hm_u = heatmap!(ax_u, xw, yw, uₙ; colormap = :balance, colorrange = ulims) +Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") + +fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) + +current_figure() # hide +fig + +frames = intro:length(times) + +@info "Making a motion picture of ocean wind mixing and convection..." + +record(fig, filename * ".mp4", frames, framerate=8) do i + n[] = i +end \ No newline at end of file From 8a216e8b9774174b2a01b29de846c07a14c2a2a8 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Fri, 14 Jul 2023 13:52:42 -0500 Subject: [PATCH 21/49] Modified script for 5 day run --- examples/antarctic_slope_current_model.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 8e8732ae..d73b30fc 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -2,7 +2,7 @@ #using CairoMakie using Oceananigans using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule with CATKE -using Oceananigans.Units: minute, minutes, hour, kilometers +using Oceananigans.Units: minute, minutes, hour, days, kilometers using SeawaterPolynomials.TEOS10 using GLMakie @@ -191,7 +191,7 @@ fig # # Now create a simulation and run the model # -simulation = Simulation(model; Δt=100.0, stop_iteration=100) +simulation = Simulation(model; Δt=100.0, stop_time=5days) # Create a NamedTuple @@ -255,22 +255,23 @@ ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) title = @lift @sprintf("t = %s", prettytime(times[$n])) -wlims = (-3e-6, 3e-6) -Tlims = (-1.0, 1.0) +wlims = (-0.002, 0.002) +Tlims = (-0.02, 0.02) Slims = (35, 35.005) -ulims = (-0.0005, 0.0005) +ulims = (-0.08, 0.08) +vlims = (-0.08, 0.08) hm_w = heatmap!(ax_w, xw, yw, wₙ; colormap = :balance, colorrange = wlims) Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") -hm_T = heatmap!(ax_T, xT, yT, Tₙ; colormap = :thermal)#, colorrange = Tlims) +hm_T = heatmap!(ax_T, xT, yT, Tₙ; colormap = :thermal, colorrange = Tlims) Colorbar(fig[2, 4], hm_T; label = "ᵒC") #hm_S = heatmap!(ax_S, xT, yT, Sₙ; colormap = :haline)#, colorrange = Slims) #Colorbar(fig[3, 2], hm_S; label = "g / kg") -hm_v = heatmap!(ax_v, xw, yw, vₙ; colormap = :balance, colorrange = ulims) +hm_v = heatmap!(ax_v, xw, yw, vₙ; colormap = :balance, colorrange = vlims) Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") hm_u = heatmap!(ax_u, xw, yw, uₙ; colormap = :balance, colorrange = ulims) From 90bb0f4d151320499aff16fd87ba6d96ef3efb18 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 19 Jul 2023 16:04:11 -0500 Subject: [PATCH 22/49] Added wind stress as a surface boundary condition --- examples/antarctic_slope_current_model.jl | 97 ++++++++++++----------- 1 file changed, 50 insertions(+), 47 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index d73b30fc..b88824b4 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -11,7 +11,7 @@ using Printf # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman -arch = CPU() +arch = GPU() g_Earth = 9.80665 @@ -24,18 +24,10 @@ Ny = 64 #450 Nz = 32 #70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20kilometers -#= + # # Setting up the grid and bathymetry: # -# Using a linear slope that approximates the min and max spacings in the paper: -init_adjustment = (0.5*(100-10)/(Nz-1)) + (10 - (100-10)/(Nz-1)) -linear_slope(k) = (0.5*(100-10)/(Nz-1))*(k)^2 + (10 - (100-10)/(Nz-1))*(k) - init_adjustment - -spacing_adjustment = Lz / linear_slope(Nz+1) -linear_slope_z_faces(k) = -spacing_adjustment * linear_slope(k) -# Can reverse this to get grid from -4000 to 0 later -=# refinement = 20.0 # controls spacing near surface (higher means finer spaced), 12.0 stretching = 3 # controls rate of stretching at bottom, 3 @@ -93,7 +85,9 @@ u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) # Not actual forcing, just exam # of the N/S boundaries, they impose a cross-slope buoyancy gradient, and the relaxation # tiem scales decrease linearly with distance from the interior termination of the sponge layers +# # We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity +# damping_rate = 1 / 43200 # Relaxation time scale in seconds, need to make this decrease linearly toward outermost boundary south_mask = GaussianMask{:y}(center=0, width=sponge_width) north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) @@ -102,15 +96,40 @@ north_sponge_layer = Relaxation(; rate=damping_rate, mask=north_mask) sponge_layers = (south_sponge_layer, north_sponge_layer) # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers -#= Currently not used +# # Boundary Conditions: -# Eventually want open boundary conditions with tides, but can add that in later -no_slip_bc = ValueBoundaryCondition(0.0) -free_slip_bc = FluxBoundaryCondition(nothing) +# (from ocean wind mixing and convection example, +# https://clima.github.io/OceananigansDocumentation/stable/generated/ocean_wind_mixing_and_convection/) +# + +""" + boundary_ramp(y, Δy) + +Linear ramp from 1 to 0 between y = 0 and y = δy + +For example: +``` + y=0 => boundary_ramp = 1 + 0 < y < δy => boundary_ramp = (δy - y) / δy + y > δy => boundary_ramp = 0 +``` +""" +δy = 10kilometers +boundary_ramp(y, δy) = min(max(0, (δy - y)/δy), 1) + +u₁₀ = -6 # m s⁻¹, average zonal wind velocity 10 meters above the ocean at the southern boundary +v₁₀ = 6 # m s⁻¹, average meridional wind velocity 10 meters above the ocean at the southern boundary +cᴰ = 2.5e-3 # dimensionless drag coefficient +ρₐ = 1.225 # kg m⁻³, average density of air at sea-level +ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean + +# TODO: make this only apply at Southern boundary and decay to 0 elsewhere +Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) * boundary_ramp(y, δy) # m² s⁻² +Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀ * abs(v₁₀) * boundary_ramp(y, δy) # m² s⁻² + +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) +v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) -free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)) -no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc) -=# # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 eos = TEOS10EquationOfState() # can compare to linear EOS later (linear not recommended for polar regions) @@ -122,7 +141,9 @@ coriolis = FPlane(latitude=-60) horizontal_closure = HorizontalScalarDiffusivity(ν=0.1, κ=0.1) vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) +# # Assuming no particles or biogeochemistry +# model = HydrostaticFreeSurfaceModel(; grid = grid, momentum_advection = WENO(), tracer_advection = WENO(), @@ -131,7 +152,7 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = CATKEVerticalDiffusivity(), - #boundary_conditions = (u=free_slip_surface_bcs, v=free_slip_surface_bcs), # NamedTuple(), + boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) ) @@ -161,33 +182,14 @@ Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) + ϵT * randn() set!(model, T=Tᵢ) -#= -using GLMakie - -# Build coordinates with units of kilometers -x, y, z = 1e-3 .* nodes(grid, (Center(), Center(), Center())) - -T = model.tracers.T - -fig, ax, hm = heatmap(y, z, interior(T)[1, :, :], - colormap=:deep, - axis = (xlabel = "y [km]", - ylabel = "z [km]", - title = "T(x=0, y, z, t=0)", - titlesize = 24)) - -Colorbar(fig[1, 2], hm, label = "[m s⁻²]") - -current_figure() # hide -fig -=# - # TODO: impose a temperature restoring at the surface that corresponds to Si/Stewart to help get baroclinic eddies # Add initial condition to get turbulence running (should be doable on a laptop), then: # Buoyancy gradient imposed at North/South boundary, surface temp as freezing temperature imposed @show model +@show u_bcs + # # Now create a simulation and run the model # @@ -195,13 +197,13 @@ simulation = Simulation(model; Δt=100.0, stop_time=5days) # Create a NamedTuple -filename = "asc_model_run" +filename = "asc_model_run_wind10k_stress_5_days" simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * ".jld2", indices = (:, :, grid.Nz), - schedule = IterationInterval(1), + schedule = IterationInterval(10), overwrite_existing = true) run!(simulation) @@ -211,7 +213,7 @@ run!(simulation) # # Make a figure and plot it # - +#= filepath = filename * ".jld2" time_series = (u = FieldTimeSeries(filepath, "u"), @@ -262,19 +264,19 @@ ulims = (-0.08, 0.08) vlims = (-0.08, 0.08) -hm_w = heatmap!(ax_w, xw, yw, wₙ; colormap = :balance, colorrange = wlims) +hm_w = heatmap!(ax_w, xw, yw, wₙ; colormap = :balance)#, colorrange = wlims) Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") -hm_T = heatmap!(ax_T, xT, yT, Tₙ; colormap = :thermal, colorrange = Tlims) +hm_T = heatmap!(ax_T, xT, yT, Tₙ; colormap = :thermal)#, colorrange = Tlims) Colorbar(fig[2, 4], hm_T; label = "ᵒC") #hm_S = heatmap!(ax_S, xT, yT, Sₙ; colormap = :haline)#, colorrange = Slims) #Colorbar(fig[3, 2], hm_S; label = "g / kg") -hm_v = heatmap!(ax_v, xw, yw, vₙ; colormap = :balance, colorrange = vlims) +hm_v = heatmap!(ax_v, xw, yw, vₙ; colormap = :balance)#, colorrange = vlims) Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") -hm_u = heatmap!(ax_u, xw, yw, uₙ; colormap = :balance, colorrange = ulims) +hm_u = heatmap!(ax_u, xw, yw, uₙ; colormap = :balance)#, colorrange = ulims) Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) @@ -288,4 +290,5 @@ frames = intro:length(times) record(fig, filename * ".mp4", frames, framerate=8) do i n[] = i -end \ No newline at end of file +end +=# \ No newline at end of file From 9a26cb116dd7d20017ce24ea4351bb5bf810a21a Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 20 Jul 2023 22:40:18 -0500 Subject: [PATCH 23/49] Added zonal average cross-section visualization --- examples/antarctic_slope_current_model.jl | 162 +++++++++++++++++----- 1 file changed, 127 insertions(+), 35 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index b88824b4..ed76f581 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -11,7 +11,7 @@ using Printf # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman -arch = GPU() +arch = CPU() g_Earth = 9.80665 @@ -124,11 +124,11 @@ cᴰ = 2.5e-3 # dimensionless drag coefficient ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean # TODO: make this only apply at Southern boundary and decay to 0 elsewhere -Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) * boundary_ramp(y, δy) # m² s⁻² -Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀ * abs(v₁₀) * boundary_ramp(y, δy) # m² s⁻² +#Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) * boundary_ramp(y, δy) # m² s⁻² +#Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀ * abs(v₁₀) * boundary_ramp(y, δy) # m² s⁻² -u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) -v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) +#u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) +#v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 eos = TEOS10EquationOfState() # can compare to linear EOS later (linear not recommended for polar regions) @@ -152,7 +152,7 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = CATKEVerticalDiffusivity(), - boundary_conditions = (u=u_bcs, v=v_bcs), + #boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) ) @@ -193,54 +193,67 @@ set!(model, T=Tᵢ) # # Now create a simulation and run the model # -simulation = Simulation(model; Δt=100.0, stop_time=5days) +simulation = Simulation(model; Δt=100.0, stop_time=100minutes) -# Create a NamedTuple +filename = "asc_model_5_days" + +# Here we'll try also running a zonal average of the simulation: +u, v, w = model.velocities +avgT = Average(model.tracers.T, dims=1) +avgU = Average(u, dims=1) +avgV = Average(v, dims=1) +avgW = Average(w, dims=1) + +simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); + filename = filename * "_zonal_average.jld2", + schedule = IterationInterval(1), + overwrite_existing = true) -filename = "asc_model_run_wind10k_stress_5_days" simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), - filename = filename * ".jld2", + filename = filename * "_surface.jld2", indices = (:, :, grid.Nz), - schedule = IterationInterval(10), + schedule = IterationInterval(1), overwrite_existing = true) run!(simulation) +@info "Simulation completed in " * prettytime(simulation.run_wall_time) @show simulation # # Make a figure and plot it # #= -filepath = filename * ".jld2" +surface_filepath = filename * "_surface.jld2" +average_filepath = filename * "_zonal_average.jld2" -time_series = (u = FieldTimeSeries(filepath, "u"), - v = FieldTimeSeries(filepath, "v"), - w = FieldTimeSeries(filepath, "w"), - T = FieldTimeSeries(filepath, "T"), - S = FieldTimeSeries(filepath, "S")) +surface_time_series = (u = FieldTimeSeries(surface_filepath, "u"), + v = FieldTimeSeries(surface_filepath, "v"), + w = FieldTimeSeries(surface_filepath, "w"), + T = FieldTimeSeries(surface_filepath, "T"), + S = FieldTimeSeries(surface_filepath, "S")) -@show time_series.u -@show time_series.w -@show time_series.T -@show time_series.S +@show surface_time_series.u +@show surface_time_series.w +@show surface_time_series.T +@show surface_time_series.S # Coordinate arrays -xw, yw, zw = nodes(time_series.w) -xT, yT, zT = nodes(time_series.T) +srf_xw, srf_yw, srf_zw = nodes(surface_time_series.w) +srf_xT, srf_yT, srf_zT = nodes(surface_time_series.T) -times = time_series.w.times +times = surface_time_series.w.times intro = 1 n = Observable(intro) -wₙ = @lift interior(time_series.w[$n], :, :, 1) -Tₙ = @lift interior(time_series.T[$n], :, :, 1) -Sₙ = @lift interior(time_series.S[$n], :, :, 1) -uₙ = @lift interior(time_series.u[$n], :, :, 1) -vₙ = @lift interior(time_series.v[$n], :, :, 1) +srf_wₙ = @lift interior(surface_time_series.w[$n], :, :, 1) +srf_Tₙ = @lift interior(surface_time_series.T[$n], :, :, 1) +srf_Sₙ = @lift interior(surface_time_series.S[$n], :, :, 1) +srf_uₙ = @lift interior(surface_time_series.u[$n], :, :, 1) +srf_vₙ = @lift interior(surface_time_series.v[$n], :, :, 1) fig = Figure(resolution = (1000, 500)) @@ -264,19 +277,98 @@ ulims = (-0.08, 0.08) vlims = (-0.08, 0.08) -hm_w = heatmap!(ax_w, xw, yw, wₙ; colormap = :balance)#, colorrange = wlims) +hm_w = heatmap!(ax_w, srf_xw, srf_yw, srf_wₙ; colormap = :balance)#, colorrange = wlims) +Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") + +hm_T = heatmap!(ax_T, srf_xT, srf_yT, srf_Tₙ; colormap = :thermal)#, colorrange = Tlims) +Colorbar(fig[2, 4], hm_T; label = "ᵒC") + +#hm_S = heatmap!(ax_S, srf_xT, srf_yT, srf_Sₙ; colormap = :haline)#, colorrange = Slims) +#Colorbar(fig[3, 2], hm_S; label = "g / kg") + +hm_v = heatmap!(ax_v, srf_xw, srf_yw, srf_vₙ; colormap = :balance)#, colorrange = vlims) +Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") + +hm_u = heatmap!(ax_u, srf_xw, srf_yw, srf_uₙ; colormap = :balance)#, colorrange = ulims) +Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") + +fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) + +current_figure() # hide +fig + +frames = intro:length(times) + +@info "Making a motion picture of ocean wind mixing and convection..." + +record(fig, filename * "_surface.mp4", frames, framerate=8) do i + n[] = i +end + +# +# Now do all the above, but for zonal average time series data: +# + +average_time_series = (u = FieldTimeSeries(average_filepath, "u"), + v = FieldTimeSeries(average_filepath, "v"), + w = FieldTimeSeries(average_filepath, "w"), + T = FieldTimeSeries(average_filepath, "T")) + +@show average_time_series.u +@show average_time_series.w +@show average_time_series.T +@show average_time_series.v + +# Coordinate arrays +avg_xw, avg_yw, avg_zw = nodes(average_time_series.w) +avg_xT, avg_yT, avg_zT = nodes(average_time_series.T) + +times = average_time_series.w.times +intro = 1 + +n = Observable(intro) + +avg_wₙ = @lift interior(average_time_series.w[$n], 1, :, :) +avg_Tₙ = @lift interior(average_time_series.T[$n], 1, :, :) +#avg_Sₙ = @lift interior(average_time_series.S[$n], 1, :, :) +avg_uₙ = @lift interior(average_time_series.u[$n], 1, :, :) +avg_vₙ = @lift interior(average_time_series.v[$n], 1, :, :) + +fig = Figure(resolution = (1000, 500)) + +axis_kwargs = (xlabel="y (m)", + ylabel="z (m)", + aspect = AxisAspect(2.),#AxisAspect(grid.Ly/grid.Lz), + limits = ((0, grid.Ly), (-grid.Lz, 0))) + +ax_w = Axis(fig[2, 1]; title = "Vertical velocity", axis_kwargs...) +ax_T = Axis(fig[2, 3]; title = "Temperature", axis_kwargs...) +#ax_S = Axis(fig[3, 1]; title = "Salinity", axis_kwargs...) +ax_v = Axis(fig[3, 1]; title = "Meridional velocity", axis_kwargs...) +ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) + +title = @lift @sprintf("t = %s", prettytime(times[$n])) + +wlims = (-0.002, 0.002) +Tlims = (-0.02, 0.02) +Slims = (35, 35.005) +ulims = (-0.08, 0.08) +vlims = (-0.08, 0.08) + + +hm_w = heatmap!(ax_w, avg_yw, avg_zw, avg_wₙ; colormap = :balance)#, colorrange = wlims) Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") -hm_T = heatmap!(ax_T, xT, yT, Tₙ; colormap = :thermal)#, colorrange = Tlims) +hm_T = heatmap!(ax_T, avg_yT, avg_zT, avg_Tₙ; colormap = :thermal)#, colorrange = Tlims) Colorbar(fig[2, 4], hm_T; label = "ᵒC") -#hm_S = heatmap!(ax_S, xT, yT, Sₙ; colormap = :haline)#, colorrange = Slims) +#hm_S = heatmap!(ax_S, srf_xT, srf_yT, srf_Sₙ; colormap = :haline)#, colorrange = Slims) #Colorbar(fig[3, 2], hm_S; label = "g / kg") -hm_v = heatmap!(ax_v, xw, yw, vₙ; colormap = :balance)#, colorrange = vlims) +hm_v = heatmap!(ax_v, avg_yw, avg_zw, avg_vₙ; colormap = :balance)#, colorrange = vlims) Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") -hm_u = heatmap!(ax_u, xw, yw, uₙ; colormap = :balance)#, colorrange = ulims) +hm_u = heatmap!(ax_u, avg_yw, avg_zw, avg_uₙ; colormap = :balance)#, colorrange = ulims) Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) @@ -288,7 +380,7 @@ frames = intro:length(times) @info "Making a motion picture of ocean wind mixing and convection..." -record(fig, filename * ".mp4", frames, framerate=8) do i +record(fig, filename * "_average.mp4", frames, framerate=8) do i n[] = i end =# \ No newline at end of file From 788e573984b2e23a98dc5cd897357c79db130b5d Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 26 Jul 2023 11:57:53 -0400 Subject: [PATCH 24/49] Running with larger step size --- examples/antarctic_slope_current_model.jl | 46 +++++++++++------------ 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index ed76f581..10948dfb 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -188,14 +188,13 @@ set!(model, T=Tᵢ) @show model -@show u_bcs - # # Now create a simulation and run the model # -simulation = Simulation(model; Δt=100.0, stop_time=100minutes) +# Full resolution is 100 sec +simulation = Simulation(model; Δt=20minutes, stop_time=60days) -filename = "asc_model_5_days" +filename = "asc_model_60_days" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -206,7 +205,7 @@ avgW = Average(w, dims=1) simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); filename = filename * "_zonal_average.jld2", - schedule = IterationInterval(1), + schedule = IterationInterval(120), overwrite_existing = true) @@ -214,7 +213,7 @@ simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * "_surface.jld2", indices = (:, :, grid.Nz), - schedule = IterationInterval(1), + schedule = IterationInterval(120), overwrite_existing = true) run!(simulation) @@ -225,7 +224,7 @@ run!(simulation) # # Make a figure and plot it # -#= + surface_filepath = filename * "_surface.jld2" average_filepath = filename * "_zonal_average.jld2" @@ -235,10 +234,10 @@ surface_time_series = (u = FieldTimeSeries(surface_filepath, "u"), T = FieldTimeSeries(surface_filepath, "T"), S = FieldTimeSeries(surface_filepath, "S")) -@show surface_time_series.u -@show surface_time_series.w -@show surface_time_series.T -@show surface_time_series.S +#@show surface_time_series.u +#@show surface_time_series.w +#@show surface_time_series.T +#@show surface_time_series.S # Coordinate arrays srf_xw, srf_yw, srf_zw = nodes(surface_time_series.w) @@ -314,10 +313,10 @@ average_time_series = (u = FieldTimeSeries(average_filepath, "u"), w = FieldTimeSeries(average_filepath, "w"), T = FieldTimeSeries(average_filepath, "T")) -@show average_time_series.u -@show average_time_series.w -@show average_time_series.T -@show average_time_series.v +#@show average_time_series.u +#@show average_time_series.w +#@show average_time_series.T +#@show average_time_series.v # Coordinate arrays avg_xw, avg_yw, avg_zw = nodes(average_time_series.w) @@ -350,25 +349,25 @@ ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) title = @lift @sprintf("t = %s", prettytime(times[$n])) wlims = (-0.002, 0.002) -Tlims = (-0.02, 0.02) +Tlims = (-0.04, 0.04) Slims = (35, 35.005) -ulims = (-0.08, 0.08) -vlims = (-0.08, 0.08) +ulims = (-0.004, 0.004) +vlims = (-0.004, 0.004) -hm_w = heatmap!(ax_w, avg_yw, avg_zw, avg_wₙ; colormap = :balance)#, colorrange = wlims) +hm_w = heatmap!(ax_w, avg_yw, avg_zw, avg_wₙ; colormap = :balance) #, colorrange = wlims) Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") -hm_T = heatmap!(ax_T, avg_yT, avg_zT, avg_Tₙ; colormap = :thermal)#, colorrange = Tlims) +hm_T = heatmap!(ax_T, avg_yT, avg_zT, avg_Tₙ; colormap = :thermal) #, colorrange = Tlims) Colorbar(fig[2, 4], hm_T; label = "ᵒC") #hm_S = heatmap!(ax_S, srf_xT, srf_yT, srf_Sₙ; colormap = :haline)#, colorrange = Slims) #Colorbar(fig[3, 2], hm_S; label = "g / kg") -hm_v = heatmap!(ax_v, avg_yw, avg_zw, avg_vₙ; colormap = :balance)#, colorrange = vlims) +hm_v = heatmap!(ax_v, avg_yw, avg_zw, avg_vₙ; colormap = :balance) #, colorrange = vlims) Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") -hm_u = heatmap!(ax_u, avg_yw, avg_zw, avg_uₙ; colormap = :balance)#, colorrange = ulims) +hm_u = heatmap!(ax_u, avg_yw, avg_zw, avg_uₙ; colormap = :balance) #, colorrange = ulims) Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) @@ -380,7 +379,6 @@ frames = intro:length(times) @info "Making a motion picture of ocean wind mixing and convection..." -record(fig, filename * "_average.mp4", frames, framerate=8) do i +record(fig, filename * "_average_steady_bar.mp4", frames, framerate=8) do i n[] = i end -=# \ No newline at end of file From 9836e4442164f7c87a0e6af0d3126ed59cc9efb1 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 27 Jul 2023 10:44:38 -0400 Subject: [PATCH 25/49] Using smaller vertical stratification --- examples/antarctic_slope_current_model.jl | 31 ++++++++++++----------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 10948dfb..94e1d401 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -171,7 +171,7 @@ For example: """ ramp(y, Δy) = min(max(0, (y - Ly/2)/Δy + 1/2), 1) -N² = 1e-5 # [s⁻²] vertical stratification +N² = 1e-8 # [s⁻²] vertical stratification, was 1e-5 M² = 1e-7 # [s⁻²] meridional temperature gradient Δy = 100kilometers # width of the region of the front @@ -192,9 +192,9 @@ set!(model, T=Tᵢ) # Now create a simulation and run the model # # Full resolution is 100 sec -simulation = Simulation(model; Δt=20minutes, stop_time=60days) +simulation = Simulation(model; Δt=20minutes, stop_time=30days) -filename = "asc_model_60_days" +filename = "asc_model_30_days_Nsq_is_e-8" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -205,7 +205,7 @@ avgW = Average(w, dims=1) simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); filename = filename * "_zonal_average.jld2", - schedule = IterationInterval(120), + schedule = IterationInterval(10), overwrite_existing = true) @@ -213,7 +213,7 @@ simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * "_surface.jld2", indices = (:, :, grid.Nz), - schedule = IterationInterval(120), + schedule = IterationInterval(10), overwrite_existing = true) run!(simulation) @@ -224,7 +224,7 @@ run!(simulation) # # Make a figure and plot it # - +#= surface_filepath = filename * "_surface.jld2" average_filepath = filename * "_zonal_average.jld2" @@ -254,7 +254,7 @@ srf_Sₙ = @lift interior(surface_time_series.S[$n], :, :, 1) srf_uₙ = @lift interior(surface_time_series.u[$n], :, :, 1) srf_vₙ = @lift interior(surface_time_series.v[$n], :, :, 1) -fig = Figure(resolution = (1000, 500)) +fig = Figure(resolution = (4000, 2000)) axis_kwargs = (xlabel="x (m)", ylabel="y (m)", @@ -272,8 +272,8 @@ title = @lift @sprintf("t = %s", prettytime(times[$n])) wlims = (-0.002, 0.002) Tlims = (-0.02, 0.02) Slims = (35, 35.005) -ulims = (-0.08, 0.08) -vlims = (-0.08, 0.08) +ulims = (-0.12, 0.12) +vlims = (-0.12, 0.12) hm_w = heatmap!(ax_w, srf_xw, srf_yw, srf_wₙ; colormap = :balance)#, colorrange = wlims) @@ -300,7 +300,7 @@ frames = intro:length(times) @info "Making a motion picture of ocean wind mixing and convection..." -record(fig, filename * "_surface.mp4", frames, framerate=8) do i +record(fig, filename * "_surface_steady_bar.mp4", frames, framerate=8) do i n[] = i end @@ -333,7 +333,7 @@ avg_Tₙ = @lift interior(average_time_series.T[$n], 1, :, :) avg_uₙ = @lift interior(average_time_series.u[$n], 1, :, :) avg_vₙ = @lift interior(average_time_series.v[$n], 1, :, :) -fig = Figure(resolution = (1000, 500)) +fig = Figure(resolution = (4000, 2000)) axis_kwargs = (xlabel="y (m)", ylabel="z (m)", @@ -349,10 +349,10 @@ ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) title = @lift @sprintf("t = %s", prettytime(times[$n])) wlims = (-0.002, 0.002) -Tlims = (-0.04, 0.04) +Tlims = (-0.02, 0.02) Slims = (35, 35.005) -ulims = (-0.004, 0.004) -vlims = (-0.004, 0.004) +ulims = (-0.12, 0.12) +vlims = (-0.12, 0.12) hm_w = heatmap!(ax_w, avg_yw, avg_zw, avg_wₙ; colormap = :balance) #, colorrange = wlims) @@ -379,6 +379,7 @@ frames = intro:length(times) @info "Making a motion picture of ocean wind mixing and convection..." -record(fig, filename * "_average_steady_bar.mp4", frames, framerate=8) do i +record(fig, filename * "_average.mp4", frames, framerate=8) do i n[] = i end +=# \ No newline at end of file From bb40c831f39e5a9105c19902570e08087b29d23c Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 27 Jul 2023 10:46:47 -0400 Subject: [PATCH 26/49] Minor edit --- examples/antarctic_slope_current_model.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 94e1d401..111f4f7d 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -5,9 +5,6 @@ using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule using Oceananigans.Units: minute, minutes, hour, days, kilometers using SeawaterPolynomials.TEOS10 -using GLMakie -using Printf - # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman @@ -225,6 +222,9 @@ run!(simulation) # Make a figure and plot it # #= +using GLMakie +using Printf + surface_filepath = filename * "_surface.jld2" average_filepath = filename * "_zonal_average.jld2" From c70a2ffac9bbde1d22feac890e15c03ce2a13886 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 27 Jul 2023 10:47:58 -0400 Subject: [PATCH 27/49] Minor edit --- examples/antarctic_slope_current_model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 111f4f7d..53ab44c5 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -189,9 +189,9 @@ set!(model, T=Tᵢ) # Now create a simulation and run the model # # Full resolution is 100 sec -simulation = Simulation(model; Δt=20minutes, stop_time=30days) +simulation = Simulation(model; Δt=20minutes, stop_time=60days) -filename = "asc_model_30_days_Nsq_is_e-8" +filename = "asc_model_60_days_Nsq_is_e-8" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities From 820941e0462e5bf341c574d73871a3b3cd0d8912 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 27 Jul 2023 10:51:13 -0400 Subject: [PATCH 28/49] Minor edit --- examples/antarctic_slope_current_model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 53ab44c5..5eb09e19 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -202,7 +202,7 @@ avgW = Average(w, dims=1) simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); filename = filename * "_zonal_average.jld2", - schedule = IterationInterval(10), + schedule = IterationInterval(20), overwrite_existing = true) @@ -210,7 +210,7 @@ simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), filename = filename * "_surface.jld2", indices = (:, :, grid.Nz), - schedule = IterationInterval(10), + schedule = IterationInterval(20), overwrite_existing = true) run!(simulation) From 16500721e5fd40e3eda4a4720cee151fc17a152d Mon Sep 17 00:00:00 2001 From: jlk9 Date: Fri, 28 Jul 2023 16:23:48 -0400 Subject: [PATCH 29/49] Running more tests --- examples/antarctic_slope_current_model.jl | 35 ++++++++++++----------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 5eb09e19..eb89bef9 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -8,7 +8,7 @@ using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman -arch = CPU() +arch = GPU() g_Earth = 9.80665 @@ -16,9 +16,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 64 #400 -Ny = 64 #450 -Nz = 32 #70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 400 #64 #400 +Ny = 450 #64 #450 +Nz = 70 #8 #70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20kilometers @@ -114,24 +114,25 @@ For example: δy = 10kilometers boundary_ramp(y, δy) = min(max(0, (δy - y)/δy), 1) -u₁₀ = -6 # m s⁻¹, average zonal wind velocity 10 meters above the ocean at the southern boundary -v₁₀ = 6 # m s⁻¹, average meridional wind velocity 10 meters above the ocean at the southern boundary +u₁₀(y) = -6 * (Ly - y) / Ly # m s⁻¹, average zonal wind velocity 10 meters above the ocean at the southern boundary +v₁₀(y) = 6 * (Ly - y) / Ly # m s⁻¹, average meridional wind velocity 10 meters above the ocean at the southern boundary + cᴰ = 2.5e-3 # dimensionless drag coefficient ρₐ = 1.225 # kg m⁻³, average density of air at sea-level ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean # TODO: make this only apply at Southern boundary and decay to 0 elsewhere -#Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) * boundary_ramp(y, δy) # m² s⁻² -#Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀ * abs(v₁₀) * boundary_ramp(y, δy) # m² s⁻² +Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀(x) * abs(u₁₀(y)) # m² s⁻² +Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀(x) * abs(v₁₀(y)) # m² s⁻² -#u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) -#v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) +v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 eos = TEOS10EquationOfState() # can compare to linear EOS later (linear not recommended for polar regions) # Coriolis Effect, using basic f-plane with precribed reference Coriolis parameter -coriolis = FPlane(latitude=-60) +coriolis = BetaPlane(f₀=1.3e-4, β=1e-11) #FPlane(latitude=-60) # Diffusivities as part of closure # TODO: make sure this works for biharmonic diffusivities as the horizontal, @@ -141,7 +142,7 @@ vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) # # Assuming no particles or biogeochemistry # -model = HydrostaticFreeSurfaceModel(; grid = grid, +model = HydrostaticFreeSurfaceModel(; grid = underlying_grid, momentum_advection = WENO(), tracer_advection = WENO(), buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), @@ -189,9 +190,9 @@ set!(model, T=Tᵢ) # Now create a simulation and run the model # # Full resolution is 100 sec -simulation = Simulation(model; Δt=20minutes, stop_time=60days) +simulation = Simulation(model; Δt=100.0, stop_time=60days) -filename = "asc_model_60_days_Nsq_is_e-8" +filename = "asc_model_60_days_Nsq_is_e-8_no_slope_hi_res_custom_beta_plane" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -217,11 +218,11 @@ run!(simulation) @info "Simulation completed in " * prettytime(simulation.run_wall_time) @show simulation - +#= # # Make a figure and plot it # -#= + using GLMakie using Printf @@ -300,7 +301,7 @@ frames = intro:length(times) @info "Making a motion picture of ocean wind mixing and convection..." -record(fig, filename * "_surface_steady_bar.mp4", frames, framerate=8) do i +record(fig, filename * "_surface.mp4", frames, framerate=8) do i n[] = i end From 0b43b1ec52e44f2c80c33eaa6394b88f55481a4a Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 2 Aug 2023 12:56:46 -0500 Subject: [PATCH 30/49] Trying hi-res model with custom beta plane --- examples/antarctic_slope_current_model.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index eb89bef9..1f8498f9 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -16,9 +16,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 400 #64 #400 -Ny = 450 #64 #450 -Nz = 70 #8 #70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 400 +Ny = 450 +Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20kilometers @@ -169,7 +169,7 @@ For example: """ ramp(y, Δy) = min(max(0, (y - Ly/2)/Δy + 1/2), 1) -N² = 1e-8 # [s⁻²] vertical stratification, was 1e-5 +N² = 1e-6 # [s⁻²] vertical stratification, was 1e-5 M² = 1e-7 # [s⁻²] meridional temperature gradient Δy = 100kilometers # width of the region of the front @@ -192,7 +192,7 @@ set!(model, T=Tᵢ) # Full resolution is 100 sec simulation = Simulation(model; Δt=100.0, stop_time=60days) -filename = "asc_model_60_days_Nsq_is_e-8_no_slope_hi_res_custom_beta_plane" +filename = "asc_model_60_days_no_slope_hi_res_custom_beta_plane" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities From 3e25a6a2669277a35260f75395699e6ef5b7bf18 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 15 Aug 2023 13:14:38 -0500 Subject: [PATCH 31/49] Testing different runs --- examples/antarctic_slope_current_model.jl | 25 +++++++++++------------ 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 1f8498f9..43519c65 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -8,7 +8,7 @@ using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman -arch = GPU() +arch = CPU() g_Earth = 9.80665 @@ -16,9 +16,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 400 -Ny = 450 -Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 64 +Ny = 64 +Nz = 32 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20kilometers @@ -122,8 +122,8 @@ cᴰ = 2.5e-3 # dimensionless drag coefficient ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean # TODO: make this only apply at Southern boundary and decay to 0 elsewhere -Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀(x) * abs(u₁₀(y)) # m² s⁻² -Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀(x) * abs(v₁₀(y)) # m² s⁻² +Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀(y) * abs(u₁₀(y)) # m² s⁻² +Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀(y) * abs(v₁₀(y)) # m² s⁻² u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) @@ -142,7 +142,7 @@ vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) # # Assuming no particles or biogeochemistry # -model = HydrostaticFreeSurfaceModel(; grid = underlying_grid, +model = HydrostaticFreeSurfaceModel(; grid = grid, momentum_advection = WENO(), tracer_advection = WENO(), buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), @@ -150,7 +150,7 @@ model = HydrostaticFreeSurfaceModel(; grid = underlying_grid, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = CATKEVerticalDiffusivity(), - #boundary_conditions = (u=u_bcs, v=v_bcs), + boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) ) @@ -169,7 +169,7 @@ For example: """ ramp(y, Δy) = min(max(0, (y - Ly/2)/Δy + 1/2), 1) -N² = 1e-6 # [s⁻²] vertical stratification, was 1e-5 +N² = 1e-8 # [s⁻²] vertical stratification, was 1e-5 M² = 1e-7 # [s⁻²] meridional temperature gradient Δy = 100kilometers # width of the region of the front @@ -190,9 +190,9 @@ set!(model, T=Tᵢ) # Now create a simulation and run the model # # Full resolution is 100 sec -simulation = Simulation(model; Δt=100.0, stop_time=60days) +simulation = Simulation(model; Δt=20minutes, stop_time=60days) -filename = "asc_model_60_days_no_slope_hi_res_custom_beta_plane" +filename = "asc_model_60_days_Nsq_neg8_custom_beta_plane_wind_stress" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -218,7 +218,7 @@ run!(simulation) @info "Simulation completed in " * prettytime(simulation.run_wall_time) @show simulation -#= + # # Make a figure and plot it # @@ -383,4 +383,3 @@ frames = intro:length(times) record(fig, filename * "_average.mp4", frames, framerate=8) do i n[] = i end -=# \ No newline at end of file From 380708d7e19853a339869972006bfe2508750623 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Mon, 21 Aug 2023 14:08:19 -0500 Subject: [PATCH 32/49] Added wind stress and relaxation sponge layers --- examples/antarctic_slope_current_model.jl | 25 ++++++++++++----------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 43519c65..d2e22386 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -76,7 +76,6 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) # # Forcings: -u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) # Not actual forcing, just example # TODO: need to use forcings to enact the sponge layers. Their support is within 20000 m # of the N/S boundaries, they impose a cross-slope buoyancy gradient, and the relaxation @@ -86,10 +85,12 @@ u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) # Not actual forcing, just exam # We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity # damping_rate = 1 / 43200 # Relaxation time scale in seconds, need to make this decrease linearly toward outermost boundary -south_mask = GaussianMask{:y}(center=0, width=sponge_width) -north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) +#south_mask = GaussianMask{:y}(center=0, width=sponge_width) +#north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) +south_mask(x,y,z) = y < sponge_width +north_mask(x,y,z) = y > Ly - sponge_width south_sponge_layer = Relaxation(; rate=damping_rate, mask=south_mask) -north_sponge_layer = Relaxation(; rate=damping_rate, mask=north_mask) +north_sponge_layer = Relaxation(; rate=damping_rate, target=0.1, mask=north_mask)#, target=LinearTarget{:z}(intercept=1.0, gradient=0)) sponge_layers = (south_sponge_layer, north_sponge_layer) # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers @@ -148,7 +149,7 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), - #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() + forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() closure = CATKEVerticalDiffusivity(), boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) @@ -169,8 +170,8 @@ For example: """ ramp(y, Δy) = min(max(0, (y - Ly/2)/Δy + 1/2), 1) -N² = 1e-8 # [s⁻²] vertical stratification, was 1e-5 -M² = 1e-7 # [s⁻²] meridional temperature gradient +N² = 1e-7 # [s⁻²] vertical stratification, was 1e-5 +M² = 1e-5 # [s⁻²] meridional temperature gradient, was 1e-7 Δy = 100kilometers # width of the region of the front ΔT = Δy * M² # temperature jump associated with the front @@ -192,7 +193,7 @@ set!(model, T=Tᵢ) # Full resolution is 100 sec simulation = Simulation(model; Δt=20minutes, stop_time=60days) -filename = "asc_model_60_days_Nsq_neg8_custom_beta_plane_wind_stress" +filename = "asc_model_60_days_Msq_neg5_custom_beta_plane_wind_stress_sponge" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -270,11 +271,11 @@ ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) title = @lift @sprintf("t = %s", prettytime(times[$n])) -wlims = (-0.002, 0.002) -Tlims = (-0.02, 0.02) +wlims = (-0.001, 0.001) +Tlims = (-1.1, 1.1) Slims = (35, 35.005) -ulims = (-0.12, 0.12) -vlims = (-0.12, 0.12) +ulims = (-1.2, 1.2) +vlims = (-0.2, 0.2) hm_w = heatmap!(ax_w, srf_xw, srf_yw, srf_wₙ; colormap = :balance)#, colorrange = wlims) From 12a133ea66dbdf11a58c85ae473de01f04ebd2a4 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 23 Aug 2023 11:54:44 -0500 Subject: [PATCH 33/49] Added speed data and plots --- examples/antarctic_slope_current_model.jl | 112 ++++++++++++++++++---- 1 file changed, 91 insertions(+), 21 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index d2e22386..fc7b5bd3 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,6 +1,7 @@ #using CairoMakie using Oceananigans +using Oceananigans.Fields, Oceananigans.AbstractOperations using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule with CATKE using Oceananigans.Units: minute, minutes, hour, days, kilometers using SeawaterPolynomials.TEOS10 @@ -8,7 +9,7 @@ using SeawaterPolynomials.TEOS10 # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman -arch = CPU() +arch = GPU() g_Earth = 9.80665 @@ -16,9 +17,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 64 -Ny = 64 -Nz = 32 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 400 +Ny = 450 +Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor sponge_width = 20kilometers @@ -84,14 +85,19 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) # # We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity # +#= damping_rate = 1 / 43200 # Relaxation time scale in seconds, need to make this decrease linearly toward outermost boundary #south_mask = GaussianMask{:y}(center=0, width=sponge_width) #north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) south_mask(x,y,z) = y < sponge_width -north_mask(x,y,z) = y > Ly - sponge_width +north_mask(x,y,z) = y > (Ly - sponge_width) south_sponge_layer = Relaxation(; rate=damping_rate, mask=south_mask) -north_sponge_layer = Relaxation(; rate=damping_rate, target=0.1, mask=north_mask)#, target=LinearTarget{:z}(intercept=1.0, gradient=0)) +north_sponge_layer = Relaxation(; rate=damping_rate, mask=north_mask, target=0.0)#, target=LinearTarget{:z}(intercept=1.0, gradient=0)) sponge_layers = (south_sponge_layer, north_sponge_layer) + +north_sponge_layer_T = Relaxation(; rate=damping_rate, mask=north_mask, target=1.0) +sponge_layers_T = (south_sponge_layer, north_sponge_layer_T) +=# # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers # @@ -149,7 +155,7 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), - forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers, S=sponge_layers), # NamedTuple() + #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers_T, S=sponge_layers), closure = CATKEVerticalDiffusivity(), boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) @@ -193,7 +199,7 @@ set!(model, T=Tᵢ) # Full resolution is 100 sec simulation = Simulation(model; Δt=20minutes, stop_time=60days) -filename = "asc_model_60_days_Msq_neg5_custom_beta_plane_wind_stress_sponge" +filename = "asc_model_hi_res_60_days_Msq_neg5_custom_beta_plane_wind_stress" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -202,28 +208,38 @@ avgU = Average(u, dims=1) avgV = Average(v, dims=1) avgW = Average(w, dims=1) +# Here we'll compute and record vorticity and speed: +ω = ∂x(v) - ∂y(u) +s = sqrt(u^2 + v^2 + w^2) + +iter_interval = 20 + simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); filename = filename * "_zonal_average.jld2", - schedule = IterationInterval(20), + schedule = IterationInterval(iter_interval), overwrite_existing = true) -simulation.output_writers[:slices] = - JLD2OutputWriter(model, merge(model.velocities, model.tracers), - filename = filename * "_surface.jld2", - indices = (:, :, grid.Nz), - schedule = IterationInterval(20), - overwrite_existing = true) +simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), + filename = filename * "_surface.jld2", + indices = (:, :, grid.Nz), + schedule = IterationInterval(iter_interval), + overwrite_existing = true) -run!(simulation) +simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ω, s), + filename = filename * "_fields.jld2", + indices = (:, :, grid.Nz), + schedule = IterationInterval(iter_interval), + overwrite_existing = true) +run!(simulation) @info "Simulation completed in " * prettytime(simulation.run_wall_time) @show simulation # # Make a figure and plot it # - +#= using GLMakie using Printf @@ -278,19 +294,19 @@ ulims = (-1.2, 1.2) vlims = (-0.2, 0.2) -hm_w = heatmap!(ax_w, srf_xw, srf_yw, srf_wₙ; colormap = :balance)#, colorrange = wlims) +hm_w = heatmap!(ax_w, srf_xw, srf_yw, srf_wₙ; colormap = :balance, colorrange = wlims) Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") -hm_T = heatmap!(ax_T, srf_xT, srf_yT, srf_Tₙ; colormap = :thermal)#, colorrange = Tlims) +hm_T = heatmap!(ax_T, srf_xT, srf_yT, srf_Tₙ; colormap = :thermal, colorrange = Tlims) Colorbar(fig[2, 4], hm_T; label = "ᵒC") #hm_S = heatmap!(ax_S, srf_xT, srf_yT, srf_Sₙ; colormap = :haline)#, colorrange = Slims) #Colorbar(fig[3, 2], hm_S; label = "g / kg") -hm_v = heatmap!(ax_v, srf_xw, srf_yw, srf_vₙ; colormap = :balance)#, colorrange = vlims) +hm_v = heatmap!(ax_v, srf_xw, srf_yw, srf_vₙ; colormap = :balance, colorrange = vlims) Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") -hm_u = heatmap!(ax_u, srf_xw, srf_yw, srf_uₙ; colormap = :balance)#, colorrange = ulims) +hm_u = heatmap!(ax_u, srf_xw, srf_yw, srf_uₙ; colormap = :balance, colorrange = ulims) Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) @@ -306,6 +322,59 @@ record(fig, filename * "_surface.mp4", frames, framerate=8) do i n[] = i end +# +# Here we'll plot vorticity and speed: +# + +ω_timeseries = FieldTimeSeries(filename * "_fields.jld2", "ω") +s_timeseries = FieldTimeSeries(filename * "_fields.jld2", "s") + +@show ω_timeseries +@show s_timeseries + +xω, yω, zω = nodes(ω_timeseries) +xs, ys, zs = nodes(s_timeseries) +nothing # hide + +set_theme!(Theme(fontsize = 24)) + +@info "Making a neat movie of vorticity and speed..." + +fig = Figure(resolution = (4000, 2000)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + limits = ((-Lx/2, Lx/2), (0, Ly)), + aspect = AxisAspect(1)) + +ax_ω = Axis(fig[2, 1]; title = "Vorticity", axis_kwargs...) +ax_s = Axis(fig[2, 3]; title = "Speed", axis_kwargs...) + +n = Observable(intro) + +ω = @lift interior(ω_timeseries[$n], :, :, 1) +s = @lift interior(s_timeseries[$n], :, :, 1) + +hm_ω = heatmap!(ax_ω, xω, yω, ω; colormap = :balance, colorrange = (-2e-4, 2e-4)) +Colorbar(fig[2, 2], hm_ω; label = "m s⁻²") +hm_s = heatmap!(ax_s, xs, ys, s; colormap = :speed, colorrange = (0, 1.5)) +Colorbar(fig[2, 4], hm_s; label = "m s⁻¹") + +title = @lift @sprintf("t = %s", prettytime(times[$n])) +Label(fig[1, 1:2], title, fontsize=24, tellwidth=false) + +current_figure() # hide +fig + +frames = intro:length(times) + +@info "Making a neat animation of vorticity and speed..." + +record(fig, filename * "_fields.mp4", frames, framerate=8) do i + n[] = i +end +=# +#= # # Now do all the above, but for zonal average time series data: # @@ -384,3 +453,4 @@ frames = intro:length(times) record(fig, filename * "_average.mp4", frames, framerate=8) do i n[] = i end +=# \ No newline at end of file From cf94e3ebaec91241e8db8a5a642ccf999b590d7b Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 23 Aug 2023 11:57:36 -0500 Subject: [PATCH 34/49] minor adds --- examples/antarctic_slope_current_model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index fc7b5bd3..a774667b 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -197,7 +197,7 @@ set!(model, T=Tᵢ) # Now create a simulation and run the model # # Full resolution is 100 sec -simulation = Simulation(model; Δt=20minutes, stop_time=60days) +simulation = Simulation(model; Δt=100.0, stop_time=60days) filename = "asc_model_hi_res_60_days_Msq_neg5_custom_beta_plane_wind_stress" @@ -212,7 +212,7 @@ avgW = Average(w, dims=1) ω = ∂x(v) - ∂y(u) s = sqrt(u^2 + v^2 + w^2) -iter_interval = 20 +iter_interval = 240 simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); filename = filename * "_zonal_average.jld2", From 31d174824b45e60b69b2d1592e6dd75e55cc65af Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 24 Aug 2023 12:56:00 -0500 Subject: [PATCH 35/49] Fixed boundary condition implementation --- examples/antarctic_slope_current_model.jl | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index a774667b..dc4324bb 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -121,19 +121,20 @@ For example: δy = 10kilometers boundary_ramp(y, δy) = min(max(0, (δy - y)/δy), 1) -u₁₀(y) = -6 * (Ly - y) / Ly # m s⁻¹, average zonal wind velocity 10 meters above the ocean at the southern boundary -v₁₀(y) = 6 * (Ly - y) / Ly # m s⁻¹, average meridional wind velocity 10 meters above the ocean at the southern boundary + +@inline u₁₀(y, p) = -6 * (p.Ly - y) / p.Ly # m s⁻¹, average zonal wind velocity 10 meters above the ocean at the southern boundary +@inline v₁₀(y, p) = 6 * (p.Ly - y) / p.Ly # m s⁻¹, average meridional wind velocity 10 meters above the ocean at the southern boundary cᴰ = 2.5e-3 # dimensionless drag coefficient ρₐ = 1.225 # kg m⁻³, average density of air at sea-level ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean # TODO: make this only apply at Southern boundary and decay to 0 elsewhere -Qᵘ(x, y, z) = - ρₐ / ρₒ * cᴰ * u₁₀(y) * abs(u₁₀(y)) # m² s⁻² -Qᵛ(x, y, z) = - ρₐ / ρₒ * cᴰ * v₁₀(y) * abs(v₁₀(y)) # m² s⁻² +@inline Qᵘ(x, y, t, p) = - p.ρₐ / p.ρₒ * p.cᴰ * u₁₀(y, p) * abs(u₁₀(y, p)) # m² s⁻² +@inline Qᵛ(x, y, t, p) = - p.ρₐ / p.ρₒ * p.cᴰ * v₁₀(y, p) * abs(v₁₀(y, p)) # m² s⁻² -u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) -v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ)) +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ, parameters=(; cᴰ, ρₐ, ρₒ, Ly))) +v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ, parameters=(; cᴰ, ρₐ, ρₒ, Ly))) # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 eos = TEOS10EquationOfState() # can compare to linear EOS later (linear not recommended for polar regions) @@ -183,7 +184,7 @@ M² = 1e-5 # [s⁻²] meridional temperature gradient, was 1e-7 ΔT = Δy * M² # temperature jump associated with the front ϵT = 1e-2 * ΔT # noise amplitude -Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) + ϵT * randn() +Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) #+ ϵT * randn() set!(model, T=Tᵢ) @@ -199,7 +200,7 @@ set!(model, T=Tᵢ) # Full resolution is 100 sec simulation = Simulation(model; Δt=100.0, stop_time=60days) -filename = "asc_model_hi_res_60_days_Msq_neg5_custom_beta_plane_wind_stress" +filename = "asc_model_lo_res_60_days_wind_stress_no_Tnoise" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -235,11 +236,11 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ω, s), run!(simulation) @info "Simulation completed in " * prettytime(simulation.run_wall_time) @show simulation - +#= # # Make a figure and plot it # -#= + using GLMakie using Printf From aaf894c5ec127eedaba29df0491e91d9c48a8929 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Thu, 24 Aug 2023 12:57:10 -0500 Subject: [PATCH 36/49] Changed filename --- examples/antarctic_slope_current_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index dc4324bb..3efcd418 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -200,7 +200,7 @@ set!(model, T=Tᵢ) # Full resolution is 100 sec simulation = Simulation(model; Δt=100.0, stop_time=60days) -filename = "asc_model_lo_res_60_days_wind_stress_no_Tnoise" +filename = "asc_model_hi_res_60_days_wind_stress_no_Tnoise" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities From 5cec2816a08997c4da1822757b2cd66c87ec64ad Mon Sep 17 00:00:00 2001 From: Joseph Kump Date: Fri, 20 Oct 2023 12:27:53 -0500 Subject: [PATCH 37/49] Can now run full workflow remotely --- examples/antarctic_slope_current_model.jl | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 3efcd418..186bd24f 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,5 +1,5 @@ -#using CairoMakie +using CairoMakie using Oceananigans using Oceananigans.Fields, Oceananigans.AbstractOperations using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule with CATKE @@ -177,14 +177,14 @@ For example: """ ramp(y, Δy) = min(max(0, (y - Ly/2)/Δy + 1/2), 1) -N² = 1e-7 # [s⁻²] vertical stratification, was 1e-5 +N² = 1e-8 # [s⁻²] vertical stratification, was 1e-5 M² = 1e-5 # [s⁻²] meridional temperature gradient, was 1e-7 Δy = 100kilometers # width of the region of the front ΔT = Δy * M² # temperature jump associated with the front ϵT = 1e-2 * ΔT # noise amplitude -Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) #+ ϵT * randn() +Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) + ϵT * randn() set!(model, T=Tᵢ) @@ -200,7 +200,7 @@ set!(model, T=Tᵢ) # Full resolution is 100 sec simulation = Simulation(model; Δt=100.0, stop_time=60days) -filename = "asc_model_hi_res_60_days_wind_stress_no_Tnoise" +filename = "asc_model_hi_res_60_days_wind_stress_Nsq_neg8" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities @@ -215,10 +215,10 @@ s = sqrt(u^2 + v^2 + w^2) iter_interval = 240 -simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); - filename = filename * "_zonal_average.jld2", - schedule = IterationInterval(iter_interval), - overwrite_existing = true) +#simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); +# filename = filename * "_zonal_average.jld2", +# schedule = IterationInterval(iter_interval), +# overwrite_existing = true) simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), @@ -236,12 +236,11 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ω, s), run!(simulation) @info "Simulation completed in " * prettytime(simulation.run_wall_time) @show simulation -#= + # # Make a figure and plot it # -using GLMakie using Printf surface_filepath = filename * "_surface.jld2" @@ -374,7 +373,7 @@ frames = intro:length(times) record(fig, filename * "_fields.mp4", frames, framerate=8) do i n[] = i end -=# + #= # # Now do all the above, but for zonal average time series data: From 35181a6dd621232f9ea5c242e8ac84f5350896a9 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Fri, 20 Oct 2023 13:33:03 -0500 Subject: [PATCH 38/49] Trying relaxation on boundaries again --- examples/antarctic_slope_current_model.jl | 24 +++++++---------------- 1 file changed, 7 insertions(+), 17 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 186bd24f..da3fc24d 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -1,5 +1,6 @@ using CairoMakie +using Printf using Oceananigans using Oceananigans.Fields, Oceananigans.AbstractOperations using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule with CATKE @@ -62,7 +63,6 @@ basin_depth = 4000 """ Varies smoothly between 0 when y << y₀ to 1 when y >> y₀ + Δy, over a region of width Δy. """ step(y, y₀, Δy) = (1 + tanh((y - y₀) / Δy)) / 2 -#underwater_slope(x, y) = -(slope_depth) + (slope_depth - basin_depth) * tanh((y - y₀) / Δy) underwater_slope(x, y) = (-basin_depth - slope_depth + (slope_depth - basin_depth) * tanh((y - y₀) / Δy)) / 2 # TODO: add 50km troughs to the slope, dependent on x and y. Use appendix C to add detail here @@ -76,16 +76,12 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) # Creating the Hydrostatic model: # -# Forcings: - -# TODO: need to use forcings to enact the sponge layers. Their support is within 20000 m -# of the N/S boundaries, they impose a cross-slope buoyancy gradient, and the relaxation -# tiem scales decrease linearly with distance from the interior termination of the sponge layers +# Forcings (ignored for now): # # We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity # -#= + damping_rate = 1 / 43200 # Relaxation time scale in seconds, need to make this decrease linearly toward outermost boundary #south_mask = GaussianMask{:y}(center=0, width=sponge_width) #north_mask = GaussianMask{:y}(center=Ly, width=sponge_width) @@ -97,13 +93,11 @@ sponge_layers = (south_sponge_layer, north_sponge_layer) north_sponge_layer_T = Relaxation(; rate=damping_rate, mask=north_mask, target=1.0) sponge_layers_T = (south_sponge_layer, north_sponge_layer_T) -=# + # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers # # Boundary Conditions: -# (from ocean wind mixing and convection example, -# https://clima.github.io/OceananigansDocumentation/stable/generated/ocean_wind_mixing_and_convection/) # """ @@ -156,7 +150,7 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), - #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers_T, S=sponge_layers), + forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers_T, S=sponge_layers), closure = CATKEVerticalDiffusivity(), boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) @@ -197,6 +191,7 @@ set!(model, T=Tᵢ) # # Now create a simulation and run the model # + # Full resolution is 100 sec simulation = Simulation(model; Δt=100.0, stop_time=60days) @@ -241,8 +236,6 @@ run!(simulation) # Make a figure and plot it # -using Printf - surface_filepath = filename * "_surface.jld2" average_filepath = filename * "_zonal_average.jld2" @@ -253,9 +246,6 @@ surface_time_series = (u = FieldTimeSeries(surface_filepath, "u"), S = FieldTimeSeries(surface_filepath, "S")) #@show surface_time_series.u -#@show surface_time_series.w -#@show surface_time_series.T -#@show surface_time_series.S # Coordinate arrays srf_xw, srf_yw, srf_zw = nodes(surface_time_series.w) @@ -287,13 +277,13 @@ ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) title = @lift @sprintf("t = %s", prettytime(times[$n])) +# These limits are hardcoded right now, can be changed wlims = (-0.001, 0.001) Tlims = (-1.1, 1.1) Slims = (35, 35.005) ulims = (-1.2, 1.2) vlims = (-0.2, 0.2) - hm_w = heatmap!(ax_w, srf_xw, srf_yw, srf_wₙ; colormap = :balance, colorrange = wlims) Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") From cc71c6fe84eb3ac33c3c7edb0a50f27ef7e99d42 Mon Sep 17 00:00:00 2001 From: Joseph Kump Date: Mon, 23 Oct 2023 11:47:33 -0500 Subject: [PATCH 39/49] Removed relaxation forcing --- examples/antarctic_slope_current_model.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index da3fc24d..c17116b1 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -77,7 +77,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) # # Forcings (ignored for now): - +#= # # We'll use Relaxation() to impose a sponge layer forcing on velocity, temperature, and salinity # @@ -93,7 +93,7 @@ sponge_layers = (south_sponge_layer, north_sponge_layer) north_sponge_layer_T = Relaxation(; rate=damping_rate, mask=north_mask, target=1.0) sponge_layers_T = (south_sponge_layer, north_sponge_layer_T) - +=# # TODO: compose north_mask and south_mask together into one sponge layer, OR compose north/south sponge layers # @@ -150,7 +150,7 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), coriolis = coriolis, free_surface = ImplicitFreeSurface(gravitational_acceleration=g_Earth), - forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers_T, S=sponge_layers), + #forcing = (u=sponge_layers, v=sponge_layers, w=sponge_layers, T=sponge_layers_T, S=sponge_layers), closure = CATKEVerticalDiffusivity(), boundary_conditions = (u=u_bcs, v=v_bcs), tracers = (:T, :S, :e) @@ -193,9 +193,9 @@ set!(model, T=Tᵢ) # # Full resolution is 100 sec -simulation = Simulation(model; Δt=100.0, stop_time=60days) +simulation = Simulation(model; Δt=100.0, stop_time=2days) -filename = "asc_model_hi_res_60_days_wind_stress_Nsq_neg8" +filename = "asc_model_hi_res_2_days_wind_stress_Nsq_neg8" # Here we'll try also running a zonal average of the simulation: u, v, w = model.velocities From 8fe96d1a70f273f7df53a7714300065201f4f296 Mon Sep 17 00:00:00 2001 From: Joseph Kump Date: Tue, 31 Oct 2023 14:19:51 -0500 Subject: [PATCH 40/49] Add sea ice package --- examples/antarctic_slope_current_model.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index c17116b1..67a072d0 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -7,6 +7,8 @@ using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule using Oceananigans.Units: minute, minutes, hour, days, kilometers using SeawaterPolynomials.TEOS10 +using ClimaSeaIce + # This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the # 2022 paper by Si, Stewart, and Eisenman From 9d590d8556fb17dc1be50520b727ce799e4680e5 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 31 Oct 2023 16:38:26 -0500 Subject: [PATCH 41/49] Adding architecture field to SlabIceModel --- src/SlabSeaIceModels/slab_sea_ice_model.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/SlabSeaIceModels/slab_sea_ice_model.jl b/src/SlabSeaIceModels/slab_sea_ice_model.jl index c5fdb254..00f6c5a1 100644 --- a/src/SlabSeaIceModels/slab_sea_ice_model.jl +++ b/src/SlabSeaIceModels/slab_sea_ice_model.jl @@ -46,6 +46,7 @@ struct SlabSeaIceModel{GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} < ice_consolidation_thickness :: MIT # Numerics advection :: A + architecture end const SSIM = SlabSeaIceModel @@ -156,6 +157,8 @@ function SlabSeaIceModel(grid; heat_boundary_conditions = (top = top_heat_boundary_condition, bottom = bottom_heat_boundary_condition) + architecture = architecture(grid) + return SlabSeaIceModel(grid, clock, timestepper, @@ -169,7 +172,8 @@ function SlabSeaIceModel(grid; internal_heat_flux_function, phase_transitions, ice_consolidation_thickness, - advection) + advection, + architecture) end function set!(model::SSIM; h=nothing, α=nothing) From fc13b51489d5ca28db473e6b05fa3a43b8ea7c3d Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 31 Oct 2023 16:48:23 -0500 Subject: [PATCH 42/49] Cleaning up architecture field --- src/SlabSeaIceModels/slab_sea_ice_model.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/SlabSeaIceModels/slab_sea_ice_model.jl b/src/SlabSeaIceModels/slab_sea_ice_model.jl index 00f6c5a1..e9425ef9 100644 --- a/src/SlabSeaIceModels/slab_sea_ice_model.jl +++ b/src/SlabSeaIceModels/slab_sea_ice_model.jl @@ -15,6 +15,8 @@ using Oceananigans.Fields: field, Field, Center, ZeroField, ConstantField # Simulations interface import Oceananigans: fields, prognostic_fields +import Oceananigans.Architectures: AbstractArchitecture +import Oceananigans.Grids: architecture import Oceananigans.Fields: set! import Oceananigans.Models: AbstractModel import Oceananigans.OutputWriters: default_included_properties @@ -26,7 +28,8 @@ import Oceananigans.Utils: prettytime # import Oceananigans.Fields: field # field(loc, a::Number, grid) = ConstantField(a) -struct SlabSeaIceModel{GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} <: AbstractModel{TS} +struct SlabSeaIceModel{AR<:AbstractArchitecture, GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} <: AbstractModel{TS} + architecture :: AR grid :: GR clock :: CL timestepper :: TS @@ -46,7 +49,6 @@ struct SlabSeaIceModel{GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} < ice_consolidation_thickness :: MIT # Numerics advection :: A - architecture end const SSIM = SlabSeaIceModel @@ -159,7 +161,8 @@ function SlabSeaIceModel(grid; architecture = architecture(grid) - return SlabSeaIceModel(grid, + return SlabSeaIceModel(architecture, + grid, clock, timestepper, ice_thickness, @@ -172,8 +175,7 @@ function SlabSeaIceModel(grid; internal_heat_flux_function, phase_transitions, ice_consolidation_thickness, - advection, - architecture) + advection) end function set!(model::SSIM; h=nothing, α=nothing) From a15fc560437dbec0b530e42656190abc752d71f6 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 31 Oct 2023 16:56:40 -0500 Subject: [PATCH 43/49] Still fixing grid architecture field --- src/SlabSeaIceModels/slab_sea_ice_model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SlabSeaIceModels/slab_sea_ice_model.jl b/src/SlabSeaIceModels/slab_sea_ice_model.jl index e9425ef9..2b45c6e6 100644 --- a/src/SlabSeaIceModels/slab_sea_ice_model.jl +++ b/src/SlabSeaIceModels/slab_sea_ice_model.jl @@ -159,9 +159,9 @@ function SlabSeaIceModel(grid; heat_boundary_conditions = (top = top_heat_boundary_condition, bottom = bottom_heat_boundary_condition) - architecture = architecture(grid) + arch = grid.architecture # architecture(grid) - return SlabSeaIceModel(architecture, + return SlabSeaIceModel(arch, grid, clock, timestepper, From e59a723320bc3b4113df9c2b3de6ebf67d6d1ea7 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 31 Oct 2023 17:14:15 -0500 Subject: [PATCH 44/49] Adding architecture field to ice_ocean_model struct --- validation/ice_ocean_model/ice_ocean_model.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/validation/ice_ocean_model/ice_ocean_model.jl b/validation/ice_ocean_model/ice_ocean_model.jl index 8dd9a61a..c5a8d660 100644 --- a/validation/ice_ocean_model/ice_ocean_model.jl +++ b/validation/ice_ocean_model/ice_ocean_model.jl @@ -18,7 +18,8 @@ import Oceananigans.Simulations: reset!, initialize!, iteration import Oceananigans.TimeSteppers: time_step!, update_state!, time import Oceananigans.Utils: prettytime -struct IceOceanModel{FT, C, G, I, O, S, PI, PC} <: AbstractModel{Nothing} +struct IceOceanModel{AR<:AbstractArchitecture, FT, C, G, I, O, S, PI, PC} <: AbstractModel{Nothing} + architecture :: AR clock :: C grid :: G # TODO: make it so simulation does not require this ice :: I @@ -77,7 +78,10 @@ function IceOceanModel(ice, ocean; clock = Clock{Float64}(0, 0, 1)) FT = eltype(ocean.model.grid) - return IceOceanModel(clock, + arch = ocean.model.grid.architecture + + return IceOceanModel(arch, + clock, ocean.model.grid, ice, previous_ice_thickness, From 2b1196ef8b948f2c845b40f627f174db09f0cbc3 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 31 Oct 2023 17:17:43 -0500 Subject: [PATCH 45/49] Add AbstractArchitecture --- validation/ice_ocean_model/ice_ocean_model.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/validation/ice_ocean_model/ice_ocean_model.jl b/validation/ice_ocean_model/ice_ocean_model.jl index c5a8d660..ec28fa11 100644 --- a/validation/ice_ocean_model/ice_ocean_model.jl +++ b/validation/ice_ocean_model/ice_ocean_model.jl @@ -12,6 +12,8 @@ using KernelAbstractions.Extras.LoopInfo: @unroll # Simulations interface import Oceananigans: fields, prognostic_fields import Oceananigans.Fields: set! +import Oceananigans.Architectures: AbstractArchitecture +import Oceananigans.Grids: architecture import Oceananigans.Models: timestepper, NaNChecker, default_nan_checker import Oceananigans.OutputWriters: default_included_properties import Oceananigans.Simulations: reset!, initialize!, iteration From 1c5790fe1ac478968396a14446562583bda3809a Mon Sep 17 00:00:00 2001 From: jlk9 Date: Tue, 31 Oct 2023 21:37:37 -0500 Subject: [PATCH 46/49] Set up model run with ice --- examples/antarctic_slope_current_model.jl | 369 ++++++++++++------ validation/ice_ocean_model/ice_ocean_model.jl | 5 +- 2 files changed, 243 insertions(+), 131 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 67a072d0..6829c868 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -2,17 +2,26 @@ using CairoMakie using Printf using Oceananigans -using Oceananigans.Fields, Oceananigans.AbstractOperations +using Oceananigans.Architectures: arch_array +using Oceananigans.Fields: ZeroField, ConstantField +using Oceananigans.AbstractOperations using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity # FIND submodule with CATKE using Oceananigans.Units: minute, minutes, hour, days, kilometers -using SeawaterPolynomials.TEOS10 + +using SeawaterPolynomials: TEOS10EquationOfState, haline_contraction + +using Statistics using ClimaSeaIce +using ClimaSeaIce: melting_temperature +using ClimaSeaIce.HeatBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium + +include("../validation/ice_ocean_model/ice_ocean_model.jl") -# This file sets up a model that resembles the Antarctic Slope Current (ASC) model in the +# This file sets up a ocean_model that resembles the Antarctic Slope Current (ASC) ocean_model in the # 2022 paper by Si, Stewart, and Eisenman -arch = GPU() +arch = CPU() g_Earth = 9.80665 @@ -24,10 +33,16 @@ Nx = 400 Ny = 450 Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +x = (-Lx/2, Lx/2) +y = (0, Ly) + +topology = (Periodic, Bounded, Bounded) +halo = (4, 4, 4) + sponge_width = 20kilometers # -# Setting up the grid and bathymetry: +# Setting up the ocean_grid and bathymetry: # refinement = 20.0 # controls spacing near surface (higher means finer spaced), 12.0 stretching = 3 # controls rate of stretching at bottom, 3 @@ -44,15 +59,23 @@ h(k) = (k - 1) / Nz # Generating function z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1) -underlying_grid = RectilinearGrid(arch, +ice_grid = RectilinearGrid(arch, + size = (Nx, Ny), + topology = (topology[1], topology[2], Flat), + x = x, + y = y, + halo = halo[1:2]) + +underlying_ocean_grid = RectilinearGrid(arch, size = (Nx, Ny, Nz), - topology = (Periodic, Bounded, Bounded), - x = (-Lx/2, Lx/2), - y = (0, Ly), + topology = topology, + x = x, + y = y, z = z_faces, - halo = (4, 4, 4)) + halo = halo) -@show underlying_grid +@show ice_grid +@show underlying_ocean_grid # We want the underwater slope to provide a depth of 500 m at y = 0 and the full 4 km by y =200. It follows # a hyperbolic tangent curve with its center at y ~= 150 at a depth of ~ -4000 + (4000 - 500)/2 = -2250 m @@ -68,14 +91,14 @@ step(y, y₀, Δy) = (1 + tanh((y - y₀) / Δy)) / 2 underwater_slope(x, y) = (-basin_depth - slope_depth + (slope_depth - basin_depth) * tanh((y - y₀) / Δy)) / 2 # TODO: add 50km troughs to the slope, dependent on x and y. Use appendix C to add detail here -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(underwater_slope)) +ocean_grid = ImmersedBoundaryGrid(underlying_ocean_grid, GridFittedBottom(underwater_slope)) -@show grid +@show ocean_grid # TODO: add underwater slope at y = 0 with troughs # -# Creating the Hydrostatic model: +# Creating the Hydrostatic ocean_model: # # Forcings (ignored for now): @@ -102,6 +125,14 @@ sponge_layers_T = (south_sponge_layer, north_sponge_layer_T) # Boundary Conditions: # +# Top boundary conditions for ice/ocean: +# - outgoing radiative fluxes emitted from surface +# - incoming shortwave radiation starting after 40 days + +ice_ocean_heat_flux = Field{Center, Center, Nothing}(ice_grid) +top_ocean_heat_flux = Qᵀ = Field{Center, Center, Nothing}(ice_grid) +top_salt_flux = Qˢ = Field{Center, Center, Nothing}(ice_grid) + """ boundary_ramp(y, Δy) @@ -132,21 +163,20 @@ cᴰ = 2.5e-3 # dimensionless drag coefficient u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ, parameters=(; cᴰ, ρₐ, ρₒ, Ly))) v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ, parameters=(; cᴰ, ρₐ, ρₒ, Ly))) +# Generate a zero-dimensional grid for a single column slab model boundary conditions: +boundary_conditions = (T = FieldBoundaryConditions(top=FluxBoundaryCondition(Qᵀ)), + S = FieldBoundaryConditions(top=FluxBoundaryCondition(Qˢ))) + # Buoyancy Equations of State - we want high order polynomials, so we'll use TEOS-10 eos = TEOS10EquationOfState() # can compare to linear EOS later (linear not recommended for polar regions) # Coriolis Effect, using basic f-plane with precribed reference Coriolis parameter coriolis = BetaPlane(f₀=1.3e-4, β=1e-11) #FPlane(latitude=-60) -# Diffusivities as part of closure -# TODO: make sure this works for biharmonic diffusivities as the horizontal, -horizontal_closure = HorizontalScalarDiffusivity(ν=0.1, κ=0.1) -vertical_closure = VerticalScalarDiffusivity(ν=3e-4, κ=1e-5) - # # Assuming no particles or biogeochemistry # -model = HydrostaticFreeSurfaceModel(; grid = grid, +ocean_model = HydrostaticFreeSurfaceModel(; grid = ocean_grid, momentum_advection = WENO(), tracer_advection = WENO(), buoyancy = SeawaterBuoyancy(equation_of_state=eos, gravitational_acceleration=g_Earth), @@ -158,7 +188,33 @@ model = HydrostaticFreeSurfaceModel(; grid = grid, tracers = (:T, :S, :e) ) +Nz = size(ocean_grid, 3) +So = ocean_model.tracers.S +ocean_surface_salinity = view(So, :, :, Nz) +bottom_bc = IceWaterThermalEquilibrium(ConstantField(30)) #ocean_surface_salinity) + +u, v, w = ocean_model.velocities +ocean_surface_velocities = (u = view(u, :, :, Nz), #interior(u, :, :, Nz), + v = view(v, :, :, Nz), #interior(v, :, :, Nz), + w = ZeroField()) + +ice_model = SlabSeaIceModel(ice_grid; + velocities = ocean_surface_velocities, + advection = nothing, #WENO(), + ice_consolidation_thickness = 0.05, + ice_salinity = 4, + internal_heat_flux = ConductiveFlux(conductivity=2), + #top_heat_flux = ConstantField(-100), # W m⁻² + top_heat_flux = ConstantField(0), # W m⁻² + top_heat_boundary_condition = PrescribedTemperature(0), + bottom_heat_boundary_condition = bottom_bc, + bottom_heat_flux = ice_ocean_heat_flux) + + + +# # INITIAL CONDITION FOR TEMPERATURE: a baroclinically unstable temperature distribution +# """ ramp(y, Δy) @@ -182,57 +238,192 @@ M² = 1e-5 # [s⁻²] meridional temperature gradient, was 1e-7 Tᵢ(x, y, z) = N² * z + ΔT * ramp(y, Δy) + ϵT * randn() -set!(model, T=Tᵢ) +# INITIAL CONDITION FOR ICE SALINITY: +function hᵢ(x, y) + if sqrt(x^2 + (y-225kilometers)^2) < 100kilometers + #return 1 + 0.1 * rand() + return 2 + else + return 0 + end +end + +set!(ocean_model, T=Tᵢ) +set!(ice_model, h=hᵢ) # TODO: impose a temperature restoring at the surface that corresponds to Si/Stewart to help get baroclinic eddies # Add initial condition to get turbulence running (should be doable on a laptop), then: # Buoyancy gradient imposed at North/South boundary, surface temp as freezing temperature imposed -@show model +@show ice_model +@show ocean_model # -# Now create a simulation and run the model +# Now create a ocean_simulation and run the ocean_model # # Full resolution is 100 sec -simulation = Simulation(model; Δt=100.0, stop_time=2days) +time_step = 100 +stop_time = 60days + +ice_simulation = Simulation(ice_model; Δt=time_step, stop_time=stop_time, verbose=false) +ocean_simulation = Simulation(ocean_model; Δt=time_step, stop_time=stop_time, verbose=false) + +filename = "asc_model_hi_res_60_days_ice" + +coupled_model = IceOceanModel(ice_simulation, ocean_simulation) +coupled_simulation = Simulation(coupled_model, Δt=time_step, stop_time=stop_time) + +S = ocean_model.tracers.S +# Initial condition +S₀ = 30 +T₀ = melting_temperature(ice_model.phase_transitions.liquidus, S₀) + 2.0 + +N²S = 1e-6 +β = haline_contraction(T₀, S₀, 0, eos) +g = ocean_model.buoyancy.model.gravitational_acceleration +by = - g * β * ∂y(S) + +function progress(sim) + h = sim.model.ice.model.ice_thickness + S = sim.model.ocean.model.tracers.S + T = sim.model.ocean.model.tracers.T + u = sim.model.ocean.model.velocities.u + msg1 = @sprintf("Iter: % 6d, time: % 12s", iteration(sim), prettytime(sim)) + msg2 = @sprintf(", max(h): %.2f", maximum(h)) + msg3 = @sprintf(", min(S): %.2f", minimum(S)) + msg4 = @sprintf(", extrema(T): (%.2f, %.2f)", minimum(T), maximum(T)) + msg5 = @sprintf(", max|∂y b|: %.2e", maximum(abs, by)) + msg6 = @sprintf(", max|u|: %.2e", maximum(abs, u)) + @info msg1 * msg2 * msg3 * msg4 * msg5 * msg6 + return nothing +end + +coupled_simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) + +H = ice_model.ice_thickness +T = ocean_model.tracers.T +S = ocean_model.tracers.S +u, v, w = ocean_model.velocities +η = ocean_model.free_surface.η + +Ht = [] +Tt = [] +Ft = [] +Qt = [] +St = [] +ut = [] +vt = [] +ηt = [] +ζt = [] +tt = [] + +ζ = Field(∂x(v) - ∂y(u)) + +function saveoutput(sim) + compute!(ζ) + Hn = Array(interior(H, :, :, 1)) + Fn = Array(interior(Qˢ, :, :, 1)) + Qn = Array(interior(Qᵀ, :, :, 1)) + Tn = Array(interior(T, :, :, Nz)) + Sn = Array(interior(S, :, :, Nz)) + un = Array(interior(u, :, :, Nz)) + vn = Array(interior(v, :, :, Nz)) + ηn = Array(interior(η, :, :, 1)) + ζn = Array(interior(ζ, :, :, Nz)) + push!(Ht, Hn) + push!(Ft, Fn) + push!(Qt, Qn) + push!(Tt, Tn) + push!(St, Sn) + push!(ut, un) + push!(vt, vn) + push!(ηt, ηn) + push!(ζt, ζn) + push!(tt, time(sim)) +end -filename = "asc_model_hi_res_2_days_wind_stress_Nsq_neg8" +coupled_simulation.callbacks[:output] = Callback(saveoutput, IterationInterval(10)) -# Here we'll try also running a zonal average of the simulation: -u, v, w = model.velocities -avgT = Average(model.tracers.T, dims=1) -avgU = Average(u, dims=1) -avgV = Average(v, dims=1) -avgW = Average(w, dims=1) +run!(coupled_simulation) -# Here we'll compute and record vorticity and speed: -ω = ∂x(v) - ∂y(u) -s = sqrt(u^2 + v^2 + w^2) +##### +##### Viz for this model with ice +##### -iter_interval = 240 +set_theme!(Theme(fontsize=24)) -#simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; T=avgT, u=avgU, v=avgV, w=avgW); -# filename = filename * "_zonal_average.jld2", -# schedule = IterationInterval(iter_interval), -# overwrite_existing = true) +x = xnodes(ocean_grid, Center()) +y = ynodes(ocean_grid, Center()) +fig = Figure(resolution=(2400, 700)) -simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers), +axh = Axis(fig[1, 1], xlabel="x (km)", ylabel="y (km)", title="Ice thickness") +axT = Axis(fig[1, 2], xlabel="x (km)", ylabel="y (km)", title="Ocean surface temperature") +axS = Axis(fig[1, 3], xlabel="x (km)", ylabel="y (km)", title="Ocean surface salinity") +axZ = Axis(fig[1, 4], xlabel="x (km)", ylabel="y (km)", title="Ocean vorticity") + +Nt = length(tt) +slider = Slider(fig[2, 1:4], range=1:Nt, startvalue=Nt) +n = slider.value + +title = @lift string("Melt-driven baroclinic instability after ", prettytime(tt[$n])) +Label(fig[0, 1:3], title) + +Hn = @lift Ht[$n] +Fn = @lift Ft[$n] +Tn = @lift Tt[$n] +Sn = @lift St[$n] +un = @lift ut[$n] +vn = @lift vt[$n] +ηn = @lift ηt[$n] +ζn = @lift ζt[$n] +Un = @lift mean(ut[$n], dims=1)[:] + +x = x ./ 1e3 +y = y ./ 1e3 + +Stop = view(S, :, :, Nz) +Smax = maximum(Stop) +Smin = minimum(Stop) + +compute!(ζ) +ζtop = view(ζ, :, :, Nz) +ζmax = maximum(abs, ζtop) +ζlim = 2e-4 #ζmax / 2 + +heatmap!(axh, x, y, Hn, colorrange=(0, 1), colormap=:grays) +heatmap!(axT, x, y, Tn, colormap=:heat) +heatmap!(axS, x, y, Sn, colorrange = (29, 30), colormap=:haline) +heatmap!(axZ, x, y, ζn, colorrange=(-ζlim, ζlim), colormap=:redblue) + +#heatmap!(axZ, x, y, Tn, colormap=:heat) +#heatmap!(axF, x, y, Fn) + +fig #display(fig) + +record(fig, "salty_baroclinic_ice_cube.mp4", 1:Nt, framerate=48) do nn + @info string(nn) + n[] = nn +end + + +#= +ocean_simulation.output_writers[:slices] = JLD2OutputWriter(ocean_model, merge(ocean_model.velocities, ocean_model.tracers), filename = filename * "_surface.jld2", - indices = (:, :, grid.Nz), + indices = (:, :, ocean_grid.Nz), schedule = IterationInterval(iter_interval), overwrite_existing = true) -simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ω, s), +ocean_simulation.output_writers[:fields] = JLD2OutputWriter(ocean_model, (; ω, s), filename = filename * "_fields.jld2", - indices = (:, :, grid.Nz), + indices = (:, :, ocean_grid.Nz), schedule = IterationInterval(iter_interval), overwrite_existing = true) -run!(simulation) -@info "Simulation completed in " * prettytime(simulation.run_wall_time) -@show simulation +run!(ocean_simulation) +@info "Simulation completed in " * prettytime(ocean_simulation.run_wall_time) +@show ocean_simulation # # Make a figure and plot it @@ -268,8 +459,8 @@ fig = Figure(resolution = (4000, 2000)) axis_kwargs = (xlabel="x (m)", ylabel="y (m)", - aspect = AxisAspect(grid.Lx/grid.Ly), - limits = ((-grid.Lx/2, grid.Lx/2), (0, grid.Ly))) + aspect = AxisAspect(ocean_grid.Lx/ocean_grid.Ly), + limits = ((-ocean_grid.Lx/2, ocean_grid.Lx/2), (0, ocean_grid.Ly))) ax_w = Axis(fig[2, 1]; title = "Vertical velocity", axis_kwargs...) ax_T = Axis(fig[2, 3]; title = "Temperature", axis_kwargs...) @@ -365,84 +556,4 @@ frames = intro:length(times) record(fig, filename * "_fields.mp4", frames, framerate=8) do i n[] = i end - -#= -# -# Now do all the above, but for zonal average time series data: -# - -average_time_series = (u = FieldTimeSeries(average_filepath, "u"), - v = FieldTimeSeries(average_filepath, "v"), - w = FieldTimeSeries(average_filepath, "w"), - T = FieldTimeSeries(average_filepath, "T")) - -#@show average_time_series.u -#@show average_time_series.w -#@show average_time_series.T -#@show average_time_series.v - -# Coordinate arrays -avg_xw, avg_yw, avg_zw = nodes(average_time_series.w) -avg_xT, avg_yT, avg_zT = nodes(average_time_series.T) - -times = average_time_series.w.times -intro = 1 - -n = Observable(intro) - -avg_wₙ = @lift interior(average_time_series.w[$n], 1, :, :) -avg_Tₙ = @lift interior(average_time_series.T[$n], 1, :, :) -#avg_Sₙ = @lift interior(average_time_series.S[$n], 1, :, :) -avg_uₙ = @lift interior(average_time_series.u[$n], 1, :, :) -avg_vₙ = @lift interior(average_time_series.v[$n], 1, :, :) - -fig = Figure(resolution = (4000, 2000)) - -axis_kwargs = (xlabel="y (m)", - ylabel="z (m)", - aspect = AxisAspect(2.),#AxisAspect(grid.Ly/grid.Lz), - limits = ((0, grid.Ly), (-grid.Lz, 0))) - -ax_w = Axis(fig[2, 1]; title = "Vertical velocity", axis_kwargs...) -ax_T = Axis(fig[2, 3]; title = "Temperature", axis_kwargs...) -#ax_S = Axis(fig[3, 1]; title = "Salinity", axis_kwargs...) -ax_v = Axis(fig[3, 1]; title = "Meridional velocity", axis_kwargs...) -ax_u = Axis(fig[3, 3]; title = "Zonal velocity", axis_kwargs...) - -title = @lift @sprintf("t = %s", prettytime(times[$n])) - -wlims = (-0.002, 0.002) -Tlims = (-0.02, 0.02) -Slims = (35, 35.005) -ulims = (-0.12, 0.12) -vlims = (-0.12, 0.12) - - -hm_w = heatmap!(ax_w, avg_yw, avg_zw, avg_wₙ; colormap = :balance) #, colorrange = wlims) -Colorbar(fig[2, 2], hm_w; label = "m s⁻¹") - -hm_T = heatmap!(ax_T, avg_yT, avg_zT, avg_Tₙ; colormap = :thermal) #, colorrange = Tlims) -Colorbar(fig[2, 4], hm_T; label = "ᵒC") - -#hm_S = heatmap!(ax_S, srf_xT, srf_yT, srf_Sₙ; colormap = :haline)#, colorrange = Slims) -#Colorbar(fig[3, 2], hm_S; label = "g / kg") - -hm_v = heatmap!(ax_v, avg_yw, avg_zw, avg_vₙ; colormap = :balance) #, colorrange = vlims) -Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") - -hm_u = heatmap!(ax_u, avg_yw, avg_zw, avg_uₙ; colormap = :balance) #, colorrange = ulims) -Colorbar(fig[3, 4], hm_u; label = "m s⁻¹") - -fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false) - -current_figure() # hide -fig - -frames = intro:length(times) - -@info "Making a motion picture of ocean wind mixing and convection..." - -record(fig, filename * "_average.mp4", frames, framerate=8) do i - n[] = i -end -=# \ No newline at end of file +=# diff --git a/validation/ice_ocean_model/ice_ocean_model.jl b/validation/ice_ocean_model/ice_ocean_model.jl index ec28fa11..fbb8d9c3 100644 --- a/validation/ice_ocean_model/ice_ocean_model.jl +++ b/validation/ice_ocean_model/ice_ocean_model.jl @@ -177,8 +177,9 @@ function time_step!(coupled_model::IceOceanModel, Δt; callbacks=nothing) time_step!(ice) # TODO: put this in update_state! - compute_ice_ocean_salinity_flux!(coupled_model) - ice_ocean_latent_heat!(coupled_model) + #compute_ice_ocean_salinity_flux!(coupled_model) + #ice_ocean_latent_heat!(coupled_model) + # AVOID FOR NOW: #compute_solar_insolation!(coupled_model) #compute_air_sea_flux!(coupled_model) From 10c77ff0a10942957082bcc851052d3807da1797 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 8 Nov 2023 12:16:34 -0600 Subject: [PATCH 47/49] Updated ice-ocean model, added architecture functions to SlabSeaIceModel and OceanIceModel --- examples/antarctic_slope_current_model.jl | 14 +++++++------- src/SlabSeaIceModels/slab_sea_ice_model.jl | 9 +++------ 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index 6829c868..f5d9a4e7 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -29,9 +29,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 400 -Ny = 450 -Nz = 70 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 64 +Ny = 64 +Nz = 8 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor x = (-Lx/2, Lx/2) y = (0, Ly) @@ -263,13 +263,13 @@ set!(ice_model, h=hᵢ) # # Full resolution is 100 sec -time_step = 100 -stop_time = 60days +time_step = 20minutes +stop_time = 2days ice_simulation = Simulation(ice_model; Δt=time_step, stop_time=stop_time, verbose=false) ocean_simulation = Simulation(ocean_model; Δt=time_step, stop_time=stop_time, verbose=false) -filename = "asc_model_hi_res_60_days_ice" +filename = "asc_model_lo_res_2_days_ice" coupled_model = IceOceanModel(ice_simulation, ocean_simulation) coupled_simulation = Simulation(coupled_model, Δt=time_step, stop_time=stop_time) @@ -402,7 +402,7 @@ heatmap!(axZ, x, y, ζn, colorrange=(-ζlim, ζlim), colormap=:redblue) fig #display(fig) -record(fig, "salty_baroclinic_ice_cube.mp4", 1:Nt, framerate=48) do nn +record(fig, filename * ".mp4", 1:Nt, framerate=48) do nn @info string(nn) n[] = nn end diff --git a/src/SlabSeaIceModels/slab_sea_ice_model.jl b/src/SlabSeaIceModels/slab_sea_ice_model.jl index 2b45c6e6..5350aaee 100644 --- a/src/SlabSeaIceModels/slab_sea_ice_model.jl +++ b/src/SlabSeaIceModels/slab_sea_ice_model.jl @@ -28,8 +28,7 @@ import Oceananigans.Utils: prettytime # import Oceananigans.Fields: field # field(loc, a::Number, grid) = ConstantField(a) -struct SlabSeaIceModel{AR<:AbstractArchitecture, GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} <: AbstractModel{TS} - architecture :: AR +struct SlabSeaIceModel{GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} <: AbstractModel{TS} grid :: GR clock :: CL timestepper :: TS @@ -56,6 +55,7 @@ const SSIM = SlabSeaIceModel Base.summary(model::SSIM) = "SlabSeaIceModel" prettytime(model::SSIM) = prettytime(model.clock.time) iteration(model::SSIM) = model.clock.iteration +architecture(model::SSIM) = architecture(model.grid) function Base.show(io::IO, model::SSIM) grid = model.grid @@ -159,10 +159,7 @@ function SlabSeaIceModel(grid; heat_boundary_conditions = (top = top_heat_boundary_condition, bottom = bottom_heat_boundary_condition) - arch = grid.architecture # architecture(grid) - - return SlabSeaIceModel(arch, - grid, + return SlabSeaIceModel(grid, clock, timestepper, ice_thickness, From 4d6be262455b130f6625c71cdf7d901301a63d88 Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 8 Nov 2023 12:18:00 -0600 Subject: [PATCH 48/49] Added architecture functions to OceanIceModel --- validation/ice_ocean_model/ice_ocean_model.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/validation/ice_ocean_model/ice_ocean_model.jl b/validation/ice_ocean_model/ice_ocean_model.jl index fbb8d9c3..6ff11f07 100644 --- a/validation/ice_ocean_model/ice_ocean_model.jl +++ b/validation/ice_ocean_model/ice_ocean_model.jl @@ -20,8 +20,7 @@ import Oceananigans.Simulations: reset!, initialize!, iteration import Oceananigans.TimeSteppers: time_step!, update_state!, time import Oceananigans.Utils: prettytime -struct IceOceanModel{AR<:AbstractArchitecture, FT, C, G, I, O, S, PI, PC} <: AbstractModel{Nothing} - architecture :: AR +struct IceOceanModel{FT, C, G, I, O, S, PI, PC} <: AbstractModel{Nothing} clock :: C grid :: G # TODO: make it so simulation does not require this ice :: I @@ -48,6 +47,7 @@ default_included_properties(::IOM) = tuple() update_state!(::IOM) = nothing prognostic_fields(cm::IOM) = nothing fields(::IOM) = NamedTuple() +architecture(cm::IOM) = architecture(cm.grid) function IceOceanModel(ice, ocean; clock = Clock{Float64}(0, 0, 1)) @@ -80,10 +80,7 @@ function IceOceanModel(ice, ocean; clock = Clock{Float64}(0, 0, 1)) FT = eltype(ocean.model.grid) - arch = ocean.model.grid.architecture - - return IceOceanModel(arch, - clock, + return IceOceanModel(clock, ocean.model.grid, ice, previous_ice_thickness, From 64aa547109c25b11c5996dd2183f1aefb28b8dca Mon Sep 17 00:00:00 2001 From: jlk9 Date: Wed, 8 Nov 2023 13:45:29 -0600 Subject: [PATCH 49/49] Fixed vorticity plot, added heat flux computations back in --- examples/antarctic_slope_current_model.jl | 20 ++++++++++--------- validation/ice_ocean_model/ice_ocean_model.jl | 2 +- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/antarctic_slope_current_model.jl b/examples/antarctic_slope_current_model.jl index f5d9a4e7..51412ded 100644 --- a/examples/antarctic_slope_current_model.jl +++ b/examples/antarctic_slope_current_model.jl @@ -29,9 +29,9 @@ Lx = 400kilometers Ly = 450kilometers Lz = 4000 -Nx = 64 -Ny = 64 -Nz = 8 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor +Nx = 400 #64 +Ny = 450 #64 +Nz = 70 #8 # TODO: modify spacing if needed, 10 m at surface, 100m at seafloor x = (-Lx/2, Lx/2) y = (0, Ly) @@ -263,13 +263,13 @@ set!(ice_model, h=hᵢ) # # Full resolution is 100 sec -time_step = 20minutes -stop_time = 2days +time_step = 100 +stop_time = 60days ice_simulation = Simulation(ice_model; Δt=time_step, stop_time=stop_time, verbose=false) ocean_simulation = Simulation(ocean_model; Δt=time_step, stop_time=stop_time, verbose=false) -filename = "asc_model_lo_res_2_days_ice" +filename = "asc_model_hi_res_60_days_ice" coupled_model = IceOceanModel(ice_simulation, ocean_simulation) coupled_simulation = Simulation(coupled_model, Δt=time_step, stop_time=stop_time) @@ -330,7 +330,7 @@ function saveoutput(sim) un = Array(interior(u, :, :, Nz)) vn = Array(interior(v, :, :, Nz)) ηn = Array(interior(η, :, :, 1)) - ζn = Array(interior(ζ, :, :, Nz)) + ζn = Array(interior(ζ, :, :, 1)) push!(Ht, Hn) push!(Ft, Fn) push!(Qt, Qn) @@ -390,7 +390,9 @@ Smin = minimum(Stop) compute!(ζ) ζtop = view(ζ, :, :, Nz) ζmax = maximum(abs, ζtop) -ζlim = 2e-4 #ζmax / 2 +ζlim = 5e-5 #ζmax / 2 + +@show ζmax heatmap!(axh, x, y, Hn, colorrange=(0, 1), colormap=:grays) heatmap!(axT, x, y, Tn, colormap=:heat) @@ -403,7 +405,7 @@ heatmap!(axZ, x, y, ζn, colorrange=(-ζlim, ζlim), colormap=:redblue) fig #display(fig) record(fig, filename * ".mp4", 1:Nt, framerate=48) do nn - @info string(nn) + #@info string(nn) n[] = nn end diff --git a/validation/ice_ocean_model/ice_ocean_model.jl b/validation/ice_ocean_model/ice_ocean_model.jl index 6ff11f07..d73c74af 100644 --- a/validation/ice_ocean_model/ice_ocean_model.jl +++ b/validation/ice_ocean_model/ice_ocean_model.jl @@ -175,7 +175,7 @@ function time_step!(coupled_model::IceOceanModel, Δt; callbacks=nothing) # TODO: put this in update_state! #compute_ice_ocean_salinity_flux!(coupled_model) - #ice_ocean_latent_heat!(coupled_model) + ice_ocean_latent_heat!(coupled_model) # AVOID FOR NOW: #compute_solar_insolation!(coupled_model) #compute_air_sea_flux!(coupled_model)