From 3d4f21c3a01e7fe0cba7ac13024e7b943b23c97d Mon Sep 17 00:00:00 2001 From: nanophyto Date: Thu, 28 May 2026 10:02:58 +0100 Subject: [PATCH 1/2] delete old example --- 2026-05/examples/09_diffusivity.jl | 172 ----------------------------- 1 file changed, 172 deletions(-) delete mode 100644 2026-05/examples/09_diffusivity.jl diff --git a/2026-05/examples/09_diffusivity.jl b/2026-05/examples/09_diffusivity.jl deleted file mode 100644 index 3dd5ef4..0000000 --- a/2026-05/examples/09_diffusivity.jl +++ /dev/null @@ -1,172 +0,0 @@ -# # [Exercise 09: Diffusivity] (@id diffusivity_exercise) - -# This exercise introduces vertical diffusivity in a simple two-layer water-column model. - -# ## Loading dependencies -# The example uses Agate.jl, Oceananigans.jl, and OceanBioME.jl for the ocean simulations. -# CairoMakie.jl is used for plotting. - -using Agate -using Agate.Library.Light -using OceanBioME -using OceanBioME: Biogeochemistry -using Oceananigans -using Oceananigans.Units -using CairoMakie -workshop_script = let dir = @__DIR__ - while !isfile(joinpath(dir, "src", "AgateWorkshop.jl")) - parent = dirname(dir) - parent == dir && error("Could not find src/AgateWorkshop.jl") - dir = parent - end - joinpath(dir, "src", "AgateWorkshop.jl") -end -include(workshop_script) - -stop_time = 3*365day # simulate for 3 years -nothing #hide - -# ## Forcings - -# Second, we define the model physical forcings. Diffusivity is split across a 100 m interface on a two-level vertical grid, and PAR is held at its maximum surface value with a fixed attenuation coefficient. -#diffusivity -@inline function diffusivity(x, y, z, t) - κ_max = 1e-5 - layer_interface = -100meters - - if z >= layer_interface - return κ_max - else - return 0.0 - end -end - -#irradiance -function irradiance(x, y, z, t) - PAR_surface_max = 80 - layer_interface = -100meters - - return ifelse(z >= layer_interface, PAR_surface_max, 0.0) -end - -#plots -t_range = 0.0:days:(365.0 * days) # Time range from 0 to 365 days -z_range = [-150.0, -50.0] # Two 100 m layer centers -x, y, z = 0.0, 0.0, 0.0 -κₜ_values = [diffusivity(x, y, z, t) for t in t_range, z in z_range] -PAR_values = [irradiance(x, y, z, t) for t in t_range, z in z_range] - -fig_forcing = Figure(; size=(800, 600), fontsize=14) -ax1 = Axis(fig_forcing[1, 1]; xlabel="Time (days)", ylabel="Depth (m)", title="irradiance") -hm1 = CairoMakie.heatmap!(ax1, t_range ./ days, z_range, PAR_values; colormap=:viridis) -Colorbar(fig_forcing[1, 2], hm1) - -ax2 = Axis(fig_forcing[2, 1]; xlabel="Time (days)", ylabel="Depth (m)", title="diffusivity") -hm2 = CairoMakie.heatmap!(ax2, t_range ./ days, z_range, κₜ_values; colormap=:viridis) -Colorbar(fig_forcing[2, 2], hm2) - -display(fig_forcing) -fig_forcing - -# ## Physical model - -grid = RectilinearGrid(; size=(1, 1, 2), extent=(20meters, 20meters, 200meters)) -nothing #hide - -# ## Ecosystem model - -# First, we construct our ecosystem model. -# Here, we use a default 2 phytoplankton, 2 zooplankton `Agate.jl-NiPiZD` ecosystem model. -# Detritus sinks downward at 2 m/day; the closed bottom keeps sunk detritus in the lower box. - -bgc = Agate.Models.NiPiZD.construct(; -) -nothing #hide - -bgc_model = Biogeochemistry( - bgc; light_attenuation=FunctionFieldPAR(; grid, PAR_f=irradiance) -) -nothing #hide - -full_model = NonhydrostaticModel(; - grid, - clock=Clock(; time=0.0), - timestepper=:QuasiAdamsBashforth2, - closure=ScalarDiffusivity( - VerticallyImplicitTimeDiscretization(); ν=diffusivity, κ=diffusivity - ), - biogeochemistry=bgc_model, -) -nothing #hide - -# ## Initial conditions - -set!(full_model; default_initial_conditions(bgc; detritus = 0.0, total_plankton_biomass = 0.12)...) # mmol N / m³ - -# ## Simulation -filename = "N2P2ZD_column.jld2" - -simulation = Simulation(full_model; Δt=1hour, stop_time=stop_time) - -simulation.output_writers[:profiles] = JLD2Writer( - full_model, - full_model.tracers; - filename=filename, - schedule=TimeInterval(1day), - overwrite_existing=true, -) - -run!(simulation) -nothing #hide - -# ## Plotting - -#Load time series data -timeseries = NamedTuple{keys(full_model.tracers)}( - FieldTimeSeries(filename, "$field") for field in keys(full_model.tracers) -) - -timeseries_keys = keys(timeseries) -nothing #hide - -#Filter keys for P, Z, N, and D fields -P_keys = filter(k -> startswith(string(k), "P"), timeseries_keys) -Z_keys = filter(k -> startswith(string(k), "Z"), timeseries_keys) -N_key = :N -D_key = :D - -#Combine all keys into a single list for iteration -all_keys = [P_keys..., Z_keys..., N_key, D_key] - -#Create figure with appropriate size -fig = Figure(; size=(800, 1200), fontsize=16) - -#Plot all fields -for (i, key) in enumerate(all_keys) - x_nodes, y_nodes, z_nodes = nodes(timeseries[key]) - z_vals = collect(z_nodes) - times = collect(timeseries[key].times / days) - - ax = Axis( - fig[i, 1]; - title="$(key) concentration (mmol N / m³)", - xlabel="Time (days)", - ylabel="z (m)", - limits=((0, 365*3), (-200, 0)), - ) - hm = heatmap!( - ax, - times, - z_vals, - Float32.(interior(timeseries[key],1,1,:,:)'); - colormap=:viridis, - rasterize=true, - ) # Rasterize for smaller output - Colorbar(fig[i, 2], hm) -end - -#Save figure -save("N2P2ZD_column.png", fig) - -display(fig) -fig # Display the figure From cab81c969ba9d615d2a44e92b3424f3275baf3d4 Mon Sep 17 00:00:00 2001 From: nanophyto Date: Thu, 28 May 2026 10:05:31 +0100 Subject: [PATCH 2/2] delete broken example and re-org --- 2026-05/docs/make.jl | 14 +- 2026-05/examples/09_irradiance_box.jl | 136 ----------- ...ance_column.jl => 09_irradiance_column.jl} | 0 ...tified.jl => 10_diffusivity_stratified.jl} | 0 2026-05/examples/10_irradiance.jl | 211 ------------------ ...seasonal.jl => 11_diffusivity_seasonal.jl} | 0 6 files changed, 6 insertions(+), 355 deletions(-) delete mode 100644 2026-05/examples/09_irradiance_box.jl rename 2026-05/examples/{10_irradiance_column.jl => 09_irradiance_column.jl} (100%) rename 2026-05/examples/{11_diffusivity_stratified.jl => 10_diffusivity_stratified.jl} (100%) delete mode 100644 2026-05/examples/10_irradiance.jl rename 2026-05/examples/{12_diffusivity_seasonal.jl => 11_diffusivity_seasonal.jl} (100%) diff --git a/2026-05/docs/make.jl b/2026-05/docs/make.jl index 2968567..2194d36 100644 --- a/2026-05/docs/make.jl +++ b/2026-05/docs/make.jl @@ -24,10 +24,9 @@ example_files = [ "06_palatability.jl", "07_assimilation_efficiency.jl", "08_closure_terms.jl", - "09_irradiance_box.jl", - "10_irradiance_column.jl", - "11_diffusivity_stratified.jl", - "12_diffusivity_seasonal.jl", + "09_irradiance_column.jl", + "10_diffusivity_stratified.jl", + "11_diffusivity_seasonal.jl", ] function strip_jld2_warnings(content) @@ -95,10 +94,9 @@ makedocs( "06 Palatability" => "generated/06_palatability.md", "07 Assimilation Efficiency" => "generated/07_assimilation_efficiency.md", "08 Closure Terms" => "generated/08_closure_terms.md", - "09 Irradiance Box" => "generated/09_irradiance_box.md", - "10 Irradiance Column" => "generated/10_irradiance_column.md", - "11 Diffusivity Stratified" => "generated/11_diffusivity_stratified.md", - "12 Diffusivity Seasonal" => "generated/12_diffusivity_seasonal.md", + "09 Irradiance Column" => "generated/09_irradiance_column.md", + "10 Diffusivity Stratified" => "generated/10_diffusivity_stratified.md", + "11 Diffusivity Seasonal" => "generated/11_diffusivity_seasonal.md", ], ], ) diff --git a/2026-05/examples/09_irradiance_box.jl b/2026-05/examples/09_irradiance_box.jl deleted file mode 100644 index c0dde7d..0000000 --- a/2026-05/examples/09_irradiance_box.jl +++ /dev/null @@ -1,136 +0,0 @@ -# # [Exercise 09: Irradiance box model] (@id irradiance_box_exercise) -# -# This exercise introduces seasonal irradiance forcing in a zero-dimensional box model. - -# ## Loading dependencies - -using Agate -using Agate.Library.Light -using OceanBioME: BoxModelGrid -using Oceananigans.Units -using CairoMakie -workshop_script = let dir = @__DIR__ - while !isfile(joinpath(dir, "src", "AgateWorkshop.jl")) - parent = dirname(dir) - parent == dir && error("Could not find src/AgateWorkshop.jl") - dir = parent - end - joinpath(dir, "src", "AgateWorkshop.jl") -end -include(workshop_script) - -const year = years = 365day - -mkpath("outputs") -mkpath("figures") - -nothing #hide - -# ## Ecosystem model - -bgc = Agate.Models.NiPiZD.construct() -initial_conditions = default_initial_conditions(bgc) - -nothing #hide - -# ## Irradiance forcing - -@inline function seasonal_surface_PAR(t) - return 60 * - (1 - cos((t + 15days) * 2π / year)) * - (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2 -end - -@inline seasonal_PAR(t) = seasonal_surface_PAR(t) - -@inline nonseasonal_surface_PAR(t) = 100 * max(0, cos(t * π / 12hours)) -@inline nonseasonal_PAR(t) = nonseasonal_surface_PAR(t) - -t_range = 0.0:hours:(365.0days) -seasonal_PAR_values = [seasonal_PAR(t) for t in t_range] -nonseasonal_PAR_values = [nonseasonal_PAR(t) for t in t_range] - -fig_forcing = Figure(; size=(800, 350), fontsize=14) -ax = Axis(fig_forcing[1, 1]; xlabel="Time (days)", ylabel="PAR", title="Seasonal and non-seasonal irradiance") -lines!(ax, t_range ./ days, seasonal_PAR_values; linewidth=3, label="seasonal PAR") -lines!(ax, t_range ./ days, nonseasonal_PAR_values; linewidth=3, linestyle=:dash, label="non-seasonal PAR") -axislegend(ax; position=:rt) -save(joinpath("figures", "09_irradiance_box_forcing.png"), fig_forcing; px_per_unit=1) - -display(fig_forcing) -fig_forcing - -# ## Box-model simulations -# -# The diagnostic comparison uses two otherwise identical box-model runs: a -# non-seasonal reference and the seasonal-irradiance case. The reference -# uses the same default `FunctionFieldPAR(; grid=BoxModelGrid())` light forcing -# and default initial conditions as the workshop box-model helpers. - -nonseasonal_light = FunctionFieldPAR(; grid=BoxModelGrid()) -seasonal_light = FunctionFieldPAR(; grid=BoxModelGrid(), PAR_f=seasonal_PAR) - -nonseasonal_run = run_box_model( - bgc; - filename=joinpath("outputs", "09_irradiance_box_nonseasonal.jld2"), - initial_conditions, - stop_time=1year, - light_attenuation=nonseasonal_light, -) - -seasonal_run = run_box_model( - bgc; - filename=joinpath("outputs", "09_irradiance_box.jld2"), - initial_conditions, - stop_time=1year, - light_attenuation=seasonal_light, -) - -nonseasonal_timeseries = read_box_tracer_timeseries(nonseasonal_run.filename, nonseasonal_run.tracer_syms) -seasonal_timeseries = read_box_tracer_timeseries(seasonal_run.filename, seasonal_run.tracer_syms) - -nothing #hide - -# ## Tracer-concentration comparison -# -# `plot_box_timeseries` plots every shared tracer and overlays the two forcing -# cases in each panel. - -tracer_comparison_path = joinpath("figures", "09_irradiance_box_tracer_comparison.png") -fig_tracers = plot_box_timeseries( - [nonseasonal_timeseries, seasonal_timeseries]; - labels=["non-seasonal PAR", "seasonal PAR"], - figure_path=tracer_comparison_path, -) -save(joinpath("figures", "09_irradiance_box.png"), fig_tracers; px_per_unit=1) -display(fig_tracers) -fig_tracers - -# ## Nitrogen contributions -# -# The stacked-area diagnostic shows how living and non-living nitrogen pools vary -# through the seasonal-irradiance run. - -contributions_path = joinpath("figures", "09_irradiance_box_seasonal_contributions.png") -fig_contributions = plot_contributions( - seasonal_timeseries.times, - seasonal_timeseries.data; - figure_path=contributions_path, -) -display(fig_contributions) -fig_contributions - -# ## Community-weighted mean size comparison -# -# `plot_cwm_size` pairs each time series with the biogeochemistry object that -# defines the plankton diameters. Both runs use the same ecosystem model here. - -cwm_comparison_path = joinpath("figures", "09_irradiance_box_cwm_size_comparison.png") -fig_cwm = plot_cwm_size( - [nonseasonal_timeseries, seasonal_timeseries], - [bgc, bgc]; - labels=["non-seasonal PAR", "seasonal PAR"], - figure_path=cwm_comparison_path, -) -display(fig_cwm) -fig_cwm diff --git a/2026-05/examples/10_irradiance_column.jl b/2026-05/examples/09_irradiance_column.jl similarity index 100% rename from 2026-05/examples/10_irradiance_column.jl rename to 2026-05/examples/09_irradiance_column.jl diff --git a/2026-05/examples/11_diffusivity_stratified.jl b/2026-05/examples/10_diffusivity_stratified.jl similarity index 100% rename from 2026-05/examples/11_diffusivity_stratified.jl rename to 2026-05/examples/10_diffusivity_stratified.jl diff --git a/2026-05/examples/10_irradiance.jl b/2026-05/examples/10_irradiance.jl deleted file mode 100644 index 5bd17d9..0000000 --- a/2026-05/examples/10_irradiance.jl +++ /dev/null @@ -1,211 +0,0 @@ -# # [Exercise 10: Irradiance] (@id irradiance_exercise) - -# !!! info -# This example uses [Oceananigans.jl](https://clima.github.io/OceananigansDocumentation/stable/) and [OceanBioME.jl](https://oceanbiome.github.io/OceanBioME.jl/stable/). -# We recommend familiarizing yourself with their user interface if you intend to make changes to the physical model setup. - -# This exercise focuses on irradiance forcing in a simple 1D water-column model. -# The physical model setup is based on an example provided in the OceanBioME.jl documentation and represents an idealized 200m deep North Atlantic time series. - -# ## Loading dependencies -# The example uses Agate.jl, Oceananigans.jl, and OceanBioME.jl for the ocean simulations. -# CairoMakie is used for plotting. - -using Agate -using Agate.Introspection: tracer_groups -using Agate.Library.Light -using OceanBioME -using OceanBioME: Biogeochemistry -using Oceananigans -using Oceananigans.Units -using CairoMakie -workshop_script = let dir = @__DIR__ - while !isfile(joinpath(dir, "src", "AgateWorkshop.jl")) - parent = dirname(dir) - parent == dir && error("Could not find src/AgateWorkshop.jl") - dir = parent - end - joinpath(dir, "src", "AgateWorkshop.jl") -end -include(workshop_script) - -year = years = 365day -nothing #hide - -# ## Ecosystem model - -# First, we construct our ecosystem model. -# Here, we use a default 2 phytoplankton, 2 zooplankton `Agate.jl-NiPiZD` ecosystem model. - -bgc = Agate.Models.NiPiZD.construct() -groups = tracer_groups(bgc) -nothing #hide - -# ## Forcings - -# Second, we define the model physical forcings. Diffusivity is held high throughout the water column, while PAR is held at a fixed surface value with depth-dependent attenuation. - -#diffusivity -@inline diffusivity_profile(x, y, z, t) = 1e-4 - -#irradiance -@inline function constant_PAR(x, y, z, t) - PAR⁰ = 80 - attenuation = 0.04 - return PAR⁰ * exp(attenuation * z) -end - -#plots -t_range = 0.0:days:(365.0 * days) # Time range from 0 to 365 days -z_range = -200.0:10.0:0.0 # Depth range from -200m to 0m -x, y, z = 0.0, 0.0, 0.0 -κₜ_values = [diffusivity_profile(x, y, z, t) for t in t_range, z in z_range] -PAR_values = [constant_PAR(x, y, z, t) for t in t_range, z in z_range] - -fig_forcing = Figure(; size=(800, 600), fontsize=14) -ax1 = Axis(fig_forcing[1, 1]; xlabel="Time (days)", ylabel="Depth (m)", title="irradiance") -hm1 = CairoMakie.heatmap!(ax1, t_range ./ days, z_range, PAR_values; colormap=:viridis) -Colorbar(fig_forcing[1, 2], hm1) - -ax2 = Axis(fig_forcing[2, 1]; xlabel="Time (days)", ylabel="Depth (m)", title="diffusivity") -hm2 = CairoMakie.heatmap!(ax2, t_range ./ days, z_range, κₜ_values; colormap=:viridis) -Colorbar(fig_forcing[2, 2], hm2) - -display(fig_forcing) -fig_forcing - -# ## Physical model - -grid = RectilinearGrid(; size=(1, 1, 20), extent=(20meters, 20meters, 200meters)) -nothing #hide - -bgc_model = Biogeochemistry( - bgc; light_attenuation=FunctionFieldPAR(; grid, PAR_f=constant_PAR) -) -nothing #hide - -full_model = NonhydrostaticModel(; - grid, - clock=Clock(; time=0.0), - timestepper=:QuasiAdamsBashforth2, - closure=ScalarDiffusivity( - VerticallyImplicitTimeDiscretization(); ν=diffusivity_profile, κ=diffusivity_profile - ), - biogeochemistry=bgc_model, -) -nothing #hide - -# ## Initial conditions - -set!(full_model; default_initial_conditions(bgc; detritus = 0.0, total_plankton_biomass = 0.12)...) # mmol N / m³ - -# ## Simulation -filename = "N2P2ZD_column.jld2" - -simulation = Simulation(full_model; Δt=1hours, stop_time=1year) - -simulation.output_writers[:profiles] = JLD2Writer( - full_model, - full_model.tracers; - filename=filename, - schedule=TimeInterval(1day), - overwrite_existing=true, -) - -run!(simulation) -nothing #hide - -# ## Plotting - -#Load time series data -timeseries = NamedTuple{keys(full_model.tracers)}( - FieldTimeSeries(filename, "$field") for field in keys(full_model.tracers) -) - -# Use Agate's introspection helpers to recover the structural tracer layout -all_keys = [groups.plankton..., groups.nonplankton...] -nothing #hide - -#Create figure with appropriate size -fig = Figure(; size=(800, 1200), fontsize=16) - -#Plot all fields -for (i, key) in enumerate(all_keys) - x_nodes, y_nodes, z_nodes = nodes(timeseries[key]) - z_vals = collect(z_nodes) - times = collect(timeseries[key].times / days) - - ax = Axis( - fig[i, 1]; - title="$(key) concentration (mmol N / m³)", - xlabel="Time (days)", - ylabel="z (m)", - limits=((0, 365), (-200, 0)), - ) - hm = heatmap!( - ax, - times, - z_vals, - Float32.(interior(timeseries[key],1,1,:,:)'); - colormap=:viridis, - rasterize=true, - ) # Rasterize for smaller output - Colorbar(fig[i, 2], hm) -end - -#Save figure -save("N2P2ZD_column.png", fig) - -display(fig) -fig # Display the figure - -# Plot the final-time depth-bin value of every tracer as horizontal bars. -n_profiles = length(all_keys) -n_columns = min(3, n_profiles) -n_rows = cld(n_profiles, n_columns) -fig_profiles = Figure(; size=(350 * n_columns, 300 * n_rows), fontsize=16) - -profile_axes = Axis[] -for (i, key) in enumerate(all_keys) - row = cld(i, n_columns) - column = mod1(i, n_columns) - x_nodes, y_nodes, z_nodes = nodes(timeseries[key]) - z_centers = collect(z_nodes) - final_profile = vec(interior(timeseries[key], 1, 1, :, length(timeseries[key].times))) - - z_edges = similar(z_centers, length(z_centers) + 1) - z_edges[2:end-1] .= (z_centers[1:end-1] .+ z_centers[2:end]) ./ 2 - z_edges[1] = z_centers[1] - (z_centers[2] - z_centers[1]) / 2 - z_edges[end] = z_centers[end] + (z_centers[end] - z_centers[end-1]) / 2 - - ax = Axis( - fig_profiles[row, column]; - title=String(key), - xlabel="Concentration (mmol N / m³)", - ylabel=column == 1 ? "z (m)" : "", - limits=(nothing, (-200, 0)), - ) - - for j in eachindex(final_profile) - y0 = z_edges[j] - y1 = z_edges[j + 1] - x0 = min(0, final_profile[j]) - width = abs(final_profile[j]) - poly!(ax, Rect(x0, y0, width, y1 - y0); color=:dodgerblue, strokecolor=:dodgerblue) - end - - push!(profile_axes, ax) -end - -linkyaxes!(profile_axes...) -for (i, ax) in enumerate(profile_axes) - column = mod1(i, n_columns) - - if column > 1 - hideydecorations!(ax; grid=false) - end -end -save("N2P2ZD_column_final_profiles.png", fig_profiles) - -display(fig_profiles) -fig_profiles diff --git a/2026-05/examples/12_diffusivity_seasonal.jl b/2026-05/examples/11_diffusivity_seasonal.jl similarity index 100% rename from 2026-05/examples/12_diffusivity_seasonal.jl rename to 2026-05/examples/11_diffusivity_seasonal.jl