diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index d3ebf43..e45fe5c 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -28,7 +28,10 @@ jobs: - name: Set up Julia uses: julia-actions/setup-julia@v2 with: - version: '1' + version: '1.10' + + - name: Set up Quarto + uses: quarto-dev/quarto-actions/setup@v2 - name: Cache Julia artifacts uses: julia-actions/cache@v2 @@ -39,6 +42,9 @@ jobs: using Pkg Pkg.instantiate() + - name: Render Quarto slides + run: quarto render 2026-05/slides + - name: Build and deploy documentation run: julia --color=yes --project=2026-05/docs 2026-05/docs/make.jl env: diff --git a/.gitignore b/.gitignore index 0cc077a..1891b38 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,14 @@ 2026-05/docs/build/ 2026-05/docs/Manifest.toml 2026-05/docs/src/generated +2026-05/docs/src/slides + +# Quarto slides +2026-05/slides/01_intro_files/ +2026-05/slides/_site/ +2026-05/slides/.quarto + +*.html # Model output *.jld @@ -9,6 +17,5 @@ # Plots *.png -*.svg *.gif *.mp4 \ No newline at end of file diff --git a/2026-05/advances_excercises/06_diffusivity_sensitivity.jl b/2026-05/advances_excercises/06_diffusivity_sensitivity.jl new file mode 100644 index 0000000..8c548cc --- /dev/null +++ b/2026-05/advances_excercises/06_diffusivity_sensitivity.jl @@ -0,0 +1,158 @@ +# # [Exercise 06 supplement: Diffusivity sensitivity] (@id diffusivity_sensitivity_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. + +# In this example we run a default Agate.jl-NiPiZD model inside a 2 layer column model. +# The simulation is repeated across a range of vertical diffusivities, then the top-layer +# total plankton biomass is summarized as a function of κ_max. + +# ## Loading dependencies + +using Agate +using Agate.Library.Light +using OceanBioME +using OceanBioME: Biogeochemistry +using Oceananigans +using Oceananigans.Units +using CairoMakie +using Statistics + +stop_time = 3*365day +Δt = 1hour +output_interval = 1day +LAYER_INTERFACE = -100meters +PAR_SURFACE_MAX = 80 + +κ_values = [ + 3e-5, + 1e-4, + 3e-4, + 1e-3, +] + +# ## Forcings + +struct LayeredDiffusivity{T} + κ_max::T + layer_interface::T +end + +@inline function (diffusivity::LayeredDiffusivity)(x, y, z, t) + return ifelse(z >= diffusivity.layer_interface, diffusivity.κ_max, zero(diffusivity.κ_max)) +end + +function surface_irradiance(x, y, z, t) + return ifelse(z >= LAYER_INTERFACE, PAR_SURFACE_MAX, 0.0) +end + +# ## Physical and ecosystem model + +function build_model(κ_max) + grid = RectilinearGrid(; size=(1, 1, 2), extent=(20meters, 20meters, 200meters)) + diffusivity = LayeredDiffusivity(κ_max, LAYER_INTERFACE) + + bgc = Agate.Models.NiPiZD.construct() + bgc_model = Biogeochemistry( + bgc; light_attenuation=FunctionFieldPAR(; grid, PAR_f=surface_irradiance) + ) + + model = NonhydrostaticModel(; + grid, + clock=Clock(; time=0.0), + timestepper=:QuasiAdamsBashforth2, + closure=ScalarDiffusivity( + VerticallyImplicitTimeDiscretization(); ν=diffusivity, κ=diffusivity + ), + biogeochemistry=bgc_model, + ) + + set!(model; N=7.0, P1=0.01, P2=0.01, Z1=0.05, Z2=0.05, D=0.0) + + return model +end + +function output_filename(κ_max) + κ_label = replace(string(κ_max), "." => "p", "-" => "m") + return "N2P2ZD_column_k$(κ_label).jld2" +end + +function run_column(κ_max) + model = build_model(κ_max) + filename = output_filename(κ_max) + + simulation = Simulation(model; Δt, stop_time) + + simulation.output_writers[:profiles] = JLD2Writer( + model, + model.tracers; + filename, + schedule=TimeInterval(output_interval), + overwrite_existing=true, + ) + + run!(simulation) + + return filename, model +end + +function top_layer_index(timeseries_field) + _, _, z_nodes = nodes(timeseries_field) + z_vals = collect(z_nodes) + return argmax(z_vals) +end + +function top_layer_plankton_biomass(filename) + P1 = FieldTimeSeries(filename, "P1") + P2 = FieldTimeSeries(filename, "P2") + Z1 = FieldTimeSeries(filename, "Z1") + Z2 = FieldTimeSeries(filename, "Z2") + + k_top = top_layer_index(P1) + times = collect(P1.times ./ days) + + biomass = vec( + interior(P1, 1, 1, k_top, :) .+ + interior(P2, 1, 1, k_top, :) .+ + interior(Z1, 1, 1, k_top, :) .+ + interior(Z2, 1, 1, k_top, :) + ) + + final_biomass = biomass[end] + final_year = times .>= maximum(times) - 365 + mean_final_year_biomass = mean(biomass[final_year]) + + return final_biomass, mean_final_year_biomass +end + +# ## Sensitivity runs + +final_top_layer_biomass = Float64[] +mean_final_year_top_layer_biomass = Float64[] + +for κ_max in κ_values + @info "Running diffusivity sensitivity" κ_max + filename, _ = run_column(κ_max) + final_biomass, mean_final_year_biomass = top_layer_plankton_biomass(filename) + push!(final_top_layer_biomass, final_biomass) + push!(mean_final_year_top_layer_biomass, mean_final_year_biomass) +end + +# ## Plot top-layer total plankton biomass against κ_max + +fig = Figure(; size=(850, 550), fontsize=16) +ax = Axis( + fig[1, 1]; + xlabel="κ_max (m² s⁻¹)", + ylabel="Top-layer total plankton biomass (mmol N m⁻³)", + title="Diffusivity sensitivity", + xscale=log10, +) + +scatterlines!(ax, κ_values, mean_final_year_top_layer_biomass; label="Final-year mean") +axislegend(ax; position=:rb) + +save("top_layer_plankton_biomass_vs_kmax.png", fig) + +fig diff --git a/2026-05/docs/Project.toml b/2026-05/docs/Project.toml index 3ef2d1a..b91e2f2 100644 --- a/2026-05/docs/Project.toml +++ b/2026-05/docs/Project.toml @@ -10,7 +10,7 @@ OceanBioME = "a49af516-9db8-4be4-be45-1dad61c5a376" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" [compat] -Agate = "0.5.2" +Agate = "0.5.3" CairoMakie = "0.12, 0.13, 0.14, 0.15" CSV = "0.10" DataFrames = "1" diff --git a/2026-05/docs/make.jl b/2026-05/docs/make.jl index 2c4802c..2968567 100644 --- a/2026-05/docs/make.jl +++ b/2026-05/docs/make.jl @@ -1,21 +1,42 @@ using Pkg Pkg.activate(@__DIR__) + +workshop_root = normpath(joinpath(@__DIR__, "..")) Pkg.instantiate() using Literate using Documenter -workshop_root = normpath(joinpath(@__DIR__, "..")) - examples_dir = joinpath(workshop_root, "examples") docs_src_dir = joinpath(@__DIR__, "src") -docs_exercises_dir = joinpath(docs_src_dir, "exercises") +docs_exercises_dir = joinpath(docs_src_dir, "generated") mkpath(docs_exercises_dir) + + example_files = [ "00_setup_check.jl", + "01_quick_start.jl", + "02_diagnostics.jl", + "03_allometric_scaling.jl", + "04_number_size_classes.jl", + "05_size_range.jl", + "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", ] +function strip_jld2_warnings(content) + return replace( + content, + r"(?ms)^┌ Warning:.*?^└ @ JLD2 .*?/(writing_datatypes|reconstructing_datatypes)\.jl:\d+\n" => "", + ) +end + # Run from the workshop root so relative paths in examples are predictable. cd(workshop_root) do mkpath("figures") @@ -25,37 +46,72 @@ cd(workshop_root) do source = joinpath(examples_dir, file) # Generate and execute the rendered documentation page. - # - # Executing here lets CI pre-run the examples and include generated figures + # Executing here lets CI pre-run examples and include generated figures # in the deployed documentation. Literate.markdown( source, docs_exercises_dir; - documenter = true, - execute = true, - credit = false, + documenter=true, + execute=true, + credit=false, + postprocess=strip_jld2_warnings, ) end end +# Quarto renders slide decks into slides/_site before this script runs in CI. +# Copy that site into docs/src before makedocs so Documenter can validate +# local links to the rendered slide HTML. +slides_site_src = joinpath(workshop_root, "slides", "_site") +slides_docs_src = joinpath(docs_src_dir, "slides") + +if isdir(slides_site_src) + isdir(slides_docs_src) && rm(slides_docs_src; recursive=true, force=true) + cp(slides_site_src, slides_docs_src) +else + @warn "Quarto slide output not found before makedocs; slide links may be unavailable" slides_site_src +end + makedocs( - sitename = "Agate.jl workshop 2026-05", - format = Documenter.HTML( - prettyurls = get(ENV, "CI", "false") == "true", - assets = String[], + sitename="Agate.jl workshop 2026-05", + format=Documenter.HTML( + prettyurls=get(ENV, "CI", "false") == "true", + assets=String[], + size_threshold_warn=1_000_000, + size_threshold=2_500_000, ), - modules = Module[], - pages = [ + modules=Module[], + pages=[ "Home" => "index.md", "Setup" => "setup.md", + "Workshop slides" => "lectures.md", "Examples" => [ - "00 Setup Check" => "exercises/00_setup_check.md", + "00 Setup Check" => "generated/00_setup_check.md", + "01 Quick Start" => "generated/01_quick_start.md", + "02 Diagnostics" => "generated/02_diagnostics.md", + "03 Allometric Scaling" => "generated/03_allometric_scaling.md", + "04 Number of Size Classes" => "generated/04_number_size_classes.md", + "05 Size Range" => "generated/05_size_range.md", + "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", ], ], ) +# Ensure rendered slides are present in the final build output even if a +# future Documenter version changes how non-markdown files are copied. +slides_build_dst = joinpath(@__DIR__, "build", "slides") +if isdir(slides_site_src) + isdir(slides_build_dst) && rm(slides_build_dst; recursive=true, force=true) + cp(slides_site_src, slides_build_dst) +end deploydocs( - repo = "github.com/agate-model/workshops.git", - devbranch = "main", - push_preview = true, + repo="github.com/agate-model/workshops.git", + devbranch="main", + push_preview=true, ) diff --git a/2026-05/docs/src/advanced.md b/2026-05/docs/src/advanced.md new file mode 100644 index 0000000..2ed97c6 --- /dev/null +++ b/2026-05/docs/src/advanced.md @@ -0,0 +1,8 @@ +| Example | Theory | +|---|---| +| sinking + remin | POC attenuation | +| exponential mortality | viral shunt | +| kN vs umax | r/K trade-off | +| specificity vs gmax | generalist-specialist trade-off | +| alpha vs umax | LL vs HL prochlorococcus | +| interaction roles | cannibalism | diff --git a/2026-05/docs/src/exercises/00_setup_check.md b/2026-05/docs/src/exercises/00_setup_check.md deleted file mode 100644 index c3d1915..0000000 --- a/2026-05/docs/src/exercises/00_setup_check.md +++ /dev/null @@ -1,23 +0,0 @@ -```@meta -EditURL = "../../../examples/00_setup_check.jl" -``` - -# Setup check - -Run this script to check that the Docker environment is setup correctly. - -````julia -using Agate - -println("Agate.jl loaded successfully.") -println("Active project: ", Base.active_project()) -```` - -```` -Agate.jl loaded successfully. -Active project: /home/phyto/workshops/2026-05/docs/Project.toml - -```` - -The setup check is complete if this file runs without errors. - diff --git a/2026-05/docs/src/index.md b/2026-05/docs/src/index.md index c5e688e..3411c67 100644 --- a/2026-05/docs/src/index.md +++ b/2026-05/docs/src/index.md @@ -7,11 +7,44 @@ Welcome to the Agate.jl workshop. **Before the workshop** please follow the [setup instructions](setup.md) to setup the Docker environment and download the workshop materials. If you are unable to install Docker please contact the workshop organizers. + +## Schedule + +| Time | Session | +|---|---| +| 09:30–10:00 | *Docker setup (optional, this can be done before the workshop)* | +| 10:00–10:30 | Welcome and workshop overview | +| 10:30–11:10 | Session 1: Size, allometry, and plankton traits | +| 11:10–11:25 | Break | +| 11:25–12:05 | Session 2: Predation and trophic structure | +| 12:05–12:45 | Session 3: Physical forcing: light and diffusivity | +| 12:45–13:45 | Lunch | +| 13:45–14:35 | Group exercise | +| 14:35–15:10 | Own-work block 1 | +| 15:10–15:40 | Break | +| 15:40–16:30 | Own-work block 2 | +| 16:30–16:55 | Show-and-tell (optional) and feedback | +| 16:55–17:00 | Wrap-up | +| 17:00 | Pub | + + ## Examples ```@contents Pages = [ "exercises/00_setup_check.md", + "exercises/01_quick_start.md", + "exercises/02_diagnostics.md", + "exercises/03_allometric_scaling.md", + "exercises/04_number_and_size.md", + "exercises/05_palatability.md", + "exercises/06_diffusivity.md", + "exercises/07_irradiance.md", ] Depth = 2 ``` + +## Workshop slides + +The workshop lecture slides are available from the [Workshop slides](lectures.md) page. + diff --git a/2026-05/docs/src/lectures.md b/2026-05/docs/src/lectures.md new file mode 100644 index 0000000..4357029 --- /dev/null +++ b/2026-05/docs/src/lectures.md @@ -0,0 +1,3 @@ +# Workshop slides + +- [Morning overview](slides/01_intro.html) diff --git a/2026-05/docs/src/src/WorkshopDiagnostics.jl b/2026-05/docs/src/src/WorkshopDiagnostics.jl new file mode 100644 index 0000000..5623ecf --- /dev/null +++ b/2026-05/docs/src/src/WorkshopDiagnostics.jl @@ -0,0 +1 @@ +include(joinpath(@__DIR__, "..", "..", "..", "src", "WorkshopDiagnostics.jl")) diff --git a/2026-05/docs/src/src/WorkshopSetup.jl b/2026-05/docs/src/src/WorkshopSetup.jl new file mode 100644 index 0000000..8b49318 --- /dev/null +++ b/2026-05/docs/src/src/WorkshopSetup.jl @@ -0,0 +1 @@ +include(joinpath(@__DIR__, "..", "..", "..", "src", "WorkshopSetup.jl")) diff --git a/2026-05/examples/01_quick_start.jl b/2026-05/examples/01_quick_start.jl new file mode 100644 index 0000000..25a99a7 --- /dev/null +++ b/2026-05/examples/01_quick_start.jl @@ -0,0 +1,146 @@ +# # [Exercise 01: Quick start] (@id quick_start_exercise) +# +# This exercise constructs an Agate biogeochemistry model, +# wraps it in an OceanBioME `Biogeochemistry`, places it in an Oceananigans +# `BoxModel`, sets initial tracer values, and runs a simulation. + +# ## Loading dependencies + +using Agate +using Agate.Introspection: tracer_names +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) + + +mkpath(joinpath("outputs")) +mkpath(joinpath("figures")) + +nothing #hide + +# ## Agate.jl Quick start, step by step +# +# First construct the default two-phytoplankton, two-zooplankton NiPiZD model. + +bgc = Agate.Models.NiPiZD.construct() +println(tracer_names(bgc)) + +nothing #hide + +# The Quick start uses a seasonal photosynthetically active radiation field for a +# zero-dimensional box. + +light_attenuation = FunctionFieldPAR(; grid=BoxModelGrid()) + +nothing #hide + +# OceanBioME adapts the Agate biogeochemistry to an Oceananigans model. + +bgc_model = Biogeochemistry(bgc; light_attenuation) +full_model = BoxModel(; biogeochemistry=bgc_model) + +nothing #hide + +# Set the initial tracer concentrations. The helper distributes the total plankton biomass evenly across the plankton tracers in `bgc`. + +set!(full_model; default_initial_conditions(bgc)...) + +nothing #hide + +# Run the model and save daily tracer output. + +quickstart_filename = joinpath("outputs", "01_quick_start_manual.jld2") +simulation = Simulation(full_model; Δt=240minutes, stop_time=1095days) + +simulation.output_writers[:fields] = JLD2Writer( + full_model, + full_model.fields; + filename=quickstart_filename, + schedule=TimeInterval(1day), + overwrite_existing=true, +) + +run!(simulation) + +nothing #hide + +# ## Reading and plotting the output + +tracer_syms = tracer_names(bgc) +timeseries = read_box_tracer_timeseries(quickstart_filename, tracer_syms) + +fig_manual = Figure(; size=(1200, 800), fontsize=20) + +for (idx, sym) in enumerate(tracer_syms) + row = floor(Int, (idx - 1) / 2) + 1 + col = Int((idx - 1) % 2) + 1 + ax = Axis( + fig_manual[row, col]; + ylabel=string(sym), + xlabel="Days", + title="$(sym) concentration (mmol N m⁻³)", + ) + lines!(ax, timeseries.times, timeseries.data[Symbol(sym)]; linewidth=3) +end + +save(joinpath("figures", "01_quick_start_manual.png"), fig_manual) +display(fig_manual) +fig_manual + +# ## A general workshop box-model wrapper +# +# The same setup will appear repeatedly in later exercises. The workshop +# helper script therefore defines: +# +# - `default_initial_conditions(bgc; total_plankton_biomass, nutrient, detritus)` +# - `build_box_model(bgc; light_attenuation, initial_conditions)` +# - `run_box_model(bgc; filename, initial_conditions, Δt, stop_time, output_interval)` +# - `read_box_tracer_timeseries(filename, tracer_syms)` +# +# The wrapper takes `bgc` as an argument, so it can be reused with more complex +# Agate models without manually spelling out every plankton tracer. + +wrapped = run_box_model( + bgc; + filename=joinpath("outputs", "01_quick_start_wrapped.jld2"), + Δt=240minutes, + stop_time=1095days, + output_interval=1day, +) + +wrapped_timeseries = read_box_tracer_timeseries(wrapped.filename, wrapped.tracer_syms) + +nothing #hide + +# The wrapper is also useful when we change model complexity. Here we construct a +# larger NiPiZD community and use the same `run_box_model` function. + +larger_bgc = Agate.Models.NiPiZD.construct( + phyto_size_structure=(n=3, min_esd=1.0, max_esd=10.0, splitting=:log_splitting), + zoo_size_structure=[10.0, 32.0, 100.0], +) + +larger = run_box_model( + larger_bgc; + filename=joinpath("outputs", "01_quick_start_larger_community.jld2"), + initial_conditions=default_initial_conditions(larger_bgc; total_plankton_biomass=0.06), + Δt=240minutes, + stop_time=365days, + output_interval=1day, +) + +println(larger.tracer_syms) + +nothing #hide diff --git a/2026-05/examples/02_diagnostics.jl b/2026-05/examples/02_diagnostics.jl new file mode 100644 index 0000000..1e4ea9f --- /dev/null +++ b/2026-05/examples/02_diagnostics.jl @@ -0,0 +1,135 @@ +# # [Exercise 02: Diagnostics] (@id diagnostics_exercise) +# +# This exercise introduces reusable diagnostics from the workshop helper script so later +# exercises can load and reuse them directly. + +# ## Loading dependencies + +using CairoMakie +using Agate +using Oceananigans.Units: day +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) + + +mkpath(joinpath("outputs")) +mkpath(joinpath("figures")) + +nothing #hide + +# ## Run the Quick start model +# +# The diagnostics operate on saved box-model output. We reuse the wrapper from +# Exercise 01 to create a small output file. + +bgc = default_quickstart_bgc() +run = run_box_model(bgc; filename=joinpath("outputs", "02_diagnostics_quick_start.jld2")) +timeseries = read_box_tracer_timeseries(run.filename, run.tracer_syms) + +times = timeseries.times +data = timeseries.data + +nothing #hide + +# ## Tracer concentrations +# +# `plot_box_timeseries` accepts the complete named tuple returned by +# `read_box_tracer_timeseries` and plots each tracer in a separate panel. + +tracer_concentrations_figure_path = joinpath("figures", "02_diagnostic_tracer_concentrations.png") +fig_tracers = plot_box_timeseries(timeseries) +save(tracer_concentrations_figure_path, fig_tracers; px_per_unit=1) +display(fig_tracers) +fig_tracers + +# ### Comparison +# +# The same helper can also compare two or more box-model time series. The +# default detritus remineralization rate is `0.1213 / day`. Here we compare +# the default run with a higher remineralization case, `0.25 / day`, for all +# shared tracers. + +high_remineralization_bgc = Agate.Models.NiPiZD.construct(; + parameters = (detritus_remineralization = 0.25 / day,), +) +high_remineralization_run = run_box_model( + high_remineralization_bgc; + filename=joinpath("outputs", "02_diagnostics_high_detritus_remineralization.jld2"), +) +high_remineralization_timeseries = read_box_tracer_timeseries( + high_remineralization_run.filename, + high_remineralization_run.tracer_syms, +) + +comparison_figure_path = joinpath("figures", "02_diagnostic_timeseries_comparison.png") +fig_comparison = plot_box_timeseries( + [timeseries, high_remineralization_timeseries]; + labels=["default", "detritus remineralization = 0.25 / day"], +) +save(comparison_figure_path, fig_comparison; px_per_unit=1) +display(fig_comparison) +fig_comparison + +# ## Relative nitrogen contributions +# +# `plot_contributions` shows two coordinated stacked-area diagnostics: +# living versus non-living nitrogen, and phytoplankton versus zooplankton +# contributions to living biomass. + +nitrogen_contributions_figure_path = joinpath("figures", "02_diagnostic_relative_nitrogen_contributions.png") +fig_nitrogen = plot_contributions(times, data) +save(nitrogen_contributions_figure_path, fig_nitrogen; px_per_unit=1) +display(fig_nitrogen) +fig_nitrogen + +# ## Community-weighted mean size +# +# Agate.jl stores plankton equivalent spherical diameter (ESD) metadata on +# constructed biogeochemistry objects. `plot_cwm_size` accepts a box-model time +# series and the matching biogeochemistry object, then computes the +# community-weighted mean size for all plankton, phytoplankton, and zooplankton. + +cwm_size_figure_path = joinpath("figures", "02_diagnostic_cwm_size.png") +fig_size = plot_cwm_size(timeseries, bgc) +save(cwm_size_figure_path, fig_size; px_per_unit=1) +display(fig_size) +fig_size + +# ## Community-weighted mean size comparison +# +# The same helper can compare multiple runs. Each time series is paired with a +# biogeochemistry object because the plankton sizes are defined by the model, +# and different models may have different plankton groups or size structures. + +cwm_size_comparison_figure_path = joinpath("figures", "02_diagnostic_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size( + [timeseries, high_remineralization_timeseries], + [bgc, high_remineralization_bgc]; + labels=["default", "detritus remineralization = 0.25 / day"], +) +save(cwm_size_comparison_figure_path, fig_size_comparison; px_per_unit=1) +display(fig_size_comparison) +fig_size_comparison + +# ## Parameter bars +# +# The same helper can also inspect model parameters directly. Here we plot the default +# maximum phytoplankton growth-rate values for the plankton types represented in the model. + +parameter_bar_figure_path = joinpath("figures", "02_diagnostic_mumax_parameter_bars.png") +fig_mumax_parameter = plot_plankton_parameter_bars( + bgc, + :maximum_growth_rate; + ylabel = "maximum_growth_rate", + title = "Maximum growth rate by phytoplankton type", +) +save(parameter_bar_figure_path, fig_mumax_parameter; px_per_unit=1) +display(fig_mumax_parameter) +fig_mumax_parameter diff --git a/2026-05/examples/03_allometric_scaling.jl b/2026-05/examples/03_allometric_scaling.jl new file mode 100644 index 0000000..d7b129b --- /dev/null +++ b/2026-05/examples/03_allometric_scaling.jl @@ -0,0 +1,200 @@ +# # [Exercise 03: Allometric scaling] (@id allometric_scaling_exercise) +# +# This exercise changes allometric scaling in the Agate.jl NiPiZD model. +# We inspect the plankton parameter values produced by different allometric coefficients, then run each case in a well-mixed zero-dimensional box. +# +# Agate.jl represents allometric parameter rules as a power law on spherical cell volume: +# +# ```math +# \mathrm{trait} = a V^b, \qquad V = \frac{4}{3}\pi\left(\frac{d}{2}\right)^3 +# ``` +# +# where `d` is equivalent spherical diameter (ESD), `a` is the `prefactor`, and `b` is the `exponent`. +# Changing `a` shifts the whole curve up or down; changing `b` changes how strongly the trait depends on plankton size. +# +# In this exercise we focus on three size-dependent parameters. +# `mumax` is the maximum phytoplankton growth rate. +# `kN` is the nutrient half-saturation concentration: lower values mean growth saturates at lower nutrient concentration. +# `gmax` is the maximum zooplankton predation rate. + +# ## Loading dependencies +# +# The example uses Agate.jl, Oceananigans.jl, and OceanBioME.jl for the ecosystem simulation. +# CairoMakie.jl is used for plotting. + +using Agate +using Agate.Introspection: plankton_groups +using Agate.Library.Allometry: AllometricParam, PowerLaw +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) + +mkpath("outputs") +mkpath("figures") + +nothing #hide + +# ## Allometry cases +# +# We start with the default allometric coefficients, then construct two alternatives: +# a flat case where the selected parameters do not vary with plankton size, and a +# stronger small-fast case where smaller plankton have higher growth and predation rates. + +bgc_default = Agate.Models.NiPiZD.construct() + +println("Plankton groups: ", plankton_groups(bgc_default)) + +bgc_flat = Agate.Models.NiPiZD.construct(; + parameters = ( + maximum_growth_rate = AllometricParam(PowerLaw(); prefactor = 2 / day, exponent = 0.0), + nutrient_half_saturation = AllometricParam(PowerLaw(); prefactor = 0.17, exponent = 0.0), + maximum_predation_rate = AllometricParam(PowerLaw(); prefactor = 30.84 / day, exponent = 0.0), + ), +) + +bgc_strong_small_fast = Agate.Models.NiPiZD.construct(; + parameters = ( + maximum_growth_rate = AllometricParam(PowerLaw(); prefactor = 2 / day, exponent = -0.35), + nutrient_half_saturation = AllometricParam(PowerLaw(); prefactor = 0.17, exponent = 0.35), + maximum_predation_rate = AllometricParam(PowerLaw(); prefactor = 30.84 / day, exponent = -0.35), + ), +) + +bgc_cases = [bgc_default, bgc_flat, bgc_strong_small_fast] +bgc_case_labels = ["Default", "Flat", "Strong small-fast"] + +# ## Parameter bar charts +# +# `plot_plankton_parameter_bars` accepts one model or an array of models. When an +# array is provided, it draws grouped bars so the same plankton type can be compared +# across parameter sets. + +fig_mumax = plot_plankton_parameter_bars( + bgc_cases, + :maximum_growth_rate; + labels = bgc_case_labels, + ylabel = "maximum_growth_rate", + title = "Maximum growth rate by phytoplankton type", + figure_path = joinpath("figures", "03_mumax_parameter_bars.png"), +) +display(fig_mumax) +fig_mumax + +fig_kN = plot_plankton_parameter_bars( + bgc_cases, + :nutrient_half_saturation; + labels = bgc_case_labels, + ylabel = "nutrient_half_saturation", + title = "Nutrient half-saturation by phytoplankton type", + figure_path = joinpath("figures", "03_kN_parameter_bars.png"), +) +display(fig_kN) +fig_kN + +fig_gmax = plot_plankton_parameter_bars( + bgc_cases, + :maximum_predation_rate; + labels = bgc_case_labels, + ylabel = "maximum_predation_rate", + title = "Maximum predation rate by zooplankton type", + figure_path = joinpath("figures", "03_gmax_parameter_bars.png"), +) +display(fig_gmax) +fig_gmax + +# ## Zero-dimensional ecosystem simulations +# +# The parameter plots show potential rates. +# We now run the same three allometric cases in a well-mixed box model and compare the ecosystem dynamics. + +# Run the default case. +default_run = run_box_model( + bgc_default; + filename = joinpath("outputs", "03_default.jld2"), + initial_conditions = default_initial_conditions(bgc_default; nutrient = 8.0, total_plankton_biomass = 0.08), +) +default_filename = default_run.filename + +# Run the flat-allometry case. +flat_run = run_box_model( + bgc_flat; + filename = joinpath("outputs", "03_flat.jld2"), + initial_conditions = default_initial_conditions(bgc_flat; nutrient = 8.0, total_plankton_biomass = 0.08), +) +flat_filename = flat_run.filename + +# Run the strong small-fast case. +strong_run = run_box_model( + bgc_strong_small_fast; + filename = joinpath("outputs", "03_strong_small_fast.jld2"), + initial_conditions = default_initial_conditions(bgc_strong_small_fast; nutrient = 8.0, total_plankton_biomass = 0.08), +) +strong_filename = strong_run.filename + +nothing #hide + +# ## Diagnostic plots +# +# The simulations are compared with the reusable diagnostics introduced in Exercise 02. +# `plot_box_timeseries` compares tracer concentrations, `plot_contributions` shows +# relative nitrogen partitioning for each allometric configuration, and +# `plot_cwm_size` compares community-weighted mean plankton sizes. + +default_timeseries = read_box_tracer_timeseries(default_run.filename, default_run.tracer_syms) +flat_timeseries = read_box_tracer_timeseries(flat_run.filename, flat_run.tracer_syms) +strong_timeseries = read_box_tracer_timeseries(strong_run.filename, strong_run.tracer_syms) + +comparison_figure_path = joinpath("figures", "03_allometry_timeseries_comparison.png") +fig_comparison = plot_box_timeseries( + [default_timeseries, flat_timeseries, strong_timeseries]; + labels = bgc_case_labels, +) +save(comparison_figure_path, fig_comparison; px_per_unit = 1) +display(fig_comparison) +fig_comparison + +# ### Relative nitrogen contributions: default allometry + +nitrogen_default_figure_path = joinpath("figures", "03_allometry_relative_nitrogen_default.png") +fig_nitrogen_default = plot_contributions(default_timeseries.times, default_timeseries.data) +save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1) +display(fig_nitrogen_default) +fig_nitrogen_default + +# ### Relative nitrogen contributions: flat allometry + +nitrogen_flat_figure_path = joinpath("figures", "03_allometry_relative_nitrogen_flat.png") +fig_nitrogen_flat = plot_contributions(flat_timeseries.times, flat_timeseries.data) +save(nitrogen_flat_figure_path, fig_nitrogen_flat; px_per_unit = 1) +display(fig_nitrogen_flat) +fig_nitrogen_flat + +# ### Relative nitrogen contributions: strong small-fast allometry + +nitrogen_strong_figure_path = joinpath("figures", "03_allometry_relative_nitrogen_strong_small_fast.png") +fig_nitrogen_strong = plot_contributions(strong_timeseries.times, strong_timeseries.data) +save(nitrogen_strong_figure_path, fig_nitrogen_strong; px_per_unit = 1) +display(fig_nitrogen_strong) +fig_nitrogen_strong + +# ### Community-weighted mean size comparison + +cwm_size_comparison_figure_path = joinpath("figures", "03_allometry_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size( + [default_timeseries, flat_timeseries, strong_timeseries], + bgc_cases; + labels = bgc_case_labels, +) +save(cwm_size_comparison_figure_path, fig_size_comparison; px_per_unit = 1) +display(fig_size_comparison) +fig_size_comparison diff --git a/2026-05/examples/04_number_size_classes.jl b/2026-05/examples/04_number_size_classes.jl new file mode 100644 index 0000000..50f6bfb --- /dev/null +++ b/2026-05/examples/04_number_size_classes.jl @@ -0,0 +1,122 @@ +# # [Exercise 04: Number size classes] (@id number_size_classes_exercise) +# +# This exercise changes the number of phytoplankton and zooplankton size +# classes in the Agate.jl NiPiZD model. The total initial plankton biomass is +# held fixed and split evenly across the available plankton tracers, so the +# comparison isolates the effect of resolving more size classes. + +# ## Loading dependencies + +using Agate +using Agate.Introspection: tracer_names +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) + +mkpath("outputs") +mkpath("figures") + +nothing #hide + +# ## Number-of-classes cases +# +# The default model has two phytoplankton classes and two zooplankton classes. +# The alternative cases keep the same broad phytoplankton and zooplankton size +# ranges and logarithmic spacing, but increase both groups to five and ten +# classes. + +default_phyto_size_structure = (n = 2, min_esd = 2, max_esd = 10, splitting = :log_splitting) +default_zoo_size_structure = (n = 2, min_esd = 20, max_esd = 100, splitting = :linear_splitting) + +function construct_size_class_bgc(n) + return Agate.Models.NiPiZD.construct(; + phyto_size_structure = (; default_phyto_size_structure..., n), + zoo_size_structure = (; default_zoo_size_structure..., n), + ) +end + +bgc_default = Agate.Models.NiPiZD.construct(; + phyto_size_structure = default_phyto_size_structure, + zoo_size_structure = default_zoo_size_structure, +) +bgc_5_each = construct_size_class_bgc(5) +bgc_10_each = construct_size_class_bgc(10) + +bgcs = [bgc_default, bgc_5_each, bgc_10_each] +case_labels = ["default: 2 P, 2 Z", "5 P, 5 Z", "10 P, 10 Z"] + +for (label, bgc) in zip(case_labels, bgcs) + println(label) + println(tracer_names(bgc)) +end + +nothing #hide + +# ## Run zero-dimensional ecosystem simulations +# +# `default_initial_conditions` distributes the same total plankton biomass evenly +# across whichever plankton tracers are present in each configuration. + +runs = [ + run_box_model( + bgc; + filename = joinpath("outputs", "04_number_size_classes_$(i).jld2"), + initial_conditions = default_initial_conditions(bgc; nutrient = 8.0, total_plankton_biomass = 0.15), + ) + for (i, bgc) in enumerate(bgcs) +] + +timeseries = [read_box_tracer_timeseries(run.filename, run.tracer_syms) for run in runs] + +nothing #hide + +# ## Tracer concentration comparison +# +# This diagnostic compares the tracer set across the three configurations. Tracers that exist only in the higher-resolution cases appear only for those cases. + +comparison_variables = sort!(unique(vcat([collect(keys(ts.data)) for ts in timeseries]...)); by = string) +comparison_figure_path = joinpath("figures", "04_number_size_classes_timeseries_comparison.png") +fig_comparison = plot_box_timeseries(timeseries; labels = case_labels, variables = comparison_variables) +save(comparison_figure_path, fig_comparison; px_per_unit = 1) +display(fig_comparison) +fig_comparison + +# ## Relative nitrogen contributions: default size classes + +nitrogen_default_figure_path = joinpath("figures", "04_number_size_classes_relative_nitrogen_default.png") +fig_nitrogen_default = plot_contributions(timeseries[1].times, timeseries[1].data) +save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1) +display(fig_nitrogen_default) +fig_nitrogen_default + +# ## Relative nitrogen contributions: five classes each + +nitrogen_5_figure_path = joinpath("figures", "04_number_size_classes_relative_nitrogen_5_each.png") +fig_nitrogen_5 = plot_contributions(timeseries[2].times, timeseries[2].data) +save(nitrogen_5_figure_path, fig_nitrogen_5; px_per_unit = 1) +display(fig_nitrogen_5) +fig_nitrogen_5 + +# ## Relative nitrogen contributions: ten classes each + +nitrogen_10_figure_path = joinpath("figures", "04_number_size_classes_relative_nitrogen_10_each.png") +fig_nitrogen_10 = plot_contributions(timeseries[3].times, timeseries[3].data) +save(nitrogen_10_figure_path, fig_nitrogen_10; px_per_unit = 1) +display(fig_nitrogen_10) +fig_nitrogen_10 + +# ## Community-weighted mean size comparison + +size_comparison_figure_path = joinpath("figures", "04_number_size_classes_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size(timeseries, bgcs; labels = case_labels) +save(size_comparison_figure_path, fig_size_comparison; px_per_unit = 1) +display(fig_size_comparison) +fig_size_comparison diff --git a/2026-05/examples/05_size_range.jl b/2026-05/examples/05_size_range.jl new file mode 100644 index 0000000..83506a2 --- /dev/null +++ b/2026-05/examples/05_size_range.jl @@ -0,0 +1,99 @@ +# # [Exercise 05: Size range] (@id size_range_exercise) +# +# This exercise changes the size ranges represented by the Agate.jl NiPiZD +# model. We compare the default two-phytoplankton, two-zooplankton community +# with a wider but physiologically realistic range for both trophic groups. + +# ## Loading dependencies + +using Agate +using Agate.Introspection: tracer_names +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) + +mkpath("outputs") +mkpath("figures") + +nothing #hide + +# ## Size-range cases +# +# The wide-range case keeps two classes per trophic group and logarithmic spacing, +# but expands the endpoints to include small picophytoplankton through large +# microphytoplankton, and microzooplankton through small mesozooplankton. + +bgc_default = default_quickstart_bgc() + +bgc_wide_range = Agate.Models.NiPiZD.construct(; + phyto_size_structure = (n = 2, min_esd = 0.6, max_esd = 60.0, splitting = :log_splitting), + zoo_size_structure = (n = 2, min_esd = 6.0, max_esd = 600.0, splitting = :log_splitting), +) + +bgcs = [bgc_default, bgc_wide_range] +case_labels = ["default size range", "wide realistic size range"] + +for (label, bgc) in zip(case_labels, bgcs) + println(label) + println(tracer_names(bgc)) +end + +nothing #hide + +# ## Run zero-dimensional ecosystem simulations +# +# Total initial plankton biomass is the same in both cases and is split evenly +# across the plankton tracers. + +runs = [ + run_box_model( + bgc; + filename = joinpath("outputs", "05_size_range_$(i).jld2"), + initial_conditions = default_initial_conditions(bgc; nutrient = 8.0, total_plankton_biomass = 0.15), + ) + for (i, bgc) in enumerate(bgcs) +] + +timeseries = [read_box_tracer_timeseries(run.filename, run.tracer_syms) for run in runs] + +nothing #hide + +# ## Tracer concentration comparison + +comparison_figure_path = joinpath("figures", "05_size_range_timeseries_comparison.png") +fig_comparison = plot_box_timeseries(timeseries; labels = case_labels) +save(comparison_figure_path, fig_comparison; px_per_unit = 1) +display(fig_comparison) +fig_comparison + +# ## Relative nitrogen contributions: default size range + +nitrogen_default_figure_path = joinpath("figures", "05_size_range_relative_nitrogen_default.png") +fig_nitrogen_default = plot_contributions(timeseries[1].times, timeseries[1].data) +save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1) +display(fig_nitrogen_default) +fig_nitrogen_default + +# ## Relative nitrogen contributions: wide realistic size range + +nitrogen_wide_figure_path = joinpath("figures", "05_size_range_relative_nitrogen_wide.png") +fig_nitrogen_wide = plot_contributions(timeseries[2].times, timeseries[2].data) +save(nitrogen_wide_figure_path, fig_nitrogen_wide; px_per_unit = 1) +display(fig_nitrogen_wide) +fig_nitrogen_wide + +# ## Community-weighted mean size comparison + +size_comparison_figure_path = joinpath("figures", "05_size_range_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size(timeseries, bgcs; labels = case_labels) +save(size_comparison_figure_path, fig_size_comparison; px_per_unit = 1) +display(fig_size_comparison) +fig_size_comparison diff --git a/2026-05/examples/06_palatability.jl b/2026-05/examples/06_palatability.jl new file mode 100644 index 0000000..52f56c4 --- /dev/null +++ b/2026-05/examples/06_palatability.jl @@ -0,0 +1,171 @@ +# # [Exercise 06: Palatability] (@id palatability_exercise) +# +# Predator-prey palatability controls how strongly each zooplankton type grazes +# each phytoplankton type. This exercise follows the Agate.jl matrix example: +# first inspect the palatability matrices, then run each configuration in a box +# model. + +# ## Loading dependencies + +using Agate +using Agate.Introspection: interaction_matrix +using CairoMakie +using Oceananigans.Units: day +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) + +mkpath(joinpath("outputs")) +mkpath(joinpath("figures")) + +nothing #hide + +# ## Construct palatability configurations +# +# The default model derives its palatability matrix allometrically. By default, +# zooplankton prefer prey that are 10 times smaller than themselves. We compare +# that default with two alternatives: +# +# - a different preferred predator:prey size ratio, `Vopt = 5`, for both +# zooplankton types; +# - an explicit custom palatability matrix supplied directly. + +bgc = default_quickstart_bgc() + +vopt_bgc = Agate.Models.NiPiZD.construct(; + parameters = (optimum_predator_prey_ratio = (Z1 = 5.0, Z2 = 5.0),), +) + +custom_bgc = Agate.Models.NiPiZD.construct(; + palatability_matrix = [0.0 1.0; 1.0 0.0], +) + +bgcs = [bgc, vopt_bgc, custom_bgc] +case_labels = ["Vopt = 10", "Vopt = 5", "custom palatability"] + +nothing #hide + +# ## Palatability matrices +# +# `interaction_matrix(bgc, :palatability)` returns a labelled interaction table. +# Rows are consumers and columns are prey. + +function plot_interaction_matrix!(fig, position, table; title) + matrix = Matrix(table.matrix) + nrows, ncols = size(matrix) + + ax = Axis( + fig[position...]; + title, + xlabel = string(table.column_axis), + ylabel = string(table.row_axis), + xticks = (1:ncols, string.(table.columns)), + yticks = (1:nrows, string.(table.rows)), + yreversed = true, + aspect = DataAspect(), + ) + + hm = heatmap!(ax, 1:ncols, 1:nrows, matrix'; colorrange = (0, 1)) + + for row in 1:nrows, col in 1:ncols + text!( + ax, + col, + row; + text = string(round(matrix[row, col]; digits = 2)), + align = (:center, :center), + fontsize = 11, + ) + end + + return hm +end + +palatability_figure_path = joinpath("figures", "06_palatability_matrices.png") +fig_palatability = Figure(; size = (760, 280), fontsize = 12) + +palatability_tables = [interaction_matrix(bgc, :palatability) for bgc in bgcs] +hm = plot_interaction_matrix!(fig_palatability, (1, 1), palatability_tables[1]; title = case_labels[1]) +plot_interaction_matrix!(fig_palatability, (1, 2), palatability_tables[2]; title = case_labels[2]) +plot_interaction_matrix!(fig_palatability, (1, 3), palatability_tables[3]; title = case_labels[3]) +Colorbar(fig_palatability[1, 4], hm; label = "palatability") + +save(palatability_figure_path, fig_palatability; px_per_unit = 1) +display(fig_palatability) +fig_palatability + +# ## Run box models +# +# We use the same box-model wrapper as in Exercise 02 so the output can be read +# by the workshop diagnostics. + +run = run_box_model(bgc; filename = joinpath("outputs", "06_palatability_default.jld2")) +vopt_run = run_box_model(vopt_bgc; filename = joinpath("outputs", "06_palatability_vopt5.jld2")) +custom_run = run_box_model(custom_bgc; filename = joinpath("outputs", "06_palatability_custom.jld2")) + +timeseries = read_box_tracer_timeseries(run.filename, run.tracer_syms) +vopt_timeseries = read_box_tracer_timeseries(vopt_run.filename, vopt_run.tracer_syms) +custom_timeseries = read_box_tracer_timeseries(custom_run.filename, custom_run.tracer_syms) + +times = timeseries.times +data = timeseries.data + +nothing #hide + +# ## Tracer concentration comparison +# +# The concentration diagnostic compares all three palatability configurations. + +comparison_figure_path = joinpath("figures", "06_palatability_timeseries_comparison.png") +fig_comparison = plot_box_timeseries( + [timeseries, vopt_timeseries, custom_timeseries]; + labels = case_labels, +) +save(comparison_figure_path, fig_comparison; px_per_unit = 1) +display(fig_comparison) +fig_comparison + +# ## Relative nitrogen contributions: default palatability + +nitrogen_default_figure_path = joinpath("figures", "06_palatability_relative_nitrogen_default.png") +fig_nitrogen_default = plot_contributions(timeseries.times, timeseries.data) +save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1) +display(fig_nitrogen_default) +fig_nitrogen_default + +# ## Relative nitrogen contributions: Vopt = 5 + +nitrogen_vopt_figure_path = joinpath("figures", "06_palatability_relative_nitrogen_vopt5.png") +fig_nitrogen_vopt = plot_contributions(vopt_timeseries.times, vopt_timeseries.data) +save(nitrogen_vopt_figure_path, fig_nitrogen_vopt; px_per_unit = 1) +display(fig_nitrogen_vopt) +fig_nitrogen_vopt + +# ## Relative nitrogen contributions: custom palatability + +nitrogen_custom_figure_path = joinpath("figures", "06_palatability_relative_nitrogen_custom.png") +fig_nitrogen_custom = plot_contributions(custom_timeseries.times, custom_timeseries.data) +save(nitrogen_custom_figure_path, fig_nitrogen_custom; px_per_unit = 1) +display(fig_nitrogen_custom) +fig_nitrogen_custom + +# ## Community-weighted mean size comparison +# +# Each time series is paired with its matching biogeochemistry object because the +# plankton sizes are defined on the model configuration. + +cwm_size_comparison_figure_path = joinpath("figures", "06_palatability_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size( + [timeseries, vopt_timeseries, custom_timeseries], + bgcs; + labels = case_labels, +) +save(cwm_size_comparison_figure_path, fig_size_comparison; px_per_unit = 1) +display(fig_size_comparison) +fig_size_comparison diff --git a/2026-05/examples/07_assimilation_efficiency.jl b/2026-05/examples/07_assimilation_efficiency.jl new file mode 100644 index 0000000..244eed4 --- /dev/null +++ b/2026-05/examples/07_assimilation_efficiency.jl @@ -0,0 +1,180 @@ +# # [Exercise 07: Assimilation Efficiency] (@id assimilation_efficiency_exercise) +# +# Assimilation efficiency controls how efficiently consumed prey biomass is +# converted into zooplankton biomass. This exercise mirrors the palatability +# matrix exercise: first inspect the derived assimilation matrices, then run each +# configuration in a box model. + +# ## Loading dependencies + +using Agate +using Agate.Introspection: interaction_matrix +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) + +mkpath(joinpath("outputs")) +mkpath(joinpath("figures")) + +nothing #hide + +# ## Construct assimilation-efficiency configurations +# +# Palatability controls which prey are consumed. Assimilation efficiency controls +# how much consumed prey biomass becomes zooplankton biomass. In the default +# model, `assimilation_efficiency` is defined for each zooplankton type and is +# expanded into a consumer-by-prey assimilation matrix. +# +# We compare three cases: +# +# - the default model; +# - a high-assimilation case, where both zooplankton types assimilate consumed +# prey more efficiently; +# - a manual matrix case, where assimilation efficiency depends on both consumer +# and prey identity. + +bgc = default_quickstart_bgc() + +default_assimilation_efficiency = bgc.parameters.assimilation_efficiency +println("Default assimilation efficiency: ", default_assimilation_efficiency) + +high_assimilation_bgc = Agate.Models.NiPiZD.construct(; + parameters = (assimilation_efficiency = (Z1 = 0.8, Z2 = 0.8),), +) + +manual_matrix_bgc = Agate.Models.NiPiZD.construct(; + assimilation_matrix = Float32[0.9 0.25; 0.25 0.9], +) + +bgcs = [bgc, high_assimilation_bgc, manual_matrix_bgc] +case_labels = ["default", "high assimilation", "manual matrix"] + +nothing #hide + +# ## Assimilation matrices +# +# `interaction_matrix(bgc, :assimilation)` returns a labelled interaction table. +# Rows are consumers and columns are prey. The default and high-assimilation +# cases derive this matrix from the zooplankton `assimilation_efficiency` +# parameter. The manual-matrix case supplies the full matrix directly. + +function plot_interaction_matrix!(fig, position, table; title) + matrix = Matrix(table.matrix) + nrows, ncols = size(matrix) + + ax = Axis( + fig[position...]; + title, + xlabel = string(table.column_axis), + ylabel = string(table.row_axis), + xticks = (1:ncols, string.(table.columns)), + yticks = (1:nrows, string.(table.rows)), + yreversed = true, + aspect = DataAspect(), + ) + + hm = heatmap!(ax, 1:ncols, 1:nrows, matrix'; colorrange = (0, 1)) + + for row in 1:nrows, col in 1:ncols + text!( + ax, + col, + row; + text = string(round(matrix[row, col]; digits = 2)), + align = (:center, :center), + fontsize = 11, + ) + end + + return hm +end + +assimilation_figure_path = joinpath("figures", "07_assimilation_matrices.png") +fig_assimilation = Figure(; size = (760, 280), fontsize = 12) + +assimilation_tables = [interaction_matrix(bgc, :assimilation) for bgc in bgcs] +hm = plot_interaction_matrix!(fig_assimilation, (1, 1), assimilation_tables[1]; title = case_labels[1]) +plot_interaction_matrix!(fig_assimilation, (1, 2), assimilation_tables[2]; title = case_labels[2]) +plot_interaction_matrix!(fig_assimilation, (1, 3), assimilation_tables[3]; title = case_labels[3]) +Colorbar(fig_assimilation[1, 4], hm; label = "assimilation efficiency") + +save(assimilation_figure_path, fig_assimilation; px_per_unit = 1) +display(fig_assimilation) +fig_assimilation + +# ## Run box models +# +# We use the same box-model wrapper as in Exercise 02 so the outputs can be read +# by the workshop diagnostics. + +run = run_box_model(bgc; filename = joinpath("outputs", "07_assimilation_default.jld2")) +high_run = run_box_model(high_assimilation_bgc; filename = joinpath("outputs", "07_assimilation_high.jld2")) +manual_run = run_box_model(manual_matrix_bgc; filename = joinpath("outputs", "07_assimilation_manual_matrix.jld2")) + +timeseries = read_box_tracer_timeseries(run.filename, run.tracer_syms) +high_timeseries = read_box_tracer_timeseries(high_run.filename, high_run.tracer_syms) +manual_timeseries = read_box_tracer_timeseries(manual_run.filename, manual_run.tracer_syms) + +times = timeseries.times +data = timeseries.data + +nothing #hide + +# ## Tracer concentration comparison +# +# The concentration diagnostic compares all three assimilation-efficiency configurations. + +comparison_figure_path = joinpath("figures", "07_assimilation_timeseries_comparison.png") +fig_comparison = plot_box_timeseries( + [timeseries, high_timeseries, manual_timeseries]; + labels = case_labels, +) +save(comparison_figure_path, fig_comparison; px_per_unit = 1) +display(fig_comparison) +fig_comparison + +# ## Relative nitrogen contributions: default assimilation + +nitrogen_default_figure_path = joinpath("figures", "07_assimilation_relative_nitrogen_default.png") +fig_nitrogen_default = plot_contributions(timeseries.times, timeseries.data) +save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1) +display(fig_nitrogen_default) +fig_nitrogen_default + +# ## Relative nitrogen contributions: high assimilation + +nitrogen_high_figure_path = joinpath("figures", "07_assimilation_relative_nitrogen_high.png") +fig_nitrogen_high = plot_contributions(high_timeseries.times, high_timeseries.data) +save(nitrogen_high_figure_path, fig_nitrogen_high; px_per_unit = 1) +display(fig_nitrogen_high) +fig_nitrogen_high + +# ## Relative nitrogen contributions: manual matrix + +nitrogen_manual_figure_path = joinpath("figures", "07_assimilation_relative_nitrogen_manual_matrix.png") +fig_nitrogen_manual = plot_contributions(manual_timeseries.times, manual_timeseries.data) +save(nitrogen_manual_figure_path, fig_nitrogen_manual; px_per_unit = 1) +display(fig_nitrogen_manual) +fig_nitrogen_manual + +# ## Community-weighted mean size comparison +# +# Each time series is paired with its matching biogeochemistry object because +# plankton sizes are defined on the model configuration. + +cwm_size_comparison_figure_path = joinpath("figures", "07_assimilation_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size( + [timeseries, high_timeseries, manual_timeseries], + bgcs; + labels = case_labels, +) +save(cwm_size_comparison_figure_path, fig_size_comparison; px_per_unit = 1) +display(fig_size_comparison) +fig_size_comparison diff --git a/2026-05/examples/08_closure_terms.jl b/2026-05/examples/08_closure_terms.jl new file mode 100644 index 0000000..5719634 --- /dev/null +++ b/2026-05/examples/08_closure_terms.jl @@ -0,0 +1,117 @@ +# # [Exercise 08: Closure Terms] (@id closure_terms_exercise) +# +# Mortality terms close the plankton food web by routing unresolved losses back to +# nutrients and detritus. This exercise compares the default model with two +# experiments that remove one closure pathway at a time: +# +# - no linear mortality; +# - no quadratic, density-dependent mortality. + + +# ## Loading dependencies + +using Agate +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) + +mkpath(joinpath("outputs")) +mkpath(joinpath("figures")) + +nothing #hide + +# ## Construct closure-term configurations +# +# The default NiPiZD configuration includes both linear mortality for all +# plankton classes and quadratic mortality for the zooplankton classes. We build +# two sensitivity experiments by setting each mortality vector to zero while +# keeping all other parameters unchanged. + +bgc = default_quickstart_bgc() + +no_linear_bgc = Agate.Models.NiPiZD.construct(; + parameters = (linear_mortality = zero.(bgc.parameters.linear_mortality),), +) + +no_quadratic_bgc = Agate.Models.NiPiZD.construct(; + parameters = (quadratic_mortality = zero.(bgc.parameters.quadratic_mortality),), +) + +bgcs = [bgc, no_linear_bgc, no_quadratic_bgc] +case_labels = ["default", "without linear mortality", "without quadratic mortality"] + +nothing #hide + +# ## Run box models +# +# Each configuration is run with the same box-model setup and the same initial +# conditions helper. The diagnostics compare the resulting tracer trajectories. + +run = run_box_model(bgc; filename = joinpath("outputs", "08_closure_default.jld2")) +no_linear_run = run_box_model(no_linear_bgc; filename = joinpath("outputs", "08_closure_no_linear.jld2")) +no_quadratic_run = run_box_model(no_quadratic_bgc; filename = joinpath("outputs", "08_closure_no_quadratic.jld2")) + +timeseries = read_box_tracer_timeseries(run.filename, run.tracer_syms) +no_linear_timeseries = read_box_tracer_timeseries(no_linear_run.filename, no_linear_run.tracer_syms) +no_quadratic_timeseries = read_box_tracer_timeseries(no_quadratic_run.filename, no_quadratic_run.tracer_syms) + +nothing #hide + +# ## Tracer concentration comparison +# +# This plot compares all tracers across the default and closure-term experiments. + +comparison_figure_path = joinpath("figures", "08_closure_timeseries_comparison.png") +fig_comparison = plot_box_timeseries( + [timeseries, no_linear_timeseries, no_quadratic_timeseries]; + labels = case_labels, +) +save(comparison_figure_path, fig_comparison; px_per_unit = 1) +display(fig_comparison) +fig_comparison + +# ## Community-weighted mean size comparison +# +# The CWM diagnostic compares how removing mortality pathways changes the size +# structure of the plankton community. + +cwm_size_comparison_figure_path = joinpath("figures", "08_closure_cwm_size_comparison.png") +fig_size_comparison = plot_cwm_size( + [timeseries, no_linear_timeseries, no_quadratic_timeseries], + bgcs; + labels = case_labels, +) +save(cwm_size_comparison_figure_path, fig_size_comparison; px_per_unit = 1) +display(fig_size_comparison) +fig_size_comparison + +# ## Relative nitrogen contributions: default closure + +nitrogen_default_figure_path = joinpath("figures", "08_closure_relative_nitrogen_default.png") +fig_nitrogen_default = plot_contributions(timeseries.times, timeseries.data) +save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1) +display(fig_nitrogen_default) +fig_nitrogen_default + +# ## Relative nitrogen contributions: without linear mortality + +nitrogen_no_linear_figure_path = joinpath("figures", "08_closure_relative_nitrogen_no_linear.png") +fig_nitrogen_no_linear = plot_contributions(no_linear_timeseries.times, no_linear_timeseries.data) +save(nitrogen_no_linear_figure_path, fig_nitrogen_no_linear; px_per_unit = 1) +display(fig_nitrogen_no_linear) +fig_nitrogen_no_linear + +# ## Relative nitrogen contributions: without quadratic mortality + +nitrogen_no_quadratic_figure_path = joinpath("figures", "08_closure_relative_nitrogen_no_quadratic.png") +fig_nitrogen_no_quadratic = plot_contributions(no_quadratic_timeseries.times, no_quadratic_timeseries.data) +save(nitrogen_no_quadratic_figure_path, fig_nitrogen_no_quadratic; px_per_unit = 1) +display(fig_nitrogen_no_quadratic) +fig_nitrogen_no_quadratic diff --git a/2026-05/examples/09_diffusivity.jl b/2026-05/examples/09_diffusivity.jl new file mode 100644 index 0000000..3dd5ef4 --- /dev/null +++ b/2026-05/examples/09_diffusivity.jl @@ -0,0 +1,172 @@ +# # [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 diff --git a/2026-05/examples/09_irradiance_box.jl b/2026-05/examples/09_irradiance_box.jl new file mode 100644 index 0000000..c0dde7d --- /dev/null +++ b/2026-05/examples/09_irradiance_box.jl @@ -0,0 +1,136 @@ +# # [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.jl b/2026-05/examples/10_irradiance.jl new file mode 100644 index 0000000..5bd17d9 --- /dev/null +++ b/2026-05/examples/10_irradiance.jl @@ -0,0 +1,211 @@ +# # [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/10_irradiance_column.jl b/2026-05/examples/10_irradiance_column.jl new file mode 100644 index 0000000..525da94 --- /dev/null +++ b/2026-05/examples/10_irradiance_column.jl @@ -0,0 +1,232 @@ +# # [Exercise 10: Irradiance column] (@id irradiance_column_exercise) + +# This exercise focuses on constant irradiance forcing with depth-dependent attenuation 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 + +function discrete_ticks(values) + labels = [v == 0 ? "0" : string(v) for v in values] + return (collect(1:length(values)), labels) +end + +function forcing_heatmap!(fig, row, x, y, values; title, colormap=:viridis) + ax = Axis(fig[row, 1]; xlabel="Time (days)", ylabel="Depth (m)", title) + unique_values = sort(unique(vec(Float64.(values)))) + + if length(unique_values) <= 2 + indices = map(v -> findfirst(isequal(Float64(v)), unique_values), values) + discrete_colormap = length(unique_values) == 1 ? [:gray] : cgrad(colormap, length(unique_values); categorical=true) + hm = CairoMakie.heatmap!( + ax, + x, + y, + indices; + colormap=discrete_colormap, + colorrange=(0.5, length(unique_values) + 0.5), + ) + Colorbar(fig[row, 2], hm; ticks=discrete_ticks(unique_values)) + else + hm = CairoMakie.heatmap!(ax, x, y, values; colormap) + Colorbar(fig[row, 2], hm) + end + + return ax +end + +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) +forcing_heatmap!(fig_forcing, 1, t_range ./ days, z_range, PAR_values; title="irradiance") +forcing_heatmap!(fig_forcing, 2, t_range ./ days, z_range, κₜ_values; title="diffusivity") + +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 = "10_irradiance_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("10_irradiance_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("10_irradiance_column_final_profiles.png", fig_profiles) + +display(fig_profiles) +fig_profiles diff --git a/2026-05/examples/11_diffusivity_stratified.jl b/2026-05/examples/11_diffusivity_stratified.jl new file mode 100644 index 0000000..7688025 --- /dev/null +++ b/2026-05/examples/11_diffusivity_stratified.jl @@ -0,0 +1,197 @@ +# # [Exercise 11: Stratified diffusivity] (@id diffusivity_stratified_exercise) + +# This exercise introduces stratified 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 + +function discrete_ticks(values) + labels = [v == 0 ? "0" : string(v) for v in values] + return (collect(1:length(values)), labels) +end + +function forcing_heatmap!(fig, row, x, y, values; title, colormap=:viridis) + ax = Axis(fig[row, 1]; xlabel="Time (days)", ylabel="Depth (m)", title) + unique_values = sort(unique(vec(Float64.(values)))) + + if length(unique_values) <= 2 + indices = map(v -> findfirst(isequal(Float64(v)), unique_values), values) + discrete_colormap = length(unique_values) == 1 ? [:gray] : cgrad(colormap, length(unique_values); categorical=true) + hm = CairoMakie.heatmap!( + ax, + x, + y, + indices; + colormap=discrete_colormap, + colorrange=(0.5, length(unique_values) + 0.5), + ) + Colorbar(fig[row, 2], hm; ticks=discrete_ticks(unique_values)) + else + hm = CairoMakie.heatmap!(ax, x, y, values; colormap) + Colorbar(fig[row, 2], hm) + end + + return ax +end + +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) +forcing_heatmap!(fig_forcing, 1, t_range ./ days, z_range, PAR_values; title="irradiance") +forcing_heatmap!(fig_forcing, 2, t_range ./ days, z_range, κₜ_values; title="diffusivity") + +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 = "11_diffusivity_stratified.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("11_diffusivity_stratified.png", fig) + +display(fig) +fig # Display the figure diff --git a/2026-05/examples/12_diffusivity_seasonal.jl b/2026-05/examples/12_diffusivity_seasonal.jl new file mode 100644 index 0000000..7a0f406 --- /dev/null +++ b/2026-05/examples/12_diffusivity_seasonal.jl @@ -0,0 +1,172 @@ +# # [Exercise 12: Seasonal diffusivity] (@id diffusivity_seasonal_exercise) + +# This exercise uses a seasonal mixed-layer-depth diffusivity profile 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) + +const 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. Mixed layer depth (MLD) forces the physical mixing (diffusivity), while PAR is held at a fixed surface value with depth-dependent attenuation. +#diffusivity +@inline function diffusivity(x, y, z, t) + H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0) + function fmld1(t) + return H(t, 50days, year) * + (1 / (1 + exp(-(t - 100days) / 5days))) * + (1 / (1 + exp((t - 330days) / 25days))) + end + function MLD(t) + return -( + 10 + + 340 * ( + 1 - fmld1(year - eps(year)) * exp(-mod(t, year) / 25days) - + fmld1(mod(t, year)) + ) + ) + end + return 1e-2 * (1 + tanh((z - MLD(t)) / 10)) / 2 + 1e-4 +end + +#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(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, κ=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 = "12_diffusivity_seasonal.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("12_diffusivity_seasonal.png", fig) + +display(fig) +fig # Display the figure diff --git "a/2026-05/papers/Ki\303\270rboe2024.pdf" "b/2026-05/papers/Ki\303\270rboe2024.pdf" new file mode 100644 index 0000000..c1930df Binary files /dev/null and "b/2026-05/papers/Ki\303\270rboe2024.pdf" differ diff --git a/2026-05/papers/bagwell2024.pdf b/2026-05/papers/bagwell2024.pdf new file mode 100644 index 0000000..154cd6f Binary files /dev/null and b/2026-05/papers/bagwell2024.pdf differ diff --git a/2026-05/papers/follows2011.pdf b/2026-05/papers/follows2011.pdf new file mode 100644 index 0000000..cb09a96 Binary files /dev/null and b/2026-05/papers/follows2011.pdf differ diff --git a/2026-05/papers/gentleman2002.pdf b/2026-05/papers/gentleman2002.pdf new file mode 100644 index 0000000..fd0c692 Binary files /dev/null and b/2026-05/papers/gentleman2002.pdf differ diff --git a/2026-05/papers/hansen1994.pdf b/2026-05/papers/hansen1994.pdf new file mode 100644 index 0000000..39b23ec Binary files /dev/null and b/2026-05/papers/hansen1994.pdf differ diff --git a/2026-05/setup/Dockerfile b/2026-05/setup/Dockerfile index 4ae2824..7947a33 100644 --- a/2026-05/setup/Dockerfile +++ b/2026-05/setup/Dockerfile @@ -236,4 +236,4 @@ RUN chmod +x /usr/local/bin/start-code-server WORKDIR /workspace EXPOSE 8080 -CMD ["/usr/local/bin/start-code-server"] \ No newline at end of file +CMD ["/usr/local/bin/start-code-server"] diff --git a/2026-05/setup/Project.toml b/2026-05/setup/Project.toml index 1fbb337..9d8acc0 100644 --- a/2026-05/setup/Project.toml +++ b/2026-05/setup/Project.toml @@ -11,7 +11,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" [compat] -Agate = "0.5.2" +Agate = "0.5.3" CairoMakie = "0.12, 0.13, 0.14, 0.15" CSV = "0.10" DataFrames = "1" diff --git a/2026-05/slides/01_intro.qmd b/2026-05/slides/01_intro.qmd new file mode 100644 index 0000000..f7191c0 --- /dev/null +++ b/2026-05/slides/01_intro.qmd @@ -0,0 +1,254 @@ +--- +title: "Agate.jl Workshop: Motivation and Orientation" +author: "" +output-file: 01_intro.html +--- + +## Welcome {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} + +| Time | Session | +|---|---| +| 09:30–10:00 | *Docker setup (optional, this can be done before the workshop)* | +| 10:00–10:30 | Welcome and workshop overview | +| 10:30–11:10 | Session 1: Size, allometry, and plankton traits | +| 11:10–11:25 | Break | +| 11:25–12:05 | Session 2: Predation and trophic structure | +| 12:05–12:45 | Session 3: Physical forcing: light and diffusivity | +| 12:45–13:45 | Lunch | +| 13:45–14:35 | Group exercise | +| 14:35–15:10 | Own-work block 1 | +| 15:10–15:40 | Break | +| 15:40–16:30 | Own-work block 2 | +| 16:30–16:55 | Show-and-tell (optional) and feedback | +| 16:55–17:00 | Wrap-up | +| 17:00 | Pub | + +::: + +::: {.slide-figure} +![](images/01_intro/slide-00.svg) +

