diff --git a/docs/make.jl b/docs/make.jl index 3fbeb5610..432c5e59e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -17,7 +17,7 @@ to_be_literated = [ "generate_bathymetry.jl", "generate_surface_fluxes.jl", "single_column_os_papa_simulation.jl", - # "mediterranean_simulation_with_ecco_restoring.jl", + "mediterranean_simulation_with_ecco_restoring.jl", "near_global_ocean_simulation.jl" ] @@ -44,8 +44,8 @@ pages = [ "Inspect ECCO2 data" => "literated/ecco_inspect_temperature_salinity.md", "Generate bathymetry" => "literated/generate_bathymetry.md", "Surface fluxes" => "literated/generate_surface_fluxes.md", - "Single-column simulation" => "literated/single_column_os_papa_simulation.md", - # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md", + "Single column simulation" => "literated/single_column_os_papa_simulation.md", + "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md", "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md", ], diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 4cfd35ddd..be1762188 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -18,7 +18,7 @@ using CairoMakie # # Computing fluxes on the ECCO2 grid # -# We start by building the ECCO2 grid, using `ECCO_bottom_height` to identify the bottom height. +# We start by building the ECCO2 grid, using `ECCO_immersed_grid` function. grid = ECCO_immersed_grid() diff --git a/examples/mediterranean_simulation_with_ecco_restoring.jl b/examples/mediterranean_simulation_with_ecco_restoring.jl index 9d640820f..4a1253101 100644 --- a/examples/mediterranean_simulation_with_ecco_restoring.jl +++ b/examples/mediterranean_simulation_with_ecco_restoring.jl @@ -42,8 +42,8 @@ z_faces = stretched_vertical_faces(depth = 5000, stretching = PowerLawStretching(1.070), surface_layer_height = 50) -Nx = 15 * Int(λ₂ - λ₁) # 1/15 th of a degree resolution -Ny = 15 * Int(φ₂ - φ₁) # 1/15 th of a degree resolution +Nx = 15 * Int(λ₂ - λ₁) # 1/15th of a degree resolution +Ny = 15 * Int(φ₂ - φ₁) # 1/15th of a degree resolution Nz = length(z_faces) - 1 grid = LatitudeLongitudeGrid(GPU(); @@ -60,33 +60,34 @@ grid = LatitudeLongitudeGrid(GPU(); # are adjusted to refine the bathymetry representation. bottom_height = regrid_bathymetry(grid, - height_above_water = 1, minimum_depth = 10, - interpolation_passes = 25, - connected_regions_allowed = 1) + interpolation_passes = 10, + major_basins = 1) grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) -# ## Downloading ECCO data +# ## ECCO Restoring # -# The model is initialized with temperature and salinity fields from the ECCO dataset, -# using the function `ECCO_restoring_forcing` to apply restoring forcings for these tracers. +# The model is restored to at the surface to the temperature and salinity fields from the ECCO dataset. +# We build the restoring using the `ECCORestoring` functionality. # This allows us to nudge the model towards realistic temperature and salinity profiles. +# `ECCORestoring` accepts a `mask` keyword argument to restrict the restoring region. -dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 12, 1) +@inline surface_mask(x, y, z, t) = z > - 50 -temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly()) -salinity = ECCOMetadata(:salinity, dates, ECCO4Monthly()) +start = DateTimeProlepticGregorian(1993, 1, 1) +stop = DateTimeProlepticGregorian(1993, 12, 1) +dates = range(start, stop; step=Month(1)) -FT = ECCO_restoring_forcing(temperature; architecture = GPU(), timescale = 2days) -FS = ECCO_restoring_forcing(salinity; architecture = GPU(), timescale = 2days) +FT = ECCORestoring(:temperature, grid; dates, mask=surface_mask, rate=1/5days) +FS = ECCORestoring(:salinity, grid; dates, mask=surface_mask, rate=1/5days) # Constructing the Simulation # # We construct an ocean simulation that evolves two tracers, temperature (:T), salinity (:S) # and we pass the previously defined forcing that nudge these tracers -ocean = ocean_simulation(grid; forcing = (T = FT, S = FS)) +ocean = ocean_simulation(grid; forcing = (T=FT, S=FS)) # Initializing the model # @@ -94,13 +95,17 @@ ocean = ocean_simulation(grid; forcing = (T = FT, S = FS)) # In this case, our ECCO dataset has access to a temperature and a salinity # field, so we initialize temperature T and salinity S from ECCO. -set!(ocean.model, T = temperature[1], S = salinity[1]) +set!(ocean.model, T=ECCOMetadata(:temperature; dates=dates[1]), + S=ECCOMetadata(:salinity; dates=dates[1])) fig = Figure() ax = Axis(fig[1, 1]) -heatmap!(ax, interior(model.tracers.T, :, :, Nz), colorrange = (10, 20), colormap = :thermal) +heatmap!(ax, view(ocean.model.tracers.T, :, :, Nz), colorrange = (10, 20), colormap = :thermal) ax = Axis(fig[1, 2]) -heatmap!(ax, interior(model.tracers.S, :, :, Nz), colorrange = (35, 40), colormap = :haline) +heatmap!(ax, view(ocean.model.tracers.S, :, :, Nz), colorrange = (35, 40), colormap = :haline) + +save("initial_conditions.png", fig) +# ![](initial_conditions.png) function progress(sim) u, v, w = sim.model.velocities @@ -129,17 +134,13 @@ ocean.stop_iteration = 1000 run!(ocean) # ## Run the real simulation -# -# Now that the solution has adjusted to the bathymetry we can ramp up the time -# step size. We use a `TimeStepWizard` to automatically adapt to a CFL of 0.2. - -wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 10minutes, max_change = 1.1) -ocean.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) +# Let's reset the maximum number of iterations and we can now increase the time step size -# Let's reset the maximum number of iterations +ocean.Δt = 3minutes ocean.stop_iteration = Inf -ocean.stop_time = 200days +ocean.stop_time = 100days +model = ocean.model ocean.output_writers[:surface_fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers); indices = (:, :, Nz), @@ -157,52 +158,29 @@ run!(ocean) # (3) Temperature (T) # (4) Salinity (S) -u_series = FieldTimeSeries("med_surface_field.jld2", "u") -v_series = FieldTimeSeries("med_surface_field.jld2", "v") -T_series = FieldTimeSeries("med_surface_field.jld2", "T") -S_series = FieldTimeSeries("med_surface_field.jld2", "S") -c_series = FieldTimeSeries("med_surface_field.jld2", "c") +u_series = FieldTimeSeries("med_surface_field.jld2", "u"; backend=InMemory(10)) +v_series = FieldTimeSeries("med_surface_field.jld2", "v"; backend=InMemory(10)) +T_series = FieldTimeSeries("med_surface_field.jld2", "T"; backend=InMemory(10)) +S_series = FieldTimeSeries("med_surface_field.jld2", "S"; backend=InMemory(10)) iter = Observable(1) -u = @lift begin - f = interior(u_series[$iter], :, :, 1) - f[f .== 0] .= NaN - f -end -v = @lift begin - f = interior(v_series[$iter], :, :, 1) - f[f .== 0] .= NaN - f -end -T = @lift begin - f = interior(T_series[$iter], :, :, 1) - f[f .== 0] .= NaN - f -end -S = @lift begin - f = interior(S_series[$iter], :, :, 1) - f[f .== 0] .= NaN - f -end -c = @lift begin - f = interior(c_series[$iter], :, :, 1) - f[f .== 0] .= NaN - f -end +u = @lift(u_series[$iter]) +v = @lift(v_series[$iter]) +T = @lift(T_series[$iter]) +S = @lift(S_series[$iter]) fig = Figure() ax = Axis(fig[1, 1], title = "surface zonal velocity ms⁻¹") -heatmap!(u) +heatmap!(ax, u) ax = Axis(fig[1, 2], title = "surface meridional velocity ms⁻¹") -heatmap!(v) +heatmap!(ax, v) ax = Axis(fig[2, 1], title = "surface temperature ᵒC") -heatmap!(T) +heatmap!(ax, T) ax = Axis(fig[2, 2], title = "surface salinity psu") -heatmap!(S) -ax = Axis(fig[2, 3], title = "passive tracer -") -heatmap!(c) +heatmap!(ax, S) CairoMakie.record(fig, "mediterranean_video.mp4", 1:length(u_series.times); framerate = 5) do i @info "recording iteration $i" iter[] = i -end \ No newline at end of file +end +# ![](mediterranean_video.mp4)