Map of 8-10 Berkley Square (from [mazemap](https://use.mazemap.com/#v=1&config=UoBCampuses&campusid=843¢er=-2.606670,51.454875&zoom=20&zlevel=2)).

+::: + +## Current state of the art {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} + +All major Earth System Models (ESMs) use FORTRAN. +- First appeared in 1957; 69 years ago +- Highly optimized and battle tested +- Large user base / existing expertise +- Physical system well calibrated (circulation, boundary conditions) + +Downsides: + +- Ocean biogeochemistry often secondary +- Non-composable -> iron cycling for PISCES can't be used in DARWIN +- Hard to Automatically Differentiate -> makes Bayesian inference and machine learning expensive! +- Difficult to utilize GPUs (but increasingly feasible) + +::: + +::: {.slide-figure} +![](images/01_intro/slide-01.svg) +

FORTRAN on punched card from Wikimedia Commons, TIOBE May 2026 rankings from tiobe.com.

+::: + +## Julia programming language - a new frontier ? {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} + +Benefits: + +- Similar performance to FORTRAN but easier to program (like Python, R, MATLAB) +- Designed to be composable (at least in theory) +- Strong Automatic Differentiation infrastructure +- GPUs well supported +- Existing physical oceanography model (Oceananigans.jl) and an active Earth System Modelling community (NumericalEarth.jl) + +Downsides: + +- Effectively starting from scratch (ESMs are complex(!)) +- Small user base / expertise +- Plotting etc. is still behind Python/R +- Senior PI's are often hesitant to learn a new language +- FORTRAN has stood the test of time, will Julia? + +::: + +::: {.slide-figure} +![](images/01_intro/slide-02.svg) +

Figure from Silverstri et al., 2023 (ArXiv preprint)

+::: + +## From monolithic models to composable components {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Current Earth System and Ocean Models are modular (e.g. NEMO and PISCES/ MITgcm and DARWIN) but not composable with each other (no MITgcm-PISCES) +- NumericalEarth.jl aims to provide an composable alternative (at least in theory) +- For example: atmospheric physics can be resolved by Breeze.jl, SpeedyWeather.jl or data products (JRA5/ERA5) +- Only one ocean physics model (Oceananigans.jl) - so far (!) +- Automatic Differentiation a key objective + +::: + +::: {.slide-figure} +![](images/01_intro/slide-03.svg) +

NumericalEarth.jl provides a coupling frameworks to make Julia composable earth system models.

+::: + +## Oceananigans.jl and OceanBioME.jl {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +Oceananigans.jl: + +- initial commit Oct 2018 +- published Dec 2019 +- ~ 90 contributors + +OceanBioME.jl: + +- initial commit Jul 2022 +- published Oct 2023 +- ~ 14 contributors + +::: + +::: {.slide-figure} +![](images/01_intro/slide-04.svg) +

Oceans in Julia: Oceananigans.jl provides physics and infrastructure such as solvers, OceanBioME.jl provides biogeochemistry

+::: + +## Where does Agate.jl fit in? {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- **A**quatic **G**CM **A**gnostic **T**rait-based **E**cosystems +- Couples with OceanBioME.jl +- Used to implement high-complexity trait-based ecosystems of arbitrary complexity +- First commit Nov 2018, manuscript in prep... + + +::: + +::: {.slide-figure} +![](images/01_intro/slide-05.svg) +

Agate.jl couples with OceanBioME and is solved by Oceananigans.jl.

+::: + +## Intermediate vs High complexity ecosystems {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} + +- Low complexity ecosystems: NPZD models +- Intermediate complexity: multiple P and Z functional types (e.g. PISCES, MARBL) +- High complexity: size structured plus functional types (e.g. DARWIN, ecoGENIE) + + +::: + +::: {.slide-figure} +![](images/01_intro/slide-06.svg) +

Simplified representation of model complexity. Some models favour biogeochemical complexity (e.g. PISCES) while others favour ecosystems (DARWIN).

+::: + +## Trait-based approaches {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- ~ 110,000 eukaryotic OTU's (de Vargas et al., 2015) +- modelling even 1% of this is too expensive +- Large majority uncultured and can't be easily parameterized +- Trait based frameworks provide first-principles framework to capture this diversity + +::: + +::: {.slide-figure} +![](images/01_intro/slide-07.svg) +

Overview of key traits from Litchman and Klaushmeier, 2008 (Ann. Rev. Ecol. Evol. Syst. ).

+::: + +## Trait trade-offs {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Complex emergent ecossytems can be created by seeding the system with many types +- Competitive types will survive in relevant niches +- idea requires trait trade-offs, if one type is too advantageous it will dominate (Darwinian demon) +- trade-off: "every benefit comes at a cost + +::: + +::: {.slide-figure} +![](images/01_intro/slide-08.svg) +

Trait trade-offs. Only some combinations are energetically feasible, and a subset is selected by the environment

+::: + +## Allometry {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Allometry (size) is a funamental trait that applies to all organisms on ERA5 +- For plankton it determines, growth rates and nutrient requirements +- Prey selection is also largely size dependant + +::: + +::: {.slide-figure} +![](images/01_intro/slide-09.svg) +

Division rate vs size from Follows and Dutkiewicz, 2011. Prey:predator size ratios from Hansen 1994

+::: + +## Global emergent ecosystems {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Trait-based models successfully used in global context +- MITgcm-DARWIN (modern and climate change simulations of carbon cycle and diversity) +- ecoGENIE (paleo and future climate simulations of carbon cycle and diversity) + +::: + +::: {.slide-figure} +![](images/01_intro/slide-10.svg) +

Global applications of trait-based models. MITgcm-DARWIN figure from Dutkiewicz and EcoGENIE from Bagwell

+::: + +## Setup for today: box models {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Global trait-based models are too expensive to run +- We will be using "0D" box models and 1D water columns instead +- We will use these to explore the impacts of traits and physical processes (mixing and light) + +::: + +::: {.slide-figure} +![](images/01_intro/slide-11.svg) +

Overview of box-model reproduced from dar_one docs.

+::: + +## Schedule {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} + +| Time | Session | +|---|---| +| 09:30–10:00 | *Docker setup (optional, this can be done before the workshop)* | +| 10:00–10:30 | Welcome and workshop overview | +| 10:30–11:10 | Session 1: Size, allometry, and plankton traits | +| 11:10–11:25 | Break | +| 11:25–12:05 | Session 2: Predation and trophic structure | +| 12:05–12:45 | Session 3: Physical forcing: light and diffusivity | +| 12:45–13:45 | Lunch | +| 13:45–14:35 | Group exercise | +| 14:35–15:10 | Own-work block 1 | +| 15:10–15:40 | Break | +| 15:40–16:30 | Own-work block 2 | +| 16:30–16:55 | Show-and-tell (optional) and feedback | +| 16:55–17:00 | Wrap-up | +| 17:00 | Pub | + +::: + +::: {.slide-figure} +![](images/01_intro/slide-00.svg) +

Map of 8-10 Berkley Square (from [mazemap](https://use.mazemap.com/#v=1&config=UoBCampuses&campusid=843¢er=-2.606670,51.454875&zoom=20&zlevel=2)).

+::: \ No newline at end of file diff --git a/2026-05/slides/02_module_lecture.qmd b/2026-05/slides/02_module_lecture.qmd new file mode 100644 index 0000000..5c5b952 --- /dev/null +++ b/2026-05/slides/02_module_lecture.qmd @@ -0,0 +1,99 @@ +--- +title: "Module Lecture 1" +author: "" +--- + +## Size as fundamental trait {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Universal trait (bacteria to whales) +- Huge size range of plankton in ocean +::: + +::: {.slide-figure} +![](images/02_allometry/slide-01.svg) +

Modified from: Holland & McQuatters-Gollop (2024). mNCEA policy brief – The many scales of pelagic habitats. Defra mNCEA Programme – Pelagic Natural Capital. University of Plymouth.

+::: + +## Surface Area to Volume {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Allometric scaling a function of SA:V ratio +::: + +::: {.slide-figure} +![](images/02_allometry/slide-02.svg) +

Replace with the workflow or diagnostic used in this module.

+::: + + +## Maximum growth rate {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Allometric scaling a function of SA:V ratio +::: + +::: {.slide-figure} +![](images/02_allometry/slide-03.svg) +

Replace with the workflow or diagnostic used in this module.

+::: + + +## Maximum grazing rate {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- growth vs prey +::: + +::: {.slide-figure} +![](images/02_allometry/slide-04.svg) +

Replace with the workflow or diagnostic used in this module.

+::: + +## Nutrient half saturation {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- nutrients vs growth +::: + +::: {.slide-figure} +![](images/02_allometry/slide-05.svg) +

Replace with the workflow or diagnostic used in this module.

+::: + + +## Allometric scaling {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- aV^b +::: + +::: {.slide-figure} +![](images/02_allometry/slide-06.svg) +

Replace with the workflow or diagnostic used in this module.

+::: + + +## Number of groups {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- How many are needed to capture diversity? +::: + +::: {.slide-figure} +![](images/02_allometry/slide-07.svg) +

Replace with the workflow or diagnostic used in this module.

+::: + +## Size ranges {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Are there physiological limitations? +- Biogeochemical ones? +- Physics (e.g. sinking)? +::: + +::: {.slide-figure} +![](images/02_allometry/slide-08.svg) +

Replace with the workflow or diagnostic used in this module.

+::: \ No newline at end of file diff --git a/2026-05/slides/03_module_lecture.qmd b/2026-05/slides/03_module_lecture.qmd new file mode 100644 index 0000000..e95e149 --- /dev/null +++ b/2026-05/slides/03_module_lecture.qmd @@ -0,0 +1,76 @@ +--- +title: "Module Lecture 2" +author: "" +--- + +## Predation and trophic structure {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- +::: + +::: {.slide-figure} +![](images/03_predation/slide-01.svg) +

Image from Dutkiewicz et al., 2024 (L&O).

+::: + +## Palatability {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- +::: + +::: {.slide-figure} +![](images/03_predation/slide-02.svg) +

Replace with the diagnostic or comparison for this module.

+::: + + +## Size optima {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- +::: + +::: {.slide-figure} +![](images/03_predation/slide-03.svg) +

Hansen 1994, and from Garcia-Oliva and Wirtz 2022.

+::: + + + +## Assimilation efficiency {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Also known as Gross Growth Efficiency (GGE) +::: + +::: {.slide-figure} +![](images/03_predation/slide-04.svg) +

Data replotted from Straile 1997.

+::: + + +## Unresolved higher predation {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- +::: + +::: {.slide-figure} +![](images/03_predation/slide-05.svg) +

Image from Follows and Dutkiewicz, 2011 (Annu. Rev. Mar. Sci).

+::: + + +## Trade-offs {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- +::: + +::: {.slide-figure} +![](images/03_predation/slide-06.svg) +

Replace with the diagnostic or comparison for this module.

+::: + diff --git a/2026-05/slides/04_module_lecture.qmd b/2026-05/slides/04_module_lecture.qmd new file mode 100644 index 0000000..effee2a --- /dev/null +++ b/2026-05/slides/04_module_lecture.qmd @@ -0,0 +1,44 @@ +--- +title: "Physical forcings: irradiance and mixing" +author: "" +--- + +## Photosynthesis {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Replace this text with the concept lecture for the third hands-on module. +- State the modelling question clearly. +- Identify the assumption or parameterization participants will modify. +::: + +::: {.slide-figure} +![](images/03_physics/slide-01.svg) +

Replace with a module-specific figure caption.

+::: + + +## Light attenuation {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Replace this text with the concept lecture for the third hands-on module. +- State the modelling question clearly. +- Identify the assumption or parameterization participants will modify. +::: + +::: {.slide-figure} +![](images/03_physics/slide-01.svg) +

Replace with a module-specific figure caption.

+::: + +## Mixing {background-image="images/background-placeholder.svg" background-size="cover"} + +::: {.slide-copy} +- Replace this text with the concept lecture for the third hands-on module. +- State the modelling question clearly. +- Identify the assumption or parameterization participants will modify. +::: + +::: {.slide-figure} +![](images/03_physics/slide-01.svg) +

Lindemann and John, 2014 (Frontiers in Marine Science).

+::: \ No newline at end of file diff --git a/2026-05/slides/_quarto.yml b/2026-05/slides/_quarto.yml new file mode 100644 index 0000000..c60c3ff --- /dev/null +++ b/2026-05/slides/_quarto.yml @@ -0,0 +1,23 @@ +project: + type: website + output-dir: _site + render: + - 01_intro.qmd + - 02_module_lecture.qmd + - 03_module_lecture.qmd + - 04_module_lecture.qmd + +format: + revealjs: + width: 1920 + height: 1080 + theme: simple + slide-number: true + chalkboard: false + preview-links: auto + css: styles.css + title-slide-attributes: + data-background-image: images/background-placeholder.svg + data-background-size: cover +execute: + echo: false diff --git a/2026-05/slides/images/01_intro/slide-00.svg b/2026-05/slides/images/01_intro/slide-00.svg new file mode 100644 index 0000000..2dcb11f --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-00.svg @@ -0,0 +1,56 @@ + + diff --git a/2026-05/slides/images/01_intro/slide-01.svg b/2026-05/slides/images/01_intro/slide-01.svg new file mode 100644 index 0000000..b81633e --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-01.svg @@ -0,0 +1,93 @@ + +FORTRAN code on a punched cardFORTRAN 11th in popularity (TIOBE) diff --git a/2026-05/slides/images/01_intro/slide-02.svg b/2026-05/slides/images/01_intro/slide-02.svg new file mode 100644 index 0000000..8935864 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-02.svg @@ -0,0 +1,47 @@ + + diff --git a/2026-05/slides/images/01_intro/slide-03.svg b/2026-05/slides/images/01_intro/slide-03.svg new file mode 100644 index 0000000..ae86c69 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-03.svg @@ -0,0 +1,246 @@ + +Atmospheric physics(JRA5, Breeze.jl / SpeedyWeather.jl etc.)Coupling(NumericalEarth.jl)Sea ice physics(ClimaSeaIce.jl)Terrestrial ecosystems(Terrarium.jl)Ocean physics(Oceananigans.jl)Ocean biogeochemistry(OceanBioME.jl) diff --git a/2026-05/slides/images/01_intro/slide-04.svg b/2026-05/slides/images/01_intro/slide-04.svg new file mode 100644 index 0000000..41442c7 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-04.svg @@ -0,0 +1,256 @@ + +Air-sea gas exchangeOcean Physics + solver(Oceananigans.jl)Ocean Biogeochemistry(OceanBioME.jl)Individual-basedparticles (e.g. kelp)Biogeochemicaltracers (e.g. PISCES)Carbonate chemistrySedimentmodelsLight Attenuation diff --git a/2026-05/slides/images/01_intro/slide-05.svg b/2026-05/slides/images/01_intro/slide-05.svg new file mode 100644 index 0000000..66353ce --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-05.svg @@ -0,0 +1,263 @@ + +Air-sea gas exchangeOcean Physics + solver(Oceananigans.jl)Ocean Biogeochemistry(OceanBioME.jl)Individual-basedparticles (e.g. kelp)High-complexityEcosystems(Agate.jl)Carbonate chemistrySedimentmodelsLight Attenuation diff --git a/2026-05/slides/images/01_intro/slide-06.svg b/2026-05/slides/images/01_intro/slide-06.svg new file mode 100644 index 0000000..7f412cb --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-06.svg @@ -0,0 +1,470 @@ + +Phytoplankton(n_size x n_type)Zooplankton(n_size)High-complexityEcosystemDetritusNutrientsPhytoplankton(3 types)Zooplankton(2 types)DetritusNutrientsMedium-complexityBiogeochemistryHigh-complexityBiogeochemistryMedium-complexityEcosystem diff --git a/2026-05/slides/images/01_intro/slide-07.svg b/2026-05/slides/images/01_intro/slide-07.svg new file mode 100644 index 0000000..dc10fa1 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-07.svg @@ -0,0 +1,40 @@ + + diff --git a/2026-05/slides/images/01_intro/slide-08.svg b/2026-05/slides/images/01_intro/slide-08.svg new file mode 100644 index 0000000..04a2989 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-08.svg @@ -0,0 +1,53 @@ + + diff --git a/2026-05/slides/images/01_intro/slide-09.svg b/2026-05/slides/images/01_intro/slide-09.svg new file mode 100644 index 0000000..b764470 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-09.svg @@ -0,0 +1,47 @@ + + diff --git a/2026-05/slides/images/01_intro/slide-10.svg b/2026-05/slides/images/01_intro/slide-10.svg new file mode 100644 index 0000000..cc93ad4 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-10.svg @@ -0,0 +1,79 @@ + +MITgcm-DARWIN modelEcoGEnIE model diff --git a/2026-05/slides/images/01_intro/slide-11.svg b/2026-05/slides/images/01_intro/slide-11.svg new file mode 100644 index 0000000..d2a97c2 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-11.svg @@ -0,0 +1,40 @@ + + diff --git a/2026-05/slides/images/01_intro/slide-12.svg b/2026-05/slides/images/01_intro/slide-12.svg new file mode 100644 index 0000000..f7f36a1 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-12.svg @@ -0,0 +1,6 @@ + + + + Slide 12 image placeholder + Replace images/slide-12.svg + diff --git a/2026-05/slides/images/01_intro/slide-13.svg b/2026-05/slides/images/01_intro/slide-13.svg new file mode 100644 index 0000000..e799310 --- /dev/null +++ b/2026-05/slides/images/01_intro/slide-13.svg @@ -0,0 +1,6 @@ + + + + Slide 13 image placeholder + Replace images/slide-13.svg + diff --git a/2026-05/slides/images/02_allometry/slide-01.svg b/2026-05/slides/images/02_allometry/slide-01.svg new file mode 100644 index 0000000..7c82a76 --- /dev/null +++ b/2026-05/slides/images/02_allometry/slide-01.svg @@ -0,0 +1,40 @@ + + diff --git a/2026-05/slides/images/02_allometry/slide-02.svg b/2026-05/slides/images/02_allometry/slide-02.svg new file mode 100644 index 0000000..04b9858 --- /dev/null +++ b/2026-05/slides/images/02_allometry/slide-02.svg @@ -0,0 +1,40 @@ + + diff --git a/2026-05/slides/images/02_allometry/slide-03.svg b/2026-05/slides/images/02_allometry/slide-03.svg new file mode 100644 index 0000000..465cfe5 --- /dev/null +++ b/2026-05/slides/images/02_allometry/slide-03.svg @@ -0,0 +1,40 @@ + + diff --git a/2026-05/slides/images/02_allometry/slide-05.svg b/2026-05/slides/images/02_allometry/slide-05.svg new file mode 100644 index 0000000..e9dd73d --- /dev/null +++ b/2026-05/slides/images/02_allometry/slide-05.svg @@ -0,0 +1,40 @@ + + diff --git a/2026-05/slides/images/03_predation/slide-01.svg b/2026-05/slides/images/03_predation/slide-01.svg new file mode 100644 index 0000000..9fd3a5f --- /dev/null +++ b/2026-05/slides/images/03_predation/slide-01.svg @@ -0,0 +1,47 @@ + + + + diff --git a/2026-05/slides/images/03_predation/slide-02.svg b/2026-05/slides/images/03_predation/slide-02.svg new file mode 100644 index 0000000..09d5b2b --- /dev/null +++ b/2026-05/slides/images/03_predation/slide-02.svg @@ -0,0 +1,244 @@ + + + +PalatabilityHigh palatabilityP1P2ZLow palatabilitySilicificationCalcificationToxinsNoneDefense diff --git a/2026-05/slides/images/03_predation/slide-03.svg b/2026-05/slides/images/03_predation/slide-03.svg new file mode 100644 index 0000000..c7f8f0b --- /dev/null +++ b/2026-05/slides/images/03_predation/slide-03.svg @@ -0,0 +1,53 @@ + + + + diff --git a/2026-05/slides/images/03_predation/slide-04.svg b/2026-05/slides/images/03_predation/slide-04.svg new file mode 100644 index 0000000..ce02d84 --- /dev/null +++ b/2026-05/slides/images/03_predation/slide-04.svg @@ -0,0 +1,168 @@ + + + +grazing * (1-efficiency)D/NPZgrazing * efficiencyAssimilation Efficiency(aka Gross Growth Efficiency) diff --git a/2026-05/slides/images/03_predation/slide-05.svg b/2026-05/slides/images/03_predation/slide-05.svg new file mode 100644 index 0000000..2711713 --- /dev/null +++ b/2026-05/slides/images/03_predation/slide-05.svg @@ -0,0 +1,47 @@ + + + + diff --git a/2026-05/slides/images/03_predation/slide-06.svg b/2026-05/slides/images/03_predation/slide-06.svg new file mode 100644 index 0000000..652511a --- /dev/null +++ b/2026-05/slides/images/03_predation/slide-06.svg @@ -0,0 +1,46 @@ + + + + diff --git a/2026-05/slides/images/background-placeholder.svg b/2026-05/slides/images/background-placeholder.svg new file mode 100644 index 0000000..4e1e09f --- /dev/null +++ b/2026-05/slides/images/background-placeholder.svg @@ -0,0 +1,58 @@ + +background placeholder diff --git a/2026-05/slides/styles.css b/2026-05/slides/styles.css new file mode 100644 index 0000000..3366509 --- /dev/null +++ b/2026-05/slides/styles.css @@ -0,0 +1,93 @@ +/* Basic 1920x1080 Reveal.js layout for the Agate.jl opening lecture. */ + +:root { + --content-top: 155px; + --content-left: 40px; + --content-gap: 100px; + --copy-width: 700px; + --figure-max-width: 950px; + --figure-max-height: 800px; + --panel-radius: 18px; +} + +.reveal .slides section { + text-align: left; +} + +.reveal h2 { + font-size: 1.65em; + line-height: 1.08; + margin-top: 0; + margin-bottom: 0.45em; +} + +.reveal .slide-copy { + position: absolute; + left: var(--content-left); + top: var(--content-top); + width: var(--copy-width); + font-size: 0.78em; + line-height: 1.28; + padding: 0.7em 0.9em; + background: rgba(255, 255, 255, 0.82); + border-radius: var(--panel-radius); +} + +.reveal .slide-copy table { + font-size: 0.72em; + width: 100%; +} + +.reveal .slide-copy th, +.reveal .slide-copy td { + padding: 0.25em 0.35em; + vertical-align: top; +} + +.reveal .slide-figure { + position: absolute; + left: calc(var(--content-left) + var(--copy-width) + var(--content-gap)); + top: var(--content-top); + width: var(--figure-width); + text-align: center; +} + +.reveal .slide-figure img { + display: block; + max-width: var(--figure-max-width); + max-height: var(--figure-max-height); + width: auto; + height: auto; + margin: 0 auto; + object-fit: contain; + background: none; + border-radius: 0; + box-shadow: none; + padding: 0; +} + +.reveal .slide-figure .caption { + display: inline-block; + margin: 0.5em auto 0; + max-width: 95%; + padding: 0.35em 0.55em; + font-size: 0.55em; + line-height: 1.25; + color: #111; + text-align: center; + background: rgba(255, 255, 255, 0.88); + border-radius: 10px; +} + +.reveal .title-slide h1, +.reveal .title-slide .title { + max-width: 1180px; + font-size: 1.9em; + line-height: 1.05; + color: #000; +} + +.reveal .title-slide .quarto-title-authors, +.reveal .title-slide .author { + color: #000; +} diff --git a/2026-05/src/AgateWorkshop.jl b/2026-05/src/AgateWorkshop.jl new file mode 100644 index 0000000..3eadbad --- /dev/null +++ b/2026-05/src/AgateWorkshop.jl @@ -0,0 +1,451 @@ +using Agate +using Agate.Introspection: plankton_diameters, plankton_groups, plankton_tracers, tracer_names +using Agate.Library.Light +using CairoMakie +using OceanBioME: Biogeochemistry, BoxModel, BoxModelGrid +using Oceananigans +using Oceananigans.Units: day, minute + +const DEFAULT_INITIAL_NUTRIENT = 7.0 +const DEFAULT_INITIAL_DETRITUS = 0.01 +const DEFAULT_INITIAL_PLANKTON_BIOMASS = 0.13 + +function default_initial_conditions( + bgc; + nutrient=DEFAULT_INITIAL_NUTRIENT, + detritus=DEFAULT_INITIAL_DETRITUS, + total_plankton_biomass=DEFAULT_INITIAL_PLANKTON_BIOMASS, +) + plankton = Symbol.(collect(plankton_tracers(bgc))) + !isempty(plankton) || error("Cannot set initial plankton biomass: no plankton tracers found.") + + plankton_biomass = total_plankton_biomass / length(plankton) + names = (:N, plankton..., :D) + values = (nutrient, fill(plankton_biomass, length(plankton))..., detritus) + + return NamedTuple{names}(values) +end + +const WORKSHOP_COLORS = ( + living = "#1F77B4", + nonliving = "#D62728", + phytoplankton = "#6BAED6", + zooplankton = "#08519C", +) + +default_quickstart_bgc() = Agate.Models.NiPiZD.construct() + +default_quickstart_initial_conditions(bgc=default_quickstart_bgc()) = default_initial_conditions(bgc) + +function build_box_model( + bgc; + light_attenuation=FunctionFieldPAR(; grid=BoxModelGrid()), + initial_conditions=nothing, +) + initial_conditions === nothing && (initial_conditions = default_initial_conditions(bgc)) + + bgc_model = Biogeochemistry(bgc; light_attenuation) + model = BoxModel(; biogeochemistry=bgc_model) + set!(model; initial_conditions...) + return model +end + +function run_box_model( + bgc; + filename=joinpath("outputs", "quick_start.jld2"), + initial_conditions=nothing, + Δt=240minute, + stop_time=1095day, + output_interval=1day, + light_attenuation=FunctionFieldPAR(; grid=BoxModelGrid()), +) + mkpath(dirname(filename)) + + model = build_box_model(bgc; light_attenuation, initial_conditions) + simulation = Simulation(model; Δt, stop_time) + + simulation.output_writers[:fields] = JLD2Writer( + model, + model.fields; + filename, + schedule=TimeInterval(output_interval), + overwrite_existing=true, + ) + + run!(simulation) + + return ( + model=model, + simulation=simulation, + filename=filename, + tracer_syms=tracer_names(bgc), + ) +end + +function _read_field_values(field_time_series) + return [first(interior(field_time_series[n])) for n in eachindex(field_time_series.times)] +end + +function read_box_tracer_timeseries(filename, tracer_syms) + tracer_syms = collect(Symbol.(tracer_syms)) + data = Dict{Symbol, Vector{Float64}}() + times = nothing + + for tracer in tracer_syms + field_time_series = FieldTimeSeries(filename, string(tracer)) + times === nothing && (times = collect(field_time_series.times) ./ day) + data[tracer] = Float64.(_read_field_values(field_time_series)) + end + + return (times=times, data=data) +end + +_times(ts) = hasproperty(ts, :times) ? getproperty(ts, :times) : ts[1] +_data(ts) = hasproperty(ts, :data) ? getproperty(ts, :data) : ts[2] + +function _keys_from_data(data) + return sort!(collect(keys(data)); by=string) +end + +function _prefixed_keys(data, prefix) + return filter(key -> startswith(String(key), prefix), _keys_from_data(data)) +end + +function _sum_keys(data, keys) + isempty(keys) && return zeros(length(first(values(data)))) + total = zero(data[first(keys)]) + for key in keys + total = total .+ data[key] + end + return total +end + +function _plankton_keys(data) + phyto = _prefixed_keys(data, "P") + zoo = _prefixed_keys(data, "Z") + return phyto, zoo, vcat(phyto, zoo) +end + +function _plankton_diameters_by_tracer(bgc) + return Dict(Symbol.(plankton_tracers(bgc)) .=> plankton_diameters(bgc)) +end + +function _keys_with_diameters(keys, diameters) + return filter(key -> haskey(diameters, key), keys) +end + +function _cwm_or_nan(data, diameters, keys) + keys = _keys_with_diameters(keys, diameters) + isempty(keys) && return fill(NaN, length(first(values(data)))) + + total = zero(data[first(keys)]) + weighted = zero(data[first(keys)]) + + for key in keys + values = data[key] + total = total .+ values + weighted = weighted .+ diameters[key] .* values + end + + return weighted ./ ifelse.(total .> 0, total, NaN) +end + +function _save_if_requested(fig, figure_path) + figure_path !== nothing && save(figure_path, fig; px_per_unit=1) + return fig +end + +function _case_palette(n) + colors = ["#0072B2", "#E69F00", "#009E73", "#CC79A7", "#56B4E9", "#D55E00", "#F0E442"] + return [colors[mod1(i, length(colors))] for i in 1:n] +end + +function _shared_line_legend!(fig, row, columns, colors, labels) + elements = [LineElement(; color=color, linewidth=2) for color in colors] + Legend( + fig[row, columns], + elements, + labels; + orientation=:horizontal, + tellwidth=false, + tellheight=true, + framevisible=false, + ) +end + +function _parameter_tracers_and_values(bgc, parameter_name::Symbol) + groups = plankton_groups(bgc) + plankton = collect(plankton_tracers(bgc)) + values = collect(getproperty(bgc.parameters, parameter_name)) + + if length(values) == length(groups.P) + return collect(groups.P), values, "Phytoplankton" + elseif length(values) == length(groups.Z) + return collect(groups.Z), values, "Zooplankton" + elseif length(values) == length(plankton) + return plankton, values, "Plankton" + else + error("Cannot align parameter $parameter_name with plankton tracers.") + end +end + +function plot_plankton_parameter_bars( + bgcs, + parameter_name::Symbol; + labels=nothing, + ylabel=string(parameter_name), + title=string(parameter_name), + figure_path=nothing, +) + bgcs = bgcs isa Tuple || bgcs isa AbstractVector ? collect(bgcs) : [bgcs] + !isempty(bgcs) || error("At least one biogeochemistry object is required.") + + labels === nothing && (labels = ["case $i" for i in eachindex(bgcs)]) + length(labels) == length(bgcs) || error("labels must have one entry per biogeochemistry object.") + + tracers, first_values, group_label = _parameter_tracers_and_values(first(bgcs), parameter_name) + value_sets = [first_values] + + for bgc in bgcs[2:end] + case_tracers, values, case_group_label = _parameter_tracers_and_values(bgc, parameter_name) + case_tracers == tracers || error("All cases must use the same tracer names for $parameter_name.") + case_group_label == group_label || error("All cases must align $parameter_name with the same plankton group.") + push!(value_sets, values) + end + + n_tracers = length(tracers) + n_cases = length(bgcs) + centers = collect(1:n_tracers) + width = min(0.8 / n_cases, 0.28) + offsets = ((1:n_cases) .- (n_cases + 1) / 2) .* width + + fig = Figure(; size=(max(420, 80 * n_tracers), 260), fontsize=12) + ax = Axis( + fig[1, 1]; + xlabel=group_label, + ylabel, + title, + xticks=(centers, string.(tracers)), + ) + + for (case_idx, values) in enumerate(value_sets) + barplot!(ax, centers .+ offsets[case_idx], values; width, label=labels[case_idx]) + end + + n_cases > 1 && axislegend(ax; position=:rt) + + return _save_if_requested(fig, figure_path) +end + +function _as_timeseries_vector(timeseries) + if hasproperty(timeseries, :times) && hasproperty(timeseries, :data) + return [timeseries] + elseif timeseries isa AbstractVector || timeseries isa Tuple + series = collect(timeseries) + !isempty(series) || error("At least one box-model time series is required.") + all(ts -> hasproperty(ts, :times) && hasproperty(ts, :data), series) || + error("Each entry must have times and data fields.") + return series + else + error("Expected a box-model time series or a collection of box-model time series.") + end +end + +function _box_timeseries_variables(series, variables) + if variables !== nothing + return collect(Symbol.(variables)) + end + + keys_in_all = Set(keys(_data(first(series)))) + for ts in series[2:end] + intersect!(keys_in_all, keys(_data(ts))) + end + + return sort!(collect(keys_in_all); by=string) +end + +function plot_box_timeseries( + timeseries; + labels=nothing, + variables=nothing, + ylabels=nothing, + figure_path=nothing, +) + series = _as_timeseries_vector(timeseries) + variables = _box_timeseries_variables(series, variables) + + if labels === nothing + labels = length(series) == 1 ? ["box model"] : ["series $i" for i in eachindex(series)] + end + + length(labels) == length(series) || error("labels must have one entry per box-model time series.") + + if ylabels === nothing + ylabels = Dict(variable => "Concentration (mmol N m⁻³)" for variable in variables) + else + ylabels = Dict(Symbol(key) => value for (key, value) in pairs(ylabels)) + end + + n_variables = length(variables) + n_columns = min(2, n_variables) + n_rows = cld(n_variables, n_columns) + has_comparison = length(series) > 1 + legend_rows = has_comparison ? 1 : 0 + colors = _case_palette(length(series)) + fig = Figure(; size=(360 * n_columns, 180 * n_rows + 44 * legend_rows), fontsize=12) + + if has_comparison + _shared_line_legend!(fig, 1, 1:n_columns, colors, labels) + end + + for (idx, variable) in enumerate(variables) + row = cld(idx, n_columns) + legend_rows + col = mod1(idx, n_columns) + ax = Axis( + fig[row, col]; + title=string(variable), + xlabel="Days", + ylabel=get(ylabels, variable, string(variable)), + ) + + for (series_idx, ts) in enumerate(series) + data = _data(ts) + haskey(data, variable) || continue + lines!(ax, _times(ts), data[variable]; color=colors[series_idx], linewidth=2) + end + end + + return _save_if_requested(fig, figure_path) +end + +function _safe_fraction(numerator, denominator) + return numerator ./ ifelse.(denominator .> 0, denominator, NaN) +end + +function _stacked_area!(ax, times, values, labels, colors) + lower = zeros(length(times)) + + for (series, label, color) in zip(values, labels, colors) + upper = lower .+ series + band!(ax, times, lower, upper; color=color, label=label) + lower = upper + end + + return ax +end + +function plot_contributions(times, data; figure_path=nothing) + phyto, zoo, _ = _plankton_keys(data) + nutrient = haskey(data, :N) ? data[:N] : zeros(length(times)) + detritus = haskey(data, :D) ? data[:D] : zeros(length(times)) + phyto_total = _sum_keys(data, phyto) + zoo_total = _sum_keys(data, zoo) + living = phyto_total .+ zoo_total + nonliving = nutrient .+ detritus + total = living .+ nonliving + + fig = Figure(; size=(580, 410), fontsize=12) + + legend_elements = [ + PolyElement(; color=WORKSHOP_COLORS.living), + PolyElement(; color=WORKSHOP_COLORS.nonliving), + PolyElement(; color=WORKSHOP_COLORS.phytoplankton), + PolyElement(; color=WORKSHOP_COLORS.zooplankton), + ] + Legend( + fig[1, 1], + legend_elements, + ["living", "non-living", "phytoplankton", "zooplankton"]; + orientation=:horizontal, + tellwidth=false, + tellheight=true, + framevisible=false, + ) + + ax_living = Axis( + fig[2, 1]; + title="Living vs non-living nitrogen", + xlabel="Days", + ylabel="Fraction of total nitrogen", + limits=(nothing, nothing, 0, 1), + ) + + ax_plankton = Axis( + fig[3, 1]; + title="Phytoplankton vs zooplankton", + xlabel="Days", + ylabel="Fraction of living biomass", + limits=(nothing, nothing, 0, 1), + ) + + _stacked_area!( + ax_living, + times, + [_safe_fraction(living, total), _safe_fraction(nonliving, total)], + ["living", "non-living"], + [WORKSHOP_COLORS.living, WORKSHOP_COLORS.nonliving], + ) + + _stacked_area!( + ax_plankton, + times, + [_safe_fraction(phyto_total, living), _safe_fraction(zoo_total, living)], + ["phytoplankton", "zooplankton"], + [WORKSHOP_COLORS.phytoplankton, WORKSHOP_COLORS.zooplankton], + ) + + return _save_if_requested(fig, figure_path) +end + +function _as_bgc_vector(bgcs, n_series) + if bgcs isa AbstractVector || bgcs isa Tuple + bgcs = collect(bgcs) + length(bgcs) == n_series || error("bgcs must have one entry per box-model time series.") + return bgcs + else + return fill(bgcs, n_series) + end +end + +function plot_cwm_size( + timeseries, + bgcs; + labels=nothing, + figure_path=nothing, +) + series = _as_timeseries_vector(timeseries) + bgcs = _as_bgc_vector(bgcs, length(series)) + + if labels === nothing + labels = length(series) == 1 ? ["box model"] : ["series $i" for i in eachindex(series)] + end + + length(labels) == length(series) || error("labels must have one entry per box-model time series.") + + has_comparison = length(series) > 1 + legend_rows = has_comparison ? 1 : 0 + colors = _case_palette(length(series)) + fig = Figure(; size=(560, 420 + 44 * legend_rows), fontsize=12) + + if has_comparison + _shared_line_legend!(fig, 1, 1, colors, labels) + end + + axes = [ + Axis(fig[1 + legend_rows, 1]; title="All plankton", xlabel="Days", ylabel="CWM ESD (μm)"), + Axis(fig[2 + legend_rows, 1]; title="Phytoplankton", xlabel="Days", ylabel="CWM ESD (μm)"), + Axis(fig[3 + legend_rows, 1]; title="Zooplankton", xlabel="Days", ylabel="CWM ESD (μm)"), + ] + + for (idx, (ts, bgc)) in enumerate(zip(series, bgcs)) + data = _data(ts) + diameters = _plankton_diameters_by_tracer(bgc) + phyto, zoo, plankton = _plankton_keys(data) + + lines!(axes[1], _times(ts), _cwm_or_nan(data, diameters, plankton); color=colors[idx], linewidth=2) + lines!(axes[2], _times(ts), _cwm_or_nan(data, diameters, phyto); color=colors[idx], linewidth=2) + lines!(axes[3], _times(ts), _cwm_or_nan(data, diameters, zoo); color=colors[idx], linewidth=2) + end + + return _save_if_requested(fig, figure_path) +end