From 37183b13b9d1b79ddf377f5b93dfb7627188906d Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 17:17:16 +0100 Subject: [PATCH 01/33] overhauled particles... --- src/Models/Models.jl | 7 +- src/OceanBioME.jl | 2 +- src/Particles/Particles.jl | 97 ++++++++++++++++- src/Particles/advection.jl | 22 ++++ src/Particles/tendencies.jl | 35 ++++++ src/Particles/time_stepping.jl | 43 ++++++++ src/Particles/tracer_interpolation.jl | 56 ++++++++++ src/Particles/tracer_tendencies.jl | 4 - src/Particles/update_tracer_tendencies.jl | 39 +++++++ test/runtests.jl | 3 +- test/test_particles.jl | 127 ++++++++++++++++++++++ 11 files changed, 420 insertions(+), 15 deletions(-) create mode 100644 src/Particles/advection.jl create mode 100644 src/Particles/tendencies.jl create mode 100644 src/Particles/time_stepping.jl create mode 100644 src/Particles/tracer_interpolation.jl delete mode 100644 src/Particles/tracer_tendencies.jl create mode 100644 src/Particles/update_tracer_tendencies.jl create mode 100644 test/test_particles.jl diff --git a/src/Models/Models.jl b/src/Models/Models.jl index f7269e855..5d51b7aec 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -6,7 +6,7 @@ export NPZD, NutrientPhytoplanktonZooplanktonDetritus, LOBSTER -export SLatissima +export SugarKelp export CarbonChemistry @@ -22,7 +22,8 @@ export GasExchange, include("Sediments/Sediments.jl") include("AdvectedPopulations/LOBSTER/LOBSTER.jl") include("AdvectedPopulations/NPZD.jl") -include("Individuals/SLatissima.jl") +#include("Individuals/SLatissima.jl") +include("Individuals/SugarKelp.jl") include("seawater_density.jl") include("CarbonChemistry/CarbonChemistry.jl") include("GasExchange/GasExchange.jl") @@ -30,7 +31,7 @@ include("GasExchange/GasExchange.jl") using .Sediments using .LOBSTERModel using .NPZDModel -using .SLatissimaModel +using .SugarKelpModel using .CarbonChemistryModel using .GasExchangeModel diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 335a8c865..162cf027b 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -8,7 +8,7 @@ module OceanBioME export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, NPZD, redfield # Macroalgae models -export SLatissima +export SugarKelp # Box model export BoxModel, BoxModelGrid, SpeedyOutput, load_output diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 910b77d67..333307c1a 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -1,17 +1,98 @@ module Particles +using Adapt + +using KernelAbstractions: @kernel, @index + using Oceananigans: NonhydrostaticModel, HydrostaticFreeSurfaceModel + using OceanBioME: Biogeochemistry +using Oceananigans.Architectures: architecture, on_architecture + import Oceananigans.Biogeochemistry: update_tendencies! +import Oceananigans.Fields: set! import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, step_lagrangian_particles! import Oceananigans.OutputWriters: fetch_output import Base: length, size, show, summary +import Adapt: adapt_structure + +struct BiogeochemicalParticles{N, B, A, F, T, I, S, X, Y, Z} + biogeochemistry :: B + advection :: A + fields :: F + timestepper :: T + field_interpolation :: I + scalefactors :: S + x :: X # to maintain compatability with lagrangian particle advection from Oceananigans + y :: Y + z :: Z +end + +function required_particle_fields end +function required_tracers end +function coupled_tracers end + +@inline required_particle_fields(p::BiogeochemicalParticles) = required_particle_fields(p.biogeochemistry) +@inline required_tracers(p::BiogeochemicalParticles) = required_tracers(p.biogeochemistry) +@inline coupled_tracers(p::BiogeochemicalParticles) = coupled_tracers(p.biogeochemistry) + +include("advection.jl") +include("tracer_interpolation.jl") +include("tendencies.jl") +include("time_stepping.jl") +include("update_tracer_tendencies.jl") + +function BiogeochemicalParticles(number; + grid, + biogeochemistry::B, + advection::A = LagrangianAdvection(), + timestepper = ForwardEuler, + field_interpolation::I = NearestPoint(), + scalefactors = ones(number)) where {B, A, I} -abstract type BiogeochemicalParticles end + arch = architecture(grid) -# TODO: add model.particles passing + particle_fields = required_particle_fields(biogeochemistry) + fields = NamedTuple{particle_fields}(ntuple(n->on_architecture(arch, zeros(number)), Val(length(particle_fields)))) + + x = on_architecture(arch, zeros(number)) + y = on_architecture(arch, zeros(number)) + z = on_architecture(arch, zeros(number)) + + timestepper = timestepper(particle_fields, number, arch) + + scalefactors = on_architecture(arch, scalefactors) + + F = typeof(fields) + T = typeof(timestepper) + S = typeof(scalefactors) + X = typeof(x) + Y = typeof(y) + Z = typeof(z) + + return BiogeochemicalParticles{number, B, A, F, T, I, S, X, Y, Z}(biogeochemistry, + advection, + fields, + timestepper, + field_interpolation, + scalefactors, + x, y, z) +end + +Adapt.adapt_structure(to, p::BiogeochemicalParticles) = + BiogeochemicalParticles(adapt(to, p.biogeochemistry), + nothing, + nothing, + nothing, + adapt(to, p.field_interpolation), + nothing, + adapt(to, p.x), + adapt(to, p.y), + adapt(to, p.z)) + +# Type piracy...oops @inline step_lagrangian_particles!(::Nothing, model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}}, @@ -26,11 +107,15 @@ abstract type BiogeochemicalParticles end @inline update_lagrangian_particle_properties!(model, bgc::Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, Δt) = update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) -@inline update_lagrangian_particle_properties!(::BiogeochemicalParticles, model, bgc, Δt) = nothing +@inline function update_lagrangian_particle_properties!(particles::BiogeochemicalParticles, model, bgc, Δt) + advect_particles!(particles.advection, particles, model, Δt) + time_step_particle_fields!(particles.timestepper, particles, model, Δt) +end -size(particles::BiogeochemicalParticles) = size(particles.x) -length(particles::BiogeochemicalParticles) = length(particles.x) +size(particles::BiogeochemicalParticles) = (length(particles), ) +length(::BiogeochemicalParticles{N}) where N = N +# todo, fix below Base.summary(particles::BiogeochemicalParticles) = string(length(particles), " BiogeochemicalParticles with eltype ", nameof(eltype(particles)), " and properties ", propertynames(particles)) @@ -50,5 +135,5 @@ function fetch_output(particles::BiogeochemicalParticles, model) return NamedTuple{names}([getproperty(particles, name) for name in names]) end -include("tracer_tendencies.jl") +# TODO: implement set! end #module diff --git a/src/Particles/advection.jl b/src/Particles/advection.jl new file mode 100644 index 000000000..1ac64208f --- /dev/null +++ b/src/Particles/advection.jl @@ -0,0 +1,22 @@ +using Oceananigans.Architectures: on_architecture, device, architecture +using Oceananigans.Biogeochemistry: required_biogeochemical_tracers, biogeochemical_auxiliary_fields +using Oceananigans.Models.LagrangianParticleTracking: _advect_particles!, total_velocities + +# put nothing to do nothing +advect_particles!(advection, particles, model, Δt) = nothing + +struct LagrangianAdvection end + +function advect_particles!(::LagrangianAdvection, particles, model, Δt) + workgroup = min(length(particles), 256) + worksize = length(particles) + + arch = architecture(model) + + # Advect particles + advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) + + advect_particles_kernel!(particles, 1.0, model.grid, Δt, total_velocities(model)) + + return nothing +end diff --git a/src/Particles/tendencies.jl b/src/Particles/tendencies.jl new file mode 100644 index 000000000..4c3e2e2be --- /dev/null +++ b/src/Particles/tendencies.jl @@ -0,0 +1,35 @@ +using Oceananigans: fields + +function compute_particle_tendencies!(particles::BiogeochemicalParticles{N}, model) where N + tendencies = particles.timestepper.tendencies + + field_names = required_particle_fields(particles) + + dev = device(architecture(model)) + + workgroup = min(N, 256) + worksize = N + + tendency_kernel! = _tendency_kernel!(dev, workgroup, worksize) + + for name in field_names + tendency_kernel!(Val(name), particles, model.grid, fields(model), model.clock, tendencies) + end + + return nothing +end + +@kernel function _tendency_kernel!(val_field_name::Val{field_name}, particles, grid, fields, clock, tendencies) where field_name + n = @index(Global) + + t = clock.time + + field_names = required_particle_fields(particles) + nf = length(field_names) + + field_values = ntuple(n->particles.fields[n], Val(nf)) + + tracer_values = extract_tracer_values(particles.field_interpolation, particles, grid, fields, n) + + tendencies[field_name][n] = particles.biogeochemistry(val_field_name, t, field_values..., tracer_values...) +end diff --git a/src/Particles/time_stepping.jl b/src/Particles/time_stepping.jl new file mode 100644 index 000000000..0217b26fb --- /dev/null +++ b/src/Particles/time_stepping.jl @@ -0,0 +1,43 @@ +struct ForwardEuler{T} + tendencies :: T + + function ForwardEuler(field_names, number, arch) + tendencies = NamedTuple{field_names}(ntuple(n->on_architecture(arch, zeros(number)), Val(length(field_names)))) + + return new{typeof(tendencies)}(tendencies) + end +end + +# a different timestepper might want to loop over these etc. +function time_step_particle_fields!(timestepper::ForwardEuler, particles, model, Δt) + compute_particle_tendencies!(particles, model) + + step_particle_biogeochemistry!(timestepper, particles, model, Δt) + + return nothing +end + +function step_particle_biogeochemistry!(timestepper::ForwardEuler, particles::BiogeochemicalParticles{N}, model, Δt) where N + tendencies = timestepper.tendencies + + field_names = required_particle_fields(particles) + + dev = device(architecture(model)) + + workgroup = min(N, 256) + worksize = N + + step_kernel! = _euler_step!(dev, workgroup, worksize) + + for name in field_names + step_kernel!(particles.fields[name], tendencies[name], Δt) + end + + return nothing +end + +@kernel function _euler_step!(field, tendency, Δt) + n = @index(Global) + + @inbounds field[n] += tendency[n] * Δt +end \ No newline at end of file diff --git a/src/Particles/tracer_interpolation.jl b/src/Particles/tracer_interpolation.jl new file mode 100644 index 000000000..fe5504c77 --- /dev/null +++ b/src/Particles/tracer_interpolation.jl @@ -0,0 +1,56 @@ +using Oceananigans.Operators: volume +using Oceananigans.Fields: fractional_indices, _interpolate, interpolator, Center +using Oceananigans.Grids: AbstractGrid, Flat, Bounded, Periodic + +@inline get_node(::Bounded, i, N) = min(max(i, 1), N) +@inline get_node(::Periodic, i, N) = ifelse(i < 1, N, ifelse(i > N, 1, i)) + +struct NearestPoint end + +@inline function extract_tracer_values(::NearestPoint, particles, grid, fields, n) + x = @inbounds particles.x[n] + y = @inbounds particles.y[n] + z = @inbounds particles.z[n] + + i, j, k = nearest_node(x, y, z, grid) + + field_names = required_tracers(particles) + nf = length(field_names) + + field_values = ntuple(Val(nf)) do n + fields[field_names[n]][i, j, k] + end + + return field_values +end + +@inline function nearest_node(x, y, z, grid::AbstractGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} + # messy + ii, jj, kk = fractional_indices((x, y, z), grid, ifelse(isa(TX(), Flat), nothing, Center()), + ifelse(isa(TY(), Flat), nothing, Center()), + ifelse(isa(TZ(), Flat), nothing, Center())) + + ix = interpolator(ii) + iy = interpolator(jj) + iz = interpolator(kk) + + i, j, k = (get_node(TX(), Int(ifelse(ix[3] < 0.5, ix[1], ix[2])), grid.Nx), + get_node(TY(), Int(ifelse(iy[3] < 0.5, iy[1], iy[2])), grid.Ny), + get_node(TZ(), Int(ifelse(iz[3] < 0.5, iz[1], iz[2])), grid.Nz)) + + return i, j, k +end + +@inline function apply_tracer_tendency!(::NearestPoint, particles, grid, particle_tendency, tendency, n) + x = @inbounds particles.x[n] + y = @inbounds particles.y[n] + z = @inbounds particles.z[n] + + i, j, k = nearest_node(x, y, z, grid) + + node_volume = volume(i, j, k, grid, Center(), Center(), Center()) + + @inbounds tendency[i, j, k] += particle_tendency / node_volume + + return nothing +end \ No newline at end of file diff --git a/src/Particles/tracer_tendencies.jl b/src/Particles/tracer_tendencies.jl deleted file mode 100644 index 8aeaa4553..000000000 --- a/src/Particles/tracer_tendencies.jl +++ /dev/null @@ -1,4 +0,0 @@ -using Oceananigans.Grids: AbstractGrid, Bounded, Periodic - -@inline get_node(::Bounded, i, N) = min(max(i, 1), N) -@inline get_node(::Periodic, i, N) = ifelse(i < 1, N, ifelse(i > N, 1, i)) diff --git a/src/Particles/update_tracer_tendencies.jl b/src/Particles/update_tracer_tendencies.jl new file mode 100644 index 000000000..2c9f4eb9f --- /dev/null +++ b/src/Particles/update_tracer_tendencies.jl @@ -0,0 +1,39 @@ +function update_tendencies!(bgc, particles::BiogeochemicalParticles{N}, model) where N + tendencies = particles.timestepper.tendencies + + field_names = coupled_tracers(particles) + + dev = device(architecture(model)) + + workgroup = min(N, 256) + worksize = N + + tendency_kernel! = _update_tracer_tendencies!(dev, workgroup, worksize) + + for name in field_names + tendency_kernel!(Val(name), particles, model.grid, fields(model), model.clock, model.timestepper.Gⁿ, particles.scalefactors) + end + + return nothing +end + +@kernel function _update_tracer_tendencies!(val_field_name::Val{field_name}, particles, grid, fields, clock, tendencies, scalefactors) where field_name + n = @index(Global) + + t = clock.time + + field_names = required_particle_fields(particles) + nf = length(field_names) + + field_values = ntuple(n->particles.fields[n], Val(nf)) + + tracer_values = extract_tracer_values(particles.field_interpolation, particles, grid, fields, n) + + particle_tendency = particles.biogeochemistry(val_field_name, t, field_values..., tracer_values...) + + scalefactor = @inbounds scalefactors[n] + + total_tendency = @show scalefactor * particle_tendency + + apply_tracer_tendency!(particles.field_interpolation, particles, grid, total_tendency, tendencies[field_name], n) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 44cde13f6..66e2c4281 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,8 @@ include("dependencies_for_runtests.jl") include("test_utils.jl") include("test_light.jl") -include("test_slatissima.jl") +include("test_particles.jl") +#include("test_slatissima.jl") include("test_LOBSTER.jl") include("test_NPZD.jl") include("test_gasexchange_carbon_chem.jl") diff --git a/test/test_particles.jl b/test/test_particles.jl new file mode 100644 index 000000000..176ca7213 --- /dev/null +++ b/test/test_particles.jl @@ -0,0 +1,127 @@ +#include("dependencies_for_runtests.jl") + +using OceanBioME.Particles: BiogeochemicalParticles + +import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers + +struct SimpleParticleBiogeochemistry end + +@inline (::SimpleParticleBiogeochemistry)(val_name, t, A) = one(t) + +@inline required_particle_fields(::SimpleParticleBiogeochemistry) = (:A, ) + +@inline required_tracers(::SimpleParticleBiogeochemistry) = tuple() +@inline coupled_tracers(::SimpleParticleBiogeochemistry) = tuple() + +grid = RectilinearGrid(architecture; size = (3, 3, 3), extent = (3, 3, 3)) + +@testset "Testing basic particle property integration" begin + #for TX in (Bounded, Flat) # to check it works with different topologies + + particle_biogeochemistry = SimpleParticleBiogeochemistry() + + particles = BiogeochemicalParticles(3; grid, biogeochemistry = particle_biogeochemistry, advection = nothing) + + @test length(particles) == 3 + + biogeochemistry = Biogeochemistry(nothing; particles) + + model = NonhydrostaticModel(; grid, biogeochemistry) + + time_step!(model, 1) + + @test all(particles.fields.A .== 1) +end + +@inline required_tracers(::SimpleParticleBiogeochemistry) = (:B, ) + +@inline (::SimpleParticleBiogeochemistry)(val_name, t, A, B) = 0.1 * B + +@testset "Testing particle tracer detection" begin + particle_biogeochemistry = SimpleParticleBiogeochemistry() + + particles = BiogeochemicalParticles(3; grid, biogeochemistry = particle_biogeochemistry, advection = nothing) + + biogeochemistry = Biogeochemistry(nothing; particles) + + model = NonhydrostaticModel(; grid, biogeochemistry, tracers = :B) + + set!(model, B = 1) + + time_step!(model, 1) + + @test all(particles.fields.A .== 0.1) + + # also these shouldn't have moved + @test all(particles.x .== 0) +end + +@testset "Testing particle lagrangian advection" begin + particle_biogeochemistry = SimpleParticleBiogeochemistry() + + particles = BiogeochemicalParticles(3; grid, biogeochemistry = particle_biogeochemistry) + + biogeochemistry = Biogeochemistry(nothing; particles) + + model = NonhydrostaticModel(; grid, biogeochemistry, tracers = :B) + + set!(model, B = 1, u = 0.1, v = 0.2) + + time_step!(model, 1) + + @test all(particles.fields.A .== 0.1) + @test all(particles.x .== 0.1) && all(particles.y .== 0.2) && all(particles.z .== 0) +end + +coupled_tracers(::SimpleParticleBiogeochemistry) = (:B, ) + +@inline (::SimpleParticleBiogeochemistry)(::Val{:A}, t, A, B) = 0.1 +@inline (::SimpleParticleBiogeochemistry)(::Val{:B}, t, A, B) = -0.1 + +@testset "Testing particle-tracer uptake" begin + particle_biogeochemistry = SimpleParticleBiogeochemistry() + + particles = BiogeochemicalParticles(3; grid, biogeochemistry = particle_biogeochemistry, advection = nothing) + + biogeochemistry = Biogeochemistry(nothing; particles) + + model = NonhydrostaticModel(; grid, biogeochemistry, tracers = :B) + + set!(model, B = 1) + + # since the 0, 0, 0 point is ambigously closest to both 3, 3, 3 and 1, 1, 3 (and the logic makes it go to 3, 3, 3) + particles.x .= 0.5 + particles.y .= 0.5 + + time_step!(model, 1) + + @test all(particles.fields.A .== 0.1) + + @test interior(model.tracers.B, 1, 1, 3) .== 0.7 # 1 - 3 * 0.1 +end + + +@testset "Testing particle-tracer uptake with scalefactors" begin + particle_biogeochemistry = SimpleParticleBiogeochemistry() + + particles = BiogeochemicalParticles(3; grid, + biogeochemistry = particle_biogeochemistry, + advection = nothing, + scalefactors = [1, 0.5, 0.5]) + + biogeochemistry = Biogeochemistry(nothing; particles) + + model = NonhydrostaticModel(; grid, biogeochemistry, tracers = :B) + + set!(model, B = 1) + + # since the 0, 0, 0 point is ambigously closest to both 3, 3, 3 and 1, 1, 3 (and the logic makes it go to 3, 3, 3) + particles.x .= 0.5 + particles.y .= 0.5 + + time_step!(model, 1) + + @test all(particles.fields.A .== 0.1) + + @test interior(model.tracers.B, 1, 1, 3) .== 0.8 +end \ No newline at end of file From 0fc7745537898aa6a73fc8ee3d029b8e29bf5380 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 17:19:09 +0100 Subject: [PATCH 02/33] first attempt --- src/Models/Individuals/SugarKelp.jl | 212 ++++++++++++++++++++++++++++ 1 file changed, 212 insertions(+) create mode 100644 src/Models/Individuals/SugarKelp.jl diff --git a/src/Models/Individuals/SugarKelp.jl b/src/Models/Individuals/SugarKelp.jl new file mode 100644 index 000000000..f3261563b --- /dev/null +++ b/src/Models/Individuals/SugarKelp.jl @@ -0,0 +1,212 @@ +module SugarKelpModel + +using Roots + +import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers + +# seriously disgsuting number of parameters +@kwdef struct SugarKelp{FT} + growth_rate_adjustment :: FT = 4.5 + photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60) + minimum_carbon_reserve :: FT = 0.01 + structural_carbon :: FT = 0.2 + exudation :: FT = 0.5 + erosion :: FT = 0.22 + saturation_irradiance :: FT = 90 * day/ (10 ^ 6) + structural_dry_weight_per_area :: FT = 0.5 + structural_dry_to_wet_weight :: FT = 0.0785 + carbon_reserve_per_carbon :: FT = 2.1213 + nitrogen_reserve_per_nitrogen :: FT = 2.72 + minimum_nitrogen_reserve :: FT = 0.0126 + maximum_nitrogen_reserve :: FT = 0.0216 + growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) + growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2 + maximum_specific_growth_rate :: FT = 0.18 + structural_nitrogen :: FT = 0.0146 + photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24 + photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24 + photosynthesis_ref_temp_1 :: FT = 285.0 + photosynthesis_ref_temp_2 :: FT = 288.0 + photoperiod_1 :: FT = 0.85 + photoperiod_2 :: FT = 0.3 + respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24 + respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24 + respiration_ref_temp_1 :: FT = 285.0 + respiration_ref_temp_2 :: FT = 290.0 + photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1) + photosynthesis_low_temp :: FT = 271.0 + photosynthesis_high_temp :: FT = 296.0 + photosynthesis_high_arrhenius_temp :: FT = 1414.87 + photosynthesis_low_arrhenius_temp :: FT = 4547.89 + respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1) + current_speed_for_0p65_uptake :: FT = 0.03 + nitrate_half_saturation :: FT = 4.0 + ammonia_half_saturation :: FT = 1.3 + maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6) + maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6) + current_1 :: FT = 0.72 + current_2 :: FT = 0.28 + current_3 :: FT = 0.045 + base_activity_respiration_rate :: FT = 1.11e-4 * 24 + base_basal_respiration_rate :: FT = 5.57e-5 * 24 + exudation_redfield_ratio :: FT = Inf +end + +@inline required_particle_fields(::SugarKelp) = (:A, :N, :C) +@inline required_tracers(::SugarKelp) = (:u, :v, :w, :NO₃, :NH₄) +@inline coupled_tracers(::SugarKelp) = (:NO₃, :NH₄) + +@inline function (kelp::SugarKelp)(::Val{:A}, t, A, N, C, u, v, w, NO₃, NH₄) + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) + + ν = errosion(kelp, t, A, N, C, T) + + return A * (μ - ν) +end + +@inline function (kelp::SugarKelp)(::Val{:N}, t, A, N, C, u, v, w, NO₃, NH₄) + kₐ = kelp.structural_dry_weight_per_area + Nₛ = kelp.structural_nitrogen + + J = nitrate_uptake(kelp, N, NO₃, u, v, w) + + ammonia_uptake(kelp, t, N, C, T, NH₄, u, v, w) + + e = nitrogen_exudate(kelp, C, T, PAR) + + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) + + return (J - e) / kₐ - μ * (N + Nₛ) +end + +@inline function (kelp::SugarKelp)(::Val{:C}, t, A, N, C, u, v, w, NO₃, NH₄) + kₐ = kelp.structural_dry_weight_per_area + Cₛ = kelp.structural_carbon + + P = photosynthesis(kelp, T, PAR) + + R = respiration(kelp, t, N, C, T, NH₄, u, v, w, μ) + + e = specific_carbon_exudate(kelp, C) + + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) + + return (P * (1 - e) - R) / kₐ - μ * (C + Cₛ) +end + +@inline function growth(kelp, t, A, N, C, T, NH₄, u, v, w) + f = base_growth_limitation(kelp, t, A, N, C, T) + + # potential ammonia based_growth + kₐ = kelp.structural_dry_weight_per_area + Nₛ = kelp.structural_nitrogen + + j̃_NH₄ = potential_ammonia_uptake(kelp, NH₄, u, v, w) + + μNH₄ = j̃_NH₄ / kₐ / (N + Nₛ) + + # internal reserve based growth + μN = 1 - Nₘ / N + μC = 1 - Cₘ / C + + return f * min(μC, max(μN, μNH₄)) +end + +@inline function nitrate_uptake(kelp, N, NO₃, u, v, w) + k_NO₃ = kelp.nitrate_half_saturation + + N_max = kelp.maximum_nitrogen_reserve + N_min = kelp.minimum_nitrogen_reserve + + fᶜ = current_factor(kelp, u, v, w) + + return max(0, J_max_NO₃ * fᶜ * (N_max - N) / (N_max - N_min) * NO₃ / (k_NO₃ + NO₃)) +end + +@inline function ammonia_uptake(kelp, t, N, C, T, NH₄, u, v, w) + kₐ = kelp.structural_dry_weight_per_area + Nₛ = kelp.structural_nitrogen + + j̃_NH₄ = potential_ammonia_uptake(kelp, NH₄, u, v, w) + + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) + + return min(j̃_NH₄, μ * kₐ * (N + Nₛ)) +end + +@inline function potential_ammonia_uptake(kelp, NH₄, u, v, w) + fᶜ = current_factor(kelp, u, v, w) + + return J_max_NH₄ * fᶜ * NH₄ / (k_NH₄ + NH₄) +end + +@inline function photosynthesis(kelp, T, PAR) + Tk = T + 273.15 + + P₁ = kelp.photosynthesis_at_ref_temp_1 + + Tₐ = kelp.photosynthesis_arrhenius_temp + Tₐₗ = kelp.photosynthesis_low_arrhenius_temp + Tₐₕ = kelp.photosynthesis_high_arrhenius_temp + + Tₚ = kelp.photosynthesis_ref_temp_1 + Tₚₗ = kelp.photosynthesis_ref_temp_1 + Tₚₕ = kelp.photosynthesis_high_temp + + maximum_photosynthesis = P₁ * exp(Tₐ / Tₚ - Tₐ / Tk) / (1 + exp(Tₐₗ / Tk - Tₐₗ / Tₚₗ) + exp(Tₐₕ / Tₚₕ - Tₐₕ / Tk)) + + β = solve_for_light_inhibition(kelp, maximum_photosynthesis) + + pₛ = p.photosynthetic_efficiency * p.saturation_irradiance / log(1 + p.photosynthetic_efficiency / β) + + return pₛ * (1 - exp(- p.photosynthetic_efficiency * PAR / pₛ)) * exp(-β * PAR / pₛ) +end + +# solves `alkalinity_residual` for pH +@inline solve_for_light_inhibition(kelp, maximum_photosynthesis) = + find_zero(light_inhibition_residual, (0, 0.1), Bisection(); + p = (; maximum_photosynthesis, kelp.photosynthetic_efficiency, kelp.saturation_irradiance)) + +@inline function light_inhibition_residual(β, p) + pₘ = p.maximum_photosynthesis + α = p.photosynthetic_efficiency + Iₛ = p.saturation_irradiance + + return pₘ - α * Iₛ / log(1 + α / β) * (α / (α + β)) * (β / (α + β)) ^ (β / α) +end + +@inline function respiration(kelp, t, N, C, T, NH₄, u, v, w, μ) + Rᵇ = kelp.base_basal_respiration_rate + Rᵃ = kelp.base_activity_respiration_rate + + Tₐ = kelp.respiration_arrhenius_temp + T₁ = kelp.respiration_ref_temp_1 + + Tk = T + 273.15 + + f = exp(Tₐ / T₁ - Tₐ / Tk) + + μₘ = kelp.maximum_specific_growth_rate + Jₘ = kelp.maximum_nitrate_uptake + kelp.maximum_ammonia_uptake + + J = nitrate_uptake(kelp, N, NO₃, u, v, w) + + ammonia_uptake(kelp, t, N, C, T, NH₄, u, v, w) + # (basal + activity associated) * temp factor + return f * (Rᵇ + Rᵃ * (μ / μₘ + J / Jₘ)) +end + +@inline function specific_carbon_exudate(kelp, C) + γ = exudation + Cₘ = minimum_carbon_reserve + + return 1 - exp(γ * (Cₘ - C)) +end + +@inline function nitrogen_exudate(kelp, C, T, PAR) + CN = kelp.exudation_redfield_ratio + + P = photosynthesis(kelp, T, PAR) + e = specific_carbon_exudate(kelp, C) + + return P * e * 14 / 12 / CN +end +end # module \ No newline at end of file From 3453447236e514e3935663ac50ac6e9a0637ece0 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 19:59:39 +0100 Subject: [PATCH 03/33] sugar kelp steps again! --- src/Models/Individuals/SugarKelp/SugarKelp.jl | 71 +++++++ src/Models/Individuals/SugarKelp/coupling.jl | 57 ++++++ .../{SugarKelp.jl => SugarKelp/equations.jl} | 191 ++++++++++-------- src/Models/Individuals/SugarKelp/show.jl | 8 + src/Models/Models.jl | 2 +- src/OceanBioME.jl | 3 +- src/Particles/Particles.jl | 18 +- src/Particles/atomic_operations.jl | 24 +++ src/Particles/tendencies.jl | 4 +- src/Particles/tracer_interpolation.jl | 2 +- src/Particles/update_tracer_tendencies.jl | 29 ++- test/runtests.jl | 5 +- test/test_sugar_kelp.jl | 26 +++ 13 files changed, 336 insertions(+), 104 deletions(-) create mode 100644 src/Models/Individuals/SugarKelp/SugarKelp.jl create mode 100644 src/Models/Individuals/SugarKelp/coupling.jl rename src/Models/Individuals/{SugarKelp.jl => SugarKelp/equations.jl} (51%) create mode 100644 src/Models/Individuals/SugarKelp/show.jl create mode 100644 src/Particles/atomic_operations.jl create mode 100644 test/test_sugar_kelp.jl diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl new file mode 100644 index 000000000..b2f1a8bc6 --- /dev/null +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -0,0 +1,71 @@ +module SugarKelpModel + +using Roots + +using Oceananigans.Units + +import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers + +# disgsuting number of parameters +@kwdef struct SugarKelp{FT, TL} + temperature_limit :: TL = LinearOptimalTemperatureRange() # TODO: split up more parameterisations like this + growth_rate_adjustment :: FT = 4.5 + photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60) + minimum_carbon_reserve :: FT = 0.01 + structural_carbon :: FT = 0.2 + exudation :: FT = 0.5 + erosion_exponent :: FT = 0.22 + base_erosion_rate :: FT = 10^-6 + saturation_irradiance :: FT = 90 * day/ (10 ^ 6) + structural_dry_weight_per_area :: FT = 0.5 + structural_dry_to_wet_weight :: FT = 0.0785 + carbon_reserve_per_carbon :: FT = 2.1213 + nitrogen_reserve_per_nitrogen :: FT = 2.72 + minimum_nitrogen_reserve :: FT = 0.0126 + maximum_nitrogen_reserve :: FT = 0.0216 + growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) + growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2 + maximum_specific_growth_rate :: FT = 0.18 + structural_nitrogen :: FT = 0.0146 + photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24 + photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24 + photosynthesis_ref_temp_1 :: FT = 285.0 + photosynthesis_ref_temp_2 :: FT = 288.0 + photoperiod_1 :: FT = 0.85 + photoperiod_2 :: FT = 0.3 + respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24 + respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24 + respiration_ref_temp_1 :: FT = 285.0 + respiration_ref_temp_2 :: FT = 290.0 + photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1) + photosynthesis_low_temp :: FT = 271.0 + photosynthesis_high_temp :: FT = 296.0 + photosynthesis_high_arrhenius_temp :: FT = 1414.87 + photosynthesis_low_arrhenius_temp :: FT = 4547.89 + respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1) + current_speed_for_0p65_uptake :: FT = 0.03 + nitrate_half_saturation :: FT = 4.0 + ammonia_half_saturation :: FT = 1.3 + maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6) + maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6) + current_1 :: FT = 0.72 + current_2 :: FT = 0.28 + current_3 :: FT = 0.045 + base_activity_respiration_rate :: FT = 1.11e-4 * 24 + base_basal_respiration_rate :: FT = 5.57e-5 * 24 + exudation_redfield_ratio :: FT = Inf + adapted_latitude :: FT = 57.5 +end + +@inline required_particle_fields(::SugarKelp) = (:A, :N, :C) +@inline required_tracers(::SugarKelp) = (:u, :v, :w, :T, :NO₃, :NH₄, :PAR) +@inline coupled_tracers(::SugarKelp) = (:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON) +# can overload like: +# @inline coupled_tracers(::SugarKelp) = (NO₃ = :NO₃, NH₄ = :NH₄, DON = :DOM, bPON = :POM) +# if you have a simpler model ... there must be a better way todo this + +include("equations.jl") +include("coupling.jl") +include("show.jl") + +end # module \ No newline at end of file diff --git a/src/Models/Individuals/SugarKelp/coupling.jl b/src/Models/Individuals/SugarKelp/coupling.jl new file mode 100644 index 000000000..9794c22e5 --- /dev/null +++ b/src/Models/Individuals/SugarKelp/coupling.jl @@ -0,0 +1,57 @@ +#:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON + +@inline function (kelp::SugarKelp)(::Val{:NO₃}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + J = nitrate_uptake(kelp, N, NO₃, u, v, w) + + return J * A / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s +end + +@inline function (kelp::SugarKelp)(::Val{:NH₄}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + J = ammonia_uptake(kelp, t, A, N, C, T, NH₄, u, v, w) + + return J * A / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s +end + +@inline function (kelp::SugarKelp)(::Val{:DIC}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + P = photosynthesis(kelp, T, PAR) + + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) + + R = respiration(kelp, t, A, N, C, T, NO₃, NH₄, u, v, w, μ) + + return (P - R) * A / (day * 12 * 0.001) # gC/dm^2/hr to mmol C/s +end + +# I now know that this may not be correct..., it probably should vary by where the nitrogen comes from +@inline (kelp::SugarKelp)(::Val{:O₂}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) = + - kelp(Val(:DIC), t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + +@inline function (kelp::SugarKelp)(::Val{:DOC}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + P = photosynthesis(kelp, T, PAR) + + e = specific_carbon_exudate(kelp, C) + + return e * P * A / (day * 12 * 0.001) # gC/dm^2/hr to mmol C/s +end + +@inline (kelp::SugarKelp)(::Val{:DON}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) = + kelp(Val(:DOC), t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) / kelp.exudation_redfield_ratio + +@inline function (kelp::SugarKelp)(::Val{:bPON}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + kₐ = kelp.structural_dry_weight_per_area + Nₛ = kelp.structural_nitrogen + + ν = erosion(kelp, t, A, N, C, T) + + return ν * kₐ * A * (N + Nₛ) / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s +end + +# this is why the particles made by kelp can be much more carbon rich +@inline function (kelp::SugarKelp)(::Val{:bPOC}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) + kₐ = kelp.structural_dry_weight_per_area + Cₛ = kelp.structural_carbon + + ν = erosion(kelp, t, A, N, C, T) + + return ν * kₐ * A * (C + Cₛ) / (day * 14 * 0.001) # gC/dm^2/hr to mmol C/s +end diff --git a/src/Models/Individuals/SugarKelp.jl b/src/Models/Individuals/SugarKelp/equations.jl similarity index 51% rename from src/Models/Individuals/SugarKelp.jl rename to src/Models/Individuals/SugarKelp/equations.jl index f3261563b..b064f5a2a 100644 --- a/src/Models/Individuals/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/equations.jl @@ -1,96 +1,39 @@ -module SugarKelpModel - -using Roots - -import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers - -# seriously disgsuting number of parameters -@kwdef struct SugarKelp{FT} - growth_rate_adjustment :: FT = 4.5 - photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60) - minimum_carbon_reserve :: FT = 0.01 - structural_carbon :: FT = 0.2 - exudation :: FT = 0.5 - erosion :: FT = 0.22 - saturation_irradiance :: FT = 90 * day/ (10 ^ 6) - structural_dry_weight_per_area :: FT = 0.5 - structural_dry_to_wet_weight :: FT = 0.0785 - carbon_reserve_per_carbon :: FT = 2.1213 - nitrogen_reserve_per_nitrogen :: FT = 2.72 - minimum_nitrogen_reserve :: FT = 0.0126 - maximum_nitrogen_reserve :: FT = 0.0216 - growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2 - maximum_specific_growth_rate :: FT = 0.18 - structural_nitrogen :: FT = 0.0146 - photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24 - photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24 - photosynthesis_ref_temp_1 :: FT = 285.0 - photosynthesis_ref_temp_2 :: FT = 288.0 - photoperiod_1 :: FT = 0.85 - photoperiod_2 :: FT = 0.3 - respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24 - respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24 - respiration_ref_temp_1 :: FT = 285.0 - respiration_ref_temp_2 :: FT = 290.0 - photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1) - photosynthesis_low_temp :: FT = 271.0 - photosynthesis_high_temp :: FT = 296.0 - photosynthesis_high_arrhenius_temp :: FT = 1414.87 - photosynthesis_low_arrhenius_temp :: FT = 4547.89 - respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1) - current_speed_for_0p65_uptake :: FT = 0.03 - nitrate_half_saturation :: FT = 4.0 - ammonia_half_saturation :: FT = 1.3 - maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6) - maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6) - current_1 :: FT = 0.72 - current_2 :: FT = 0.28 - current_3 :: FT = 0.045 - base_activity_respiration_rate :: FT = 1.11e-4 * 24 - base_basal_respiration_rate :: FT = 5.57e-5 * 24 - exudation_redfield_ratio :: FT = Inf -end - -@inline required_particle_fields(::SugarKelp) = (:A, :N, :C) -@inline required_tracers(::SugarKelp) = (:u, :v, :w, :NO₃, :NH₄) -@inline coupled_tracers(::SugarKelp) = (:NO₃, :NH₄) - -@inline function (kelp::SugarKelp)(::Val{:A}, t, A, N, C, u, v, w, NO₃, NH₄) + +@inline function (kelp::SugarKelp)(::Val{:A}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) - ν = errosion(kelp, t, A, N, C, T) + ν = erosion(kelp, t, A, N, C, T) - return A * (μ - ν) + return A * (μ - ν) / day end -@inline function (kelp::SugarKelp)(::Val{:N}, t, A, N, C, u, v, w, NO₃, NH₄) +@inline function (kelp::SugarKelp)(::Val{:N}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) kₐ = kelp.structural_dry_weight_per_area Nₛ = kelp.structural_nitrogen J = nitrate_uptake(kelp, N, NO₃, u, v, w) + - ammonia_uptake(kelp, t, N, C, T, NH₄, u, v, w) + ammonia_uptake(kelp, t, A, N, C, T, NH₄, u, v, w) e = nitrogen_exudate(kelp, C, T, PAR) μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) - return (J - e) / kₐ - μ * (N + Nₛ) + return ((J - e) / kₐ - μ * (N + Nₛ)) / day end -@inline function (kelp::SugarKelp)(::Val{:C}, t, A, N, C, u, v, w, NO₃, NH₄) +@inline function (kelp::SugarKelp)(::Val{:C}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) kₐ = kelp.structural_dry_weight_per_area Cₛ = kelp.structural_carbon P = photosynthesis(kelp, T, PAR) + + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) - R = respiration(kelp, t, N, C, T, NH₄, u, v, w, μ) + R = respiration(kelp, t, A, N, C, T, NO₃, NH₄, u, v, w, μ) e = specific_carbon_exudate(kelp, C) - μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) - - return (P * (1 - e) - R) / kₐ - μ * (C + Cₛ) + return ((P * (1 - e) - R) / kₐ - μ * (C + Cₛ)) / day end @inline function growth(kelp, t, A, N, C, T, NH₄, u, v, w) @@ -99,19 +42,24 @@ end # potential ammonia based_growth kₐ = kelp.structural_dry_weight_per_area Nₛ = kelp.structural_nitrogen + + Nₘ = kelp.minimum_nitrogen_reserve + Cₘ = kelp.minimum_carbon_reserve j̃_NH₄ = potential_ammonia_uptake(kelp, NH₄, u, v, w) μNH₄ = j̃_NH₄ / kₐ / (N + Nₛ) # internal reserve based growth - μN = 1 - Nₘ / N - μC = 1 - Cₘ / C + μN = 1 - Nₘ / N + μC = 1 - Cₘ / C return f * min(μC, max(μN, μNH₄)) end @inline function nitrate_uptake(kelp, N, NO₃, u, v, w) + J_max_NO₃ = kelp.maximum_nitrate_uptake + k_NO₃ = kelp.nitrate_half_saturation N_max = kelp.maximum_nitrogen_reserve @@ -122,7 +70,7 @@ end return max(0, J_max_NO₃ * fᶜ * (N_max - N) / (N_max - N_min) * NO₃ / (k_NO₃ + NO₃)) end -@inline function ammonia_uptake(kelp, t, N, C, T, NH₄, u, v, w) +@inline function ammonia_uptake(kelp, t, A, N, C, T, NH₄, u, v, w) kₐ = kelp.structural_dry_weight_per_area Nₛ = kelp.structural_nitrogen @@ -134,6 +82,9 @@ end end @inline function potential_ammonia_uptake(kelp, NH₄, u, v, w) + J_max_NH₄ = kelp.maximum_ammonia_uptake + k_NH₄ = kelp.ammonia_half_saturation + fᶜ = current_factor(kelp, u, v, w) return J_max_NH₄ * fᶜ * NH₄ / (k_NH₄ + NH₄) @@ -148,17 +99,20 @@ end Tₐₗ = kelp.photosynthesis_low_arrhenius_temp Tₐₕ = kelp.photosynthesis_high_arrhenius_temp - Tₚ = kelp.photosynthesis_ref_temp_1 + Tₚ = kelp.photosynthesis_ref_temp_1 Tₚₗ = kelp.photosynthesis_ref_temp_1 Tₚₕ = kelp.photosynthesis_high_temp + α = kelp.photosynthetic_efficiency + Iₛ = kelp.saturation_irradiance + maximum_photosynthesis = P₁ * exp(Tₐ / Tₚ - Tₐ / Tk) / (1 + exp(Tₐₗ / Tk - Tₐₗ / Tₚₗ) + exp(Tₐₕ / Tₚₕ - Tₐₕ / Tk)) β = solve_for_light_inhibition(kelp, maximum_photosynthesis) - pₛ = p.photosynthetic_efficiency * p.saturation_irradiance / log(1 + p.photosynthetic_efficiency / β) + pₛ = α * Iₛ / log(1 + α / β) - return pₛ * (1 - exp(- p.photosynthetic_efficiency * PAR / pₛ)) * exp(-β * PAR / pₛ) + return pₛ * (1 - exp(- α * PAR / pₛ)) * exp(-β * PAR / pₛ) end # solves `alkalinity_residual` for pH @@ -174,7 +128,7 @@ end return pₘ - α * Iₛ / log(1 + α / β) * (α / (α + β)) * (β / (α + β)) ^ (β / α) end -@inline function respiration(kelp, t, N, C, T, NH₄, u, v, w, μ) +@inline function respiration(kelp, t, A, N, C, T, NO₃, NH₄, u, v, w, μ) Rᵇ = kelp.base_basal_respiration_rate Rᵃ = kelp.base_activity_respiration_rate @@ -189,14 +143,14 @@ end Jₘ = kelp.maximum_nitrate_uptake + kelp.maximum_ammonia_uptake J = nitrate_uptake(kelp, N, NO₃, u, v, w) + - ammonia_uptake(kelp, t, N, C, T, NH₄, u, v, w) + ammonia_uptake(kelp, t, A, N, C, T, NH₄, u, v, w) # (basal + activity associated) * temp factor return f * (Rᵇ + Rᵃ * (μ / μₘ + J / Jₘ)) end @inline function specific_carbon_exudate(kelp, C) - γ = exudation - Cₘ = minimum_carbon_reserve + γ = kelp.exudation + Cₘ = kelp.minimum_carbon_reserve return 1 - exp(γ * (Cₘ - C)) end @@ -209,4 +163,81 @@ end return P * e * 14 / 12 / CN end -end # module \ No newline at end of file + +@inline function erosion(kelp, t, A, N, C, T) + ε = kelp.erosion_exponent + ν₀ = kelp.base_erosion_rate + + return ν₀ * exp(ε * A) / (1 + ν₀ * (exp(ε * A) - 1)) +end + +@inline function current_factor(kelp, u, v, w) + f₁ = kelp.current_1 + f₂ = kelp.current_2 + f₃ = kelp.current_3 + + U = √(u^2 + v^2 + w^2) + + return f₁ * (1 - exp(-U / f₃)) + f₂ +end + +@inline function base_growth_limitation(kelp, t, A, N, C, T) + fₜ = kelp.temperature_limit(T) + fₐ = area_limitation(kelp, A) + fₚ = seasonal_limitation(kelp, t) +end + +@inline function area_limitation(kelp, A) + A₀ = kelp.growth_rate_adjustment + m₁ = kelp.growth_adjustment_1 + m₂ = kelp.growth_adjustment_2 + + return m₁ * exp(-(A / A₀)^2) + m₂ +end + +@kwdef struct LinearOptimalTemperatureRange{FT} + lower_optimal :: FT = 10.0 + upper_optimal :: FT = 15.0 + lower_gradient :: FT = 1/(lower_optimal + 1.8) #0.08 - I think its important that the minimum cut off (-1.8) is observed, + # because the literature doesn't say anything about it growing okay and then suddenly dying off + #at low temperature as the origional form suggests, although for higher temperature... + upper_gradient :: FT = -0.25 +end + +function (f::LinearOptimalTemperatureRange)(T) + Tₗ = f.lower_optimal + Tᵤ = f.upper_optimal + + αₗ = f.lower_gradient + αᵤ = f.upper_gradient + + return (max(0, αₗ * (T - Tₗ) + 1) * (T < Tₗ) + + max(0, αᵤ * (T - Tᵤ) + 1) * (T > Tᵤ) + + one(T) * (Tₗ <= T <= Tᵤ)) +end + +# the kelp magically know what day of the year it is!!! +@inline function seasonal_limitation(kelp, t) + φ = kelp.adapted_latitude + a₁ = kelp.photoperiod_1 + a₂ = kelp.photoperiod_2 + + n = floor(Int, mod(t, 364days)/day) + + λ = normed_day_length_change(φ, n) + + return a₁ * (1 + sign(λ) * abs(λ) ^ .5) + a₂ +end + +# hmmmm +normed_day_length_change(φ, n) = (day_length(φ, n) - day_length(φ, n-1)) / (day_length(φ, 76) - day_length(φ, 75)) + +@inline function day_length(φ, n) + n -= 171 + M = mod((356.5291 + 0.98560028 * n), 360) + C = 1.9148 * sin(M * π / 180) + 0.02 * sin(2 * M * π / 180) + 0.0003 * sin(3 * M * π / 180) + λ = mod(M + C + 180 + 102.9372, 360) + δ = asin(sin(λ * π / 180) * sin(23.44 * π / 180)) + ω = (sin(-0.83 * π / 180) * sin(φ * π / 180) * sin(δ)) / (cos(φ * π / 180) * cos(δ)) + return ω / 180 +end \ No newline at end of file diff --git a/src/Models/Individuals/SugarKelp/show.jl b/src/Models/Individuals/SugarKelp/show.jl new file mode 100644 index 000000000..0a8f68b35 --- /dev/null +++ b/src/Models/Individuals/SugarKelp/show.jl @@ -0,0 +1,8 @@ +import Base: summary, show + +summary(::SugarKelp{FT}) where FT = "SugarKelp{FT} biogeochemistry" + +show(io::IO, kelp::SugarKelp) = + print(io, summary(kelp), + " (Broch & Slagstad, 2012) tracking the `N`itrogen and `C`arbon in a frond of `A`rea") + diff --git a/src/Models/Models.jl b/src/Models/Models.jl index 5d51b7aec..20cb4a628 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -23,7 +23,7 @@ include("Sediments/Sediments.jl") include("AdvectedPopulations/LOBSTER/LOBSTER.jl") include("AdvectedPopulations/NPZD.jl") #include("Individuals/SLatissima.jl") -include("Individuals/SugarKelp.jl") +include("Individuals/SugarKelp/SugarKelp.jl") include("seawater_density.jl") include("CarbonChemistry/CarbonChemistry.jl") include("GasExchange/GasExchange.jl") diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 162cf027b..38ff25b1f 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -8,7 +8,7 @@ module OceanBioME export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, NPZD, redfield # Macroalgae models -export SugarKelp +export BiogeochemicalParticles, SugarKelp # Box model export BoxModel, BoxModelGrid, SpeedyOutput, load_output @@ -191,6 +191,7 @@ include("BoxModel/boxmodel.jl") include("Models/Models.jl") using .Light +using .Particles using .BoxModels using .Models diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 333307c1a..12cb1e78a 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -1,5 +1,7 @@ module Particles +export BiogeochemicalParticles + using Adapt using KernelAbstractions: @kernel, @index @@ -37,6 +39,7 @@ function coupled_tracers end @inline required_tracers(p::BiogeochemicalParticles) = required_tracers(p.biogeochemistry) @inline coupled_tracers(p::BiogeochemicalParticles) = coupled_tracers(p.biogeochemistry) +include("atomic_operations.jl") include("advection.jl") include("tracer_interpolation.jl") include("tendencies.jl") @@ -116,17 +119,16 @@ size(particles::BiogeochemicalParticles) = (length(particles), ) length(::BiogeochemicalParticles{N}) where N = N # todo, fix below -Base.summary(particles::BiogeochemicalParticles) = - string(length(particles), " BiogeochemicalParticles with eltype ", nameof(eltype(particles)), - " and properties ", propertynames(particles)) +Base.summary(particles::BiogeochemicalParticles{N}) where N = + string(N, " BiogeochemicalParticles with ", summary(particles.biogeochemistry)) function Base.show(io::IO, particles::BiogeochemicalParticles) - Tparticle = nameof(eltype(particles)) - properties = propertynames(particles) - Nparticles = length(particles) + properties = required_particle_fields(particles) + tracers = coupled_tracers(particles) - print(io, Nparticles, " BiogeochemicalParticles with eltype ", Tparticle, ":", "\n", - "└── ", length(properties), " properties: ", properties, "\n") + print(io, summary(particles), ":", "\n", + "├── fields: ", properties, "\n", + "└── coupled tracers: ", tracers, "\n") end # User may want to overload this to not output parameters over and over again diff --git a/src/Particles/atomic_operations.jl b/src/Particles/atomic_operations.jl new file mode 100644 index 000000000..a520612da --- /dev/null +++ b/src/Particles/atomic_operations.jl @@ -0,0 +1,24 @@ +# this should be in Oceananigans + +using Atomix, CUDA + +using Oceananigans: Field +using OffsetArrays: OffsetArray + +import Atomix: pointer + +@inline function pointer(ref::Atomix.Internal.IndexableRef{<:Field, Tuple{Vararg{Int64, N}}} where {N}) + i = LinearIndices(ref.data.data)[ref.indices...] + return Base.pointer(ref.data.data, i) +end + +# Field on CPU +function atomic_add!(fld::Field, i, j, k, value) + Atomix.@atomic fld[i, j, k] += value +end + +# Field on GPU which is adapted to field.data +function atomic_add!(fld::OffsetArray, i, j, k, value) + linear_index = LinearIndices(fld)[i, j, k] + CUDA.@atomic fld.parent[linear_index] += value +end diff --git a/src/Particles/tendencies.jl b/src/Particles/tendencies.jl index 4c3e2e2be..5e747f0c0 100644 --- a/src/Particles/tendencies.jl +++ b/src/Particles/tendencies.jl @@ -25,9 +25,9 @@ end t = clock.time field_names = required_particle_fields(particles) - nf = length(field_names) + Nf = length(field_names) - field_values = ntuple(n->particles.fields[n], Val(nf)) + field_values = ntuple(nf->particles.fields[nf][n], Val(Nf)) tracer_values = extract_tracer_values(particles.field_interpolation, particles, grid, fields, n) diff --git a/src/Particles/tracer_interpolation.jl b/src/Particles/tracer_interpolation.jl index fe5504c77..aefaab7f9 100644 --- a/src/Particles/tracer_interpolation.jl +++ b/src/Particles/tracer_interpolation.jl @@ -50,7 +50,7 @@ end node_volume = volume(i, j, k, grid, Center(), Center(), Center()) - @inbounds tendency[i, j, k] += particle_tendency / node_volume + atomic_add!(tendency, i, j, k, particle_tendency / node_volume) return nothing end \ No newline at end of file diff --git a/src/Particles/update_tracer_tendencies.jl b/src/Particles/update_tracer_tendencies.jl index 2c9f4eb9f..65df9b225 100644 --- a/src/Particles/update_tracer_tendencies.jl +++ b/src/Particles/update_tracer_tendencies.jl @@ -1,6 +1,4 @@ function update_tendencies!(bgc, particles::BiogeochemicalParticles{N}, model) where N - tendencies = particles.timestepper.tendencies - field_names = coupled_tracers(particles) dev = device(architecture(model)) @@ -10,30 +8,41 @@ function update_tendencies!(bgc, particles::BiogeochemicalParticles{N}, model) w tendency_kernel! = _update_tracer_tendencies!(dev, workgroup, worksize) - for name in field_names - tendency_kernel!(Val(name), particles, model.grid, fields(model), model.clock, model.timestepper.Gⁿ, particles.scalefactors) + for field_name in field_names + target_name = possible_tuple_value(field_names, field_name) # coupled_tracers can return a NamedTuple + + tendency_kernel!(Val(field_name), Val(target_name), particles, model.grid, fields(model), model.clock, model.timestepper.Gⁿ, particles.scalefactors) end return nothing end -@kernel function _update_tracer_tendencies!(val_field_name::Val{field_name}, particles, grid, fields, clock, tendencies, scalefactors) where field_name +possible_tuple_value(field_names, field_name) = field_names[field_name] +possible_tuple_value(::Tuple, field_name) = field_name + +@kernel function _update_tracer_tendencies!(val_field_name::Val{field_name}, ::Val{target_name}, + particles, + grid, + fields, + clock, + tendencies, + scalefactors) where {field_name, target_name} n = @index(Global) t = clock.time field_names = required_particle_fields(particles) - nf = length(field_names) + Nf = length(field_names) - field_values = ntuple(n->particles.fields[n], Val(nf)) + field_values = ntuple(nf->particles.fields[nf][n], Val(Nf)) tracer_values = extract_tracer_values(particles.field_interpolation, particles, grid, fields, n) - + particle_tendency = particles.biogeochemistry(val_field_name, t, field_values..., tracer_values...) scalefactor = @inbounds scalefactors[n] - total_tendency = @show scalefactor * particle_tendency + total_tendency = scalefactor * particle_tendency - apply_tracer_tendency!(particles.field_interpolation, particles, grid, total_tendency, tendencies[field_name], n) + apply_tracer_tendency!(particles.field_interpolation, particles, grid, total_tendency, tendencies[target_name], n) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 66e2c4281..4231c6ac7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,10 @@ include("test_gasexchange_carbon_chem.jl") include("test_sediments.jl") if architecture == CPU() - include("test_boxmodel.jl") # box models (probably) don't work on GPU, and it wouldn't be faster anyway + # box models (probably) don't work on GPU, and it wouldn't be faster anyway + # we would probably want to run over a grid if you were to run an enseble of box models, + # which we would want to do a different way + include("test_boxmodel.jl") end architecture == CPU() && @testset "Doctests" begin diff --git a/test/test_sugar_kelp.jl b/test/test_sugar_kelp.jl new file mode 100644 index 000000000..688542f83 --- /dev/null +++ b/test/test_sugar_kelp.jl @@ -0,0 +1,26 @@ +grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) + +particle_biogeochemistry = OceanBioME.Models.SugarKelpModel.SugarKelp() + +particles = BiogeochemicalParticles(2; grid, + biogeochemistry = particle_biogeochemistry, + advection = nothing) + +biogeochemistry = LOBSTER(; grid, + particles, + carbonates = true, + variable_redfield = true, + oxygen = true, + sinking_speeds = NamedTuple()) + +model = NonhydrostaticModel(; grid, biogeochemistry, advection = nothing, tracers = (:T, :S)) + +set!(model, NO₃ = 10.0, NH₄ = 1.0, DIC = 2000, Alk = 2000, T = 10, S = 35) + +particles.x .= 0.5 +particles.y .= 0.5 +particles.z .= -0.5 + +particles.fields.A .= 2 +particles.fields.N .= 1 +particles.fields.C .= 1 \ No newline at end of file From 3f9fd7e1e7617aa3cfa630ce29e37a0d701abe2a Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 20:45:08 +0100 Subject: [PATCH 04/33] fixed tests, project, and manifest --- Manifest.toml | 2 +- Project.toml | 2 + test/test_sugar_kelp.jl | 101 ++++++++++++++++++++++++++++++++-------- 3 files changed, 85 insertions(+), 20 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 4a090ca54..f35ff00b6 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.2" manifest_format = "2.0" -project_hash = "37e11a3bdd396c973917441db34ded05b97faa85" +project_hash = "5c813c4d993ee6a06601ef0a85ac8b817d5ad5c1" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] diff --git a/Project.toml b/Project.toml index f420f20f9..d592e24c3 100644 --- a/Project.toml +++ b/Project.toml @@ -5,11 +5,13 @@ version = "0.11.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" GibbsSeaWater = "9a22fb26-0b63-4589-b28e-8f9d0b5c3d05" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" +OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" diff --git a/test/test_sugar_kelp.jl b/test/test_sugar_kelp.jl index 688542f83..113ab1516 100644 --- a/test/test_sugar_kelp.jl +++ b/test/test_sugar_kelp.jl @@ -1,26 +1,89 @@ -grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) +using Oceananigans.Architectures: on_architecture -particle_biogeochemistry = OceanBioME.Models.SugarKelpModel.SugarKelp() +sum_tracer_nitrogen(tracers) = sum(on_architecture(CPU(), interior(tracers.NO₃))) + + sum(on_architecture(CPU(), interior(tracers.NH₄))) + + sum(on_architecture(CPU(), interior(tracers.P))) + + sum(on_architecture(CPU(), interior(tracers.Z))) + + sum(on_architecture(CPU(), interior(tracers.sPON))) + + sum(on_architecture(CPU(), interior(tracers.bPON))) + + sum(on_architecture(CPU(), interior(tracers.DON))) -particles = BiogeochemicalParticles(2; grid, - biogeochemistry = particle_biogeochemistry, - advection = nothing) +sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = + sum(on_architecture(CPU(), interior(tracers.sPOC))) + + sum(on_architecture(CPU(), interior(tracers.bPOC))) + + sum(on_architecture(CPU(), interior(tracers.DOC))) + + sum(on_architecture(CPU(), interior(tracers.DIC))) + + sum(on_architecture(CPU(), interior(tracers.P)) .* (1 + organic_carbon_calcate_ratio) .+ on_architecture(CPU(), interior(tracers.Z))) .* redfield -biogeochemistry = LOBSTER(; grid, - particles, - carbonates = true, - variable_redfield = true, - oxygen = true, - sinking_speeds = NamedTuple()) +@testset "SLatissima particle setup and conservations" begin + grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) -model = NonhydrostaticModel(; grid, biogeochemistry, advection = nothing, tracers = (:T, :S)) + particle_biogeochemistry = OceanBioME.Models.SugarKelpModel.SugarKelp() -set!(model, NO₃ = 10.0, NH₄ = 1.0, DIC = 2000, Alk = 2000, T = 10, S = 35) + particles = BiogeochemicalParticles(2; grid, + biogeochemistry = particle_biogeochemistry, + advection = nothing) -particles.x .= 0.5 -particles.y .= 0.5 -particles.z .= -0.5 + biogeochemistry = LOBSTER(; grid, + particles, + carbonates = true, + variable_redfield = true, + oxygen = true, + sinking_speeds = NamedTuple()) -particles.fields.A .= 2 -particles.fields.N .= 1 -particles.fields.C .= 1 \ No newline at end of file + model = NonhydrostaticModel(; grid, biogeochemistry, advection = nothing, tracers = (:T, :S)) + + set!(model, NO₃ = 10.0, NH₄ = 1.0, DIC = 2000, Alk = 2000, T = 10, S = 35) + + particles.x .= 0.5 + particles.y .= 0.5 + particles.z .= -0.5 + + particles.fields.A .= 2 + particles.fields.N .= 1 + particles.fields.C .= 1 + + A = on_architecture(CPU(), particles.fields.A) + N = on_architecture(CPU(), particles.fields.N) + C = on_architecture(CPU(), particles.fields.C) + + Nₛ = particles.biogeochemistry.structural_nitrogen + Cₛ = particles.biogeochemistry.structural_carbon + kₐ = particles.biogeochemistry.structural_dry_weight_per_area + + + initial_tracer_N = sum_tracer_nitrogen(model.tracers) + initial_kelp_N = sum(A .* kₐ .* (N .+ Nₛ)) / (14 * 0.001) + + initial_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield, model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) + initial_kelp_C = sum(A .* kₐ .* (C .+ Cₛ)) / (12 * 0.001) + + model.clock.time = 60days # get to a high growth phase + + for _ in 1:10 + time_step!(model, 1) + end + + # not sure we need to repeat this + A = on_architecture(CPU(), particles.fields.A) + N = on_architecture(CPU(), particles.fields.N) + C = on_architecture(CPU(), particles.fields.C) + + final_tracer_N = sum_tracer_nitrogen(model.tracers) + final_kelp_N = sum(A .* kₐ .* (N .+ Nₛ)) / (14 * 0.001) + + final_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield, model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) + final_kelp_C = sum(A .* kₐ .* (C .+ Cₛ)) / (12 * 0.001) + + # kelp is being integrated + @test initial_kelp_N != final_kelp_N + @test initial_kelp_C != final_kelp_C + + # conservaitons + # (GPU eps is much larger (~10⁻⁷) than on CPU), is this true??? And is this conservation good enough??? + rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_N + initial_kelp_N), √eps(final_tracer_N + final_kelp_N)), 2e-7) + @test isapprox(initial_tracer_N + initial_kelp_N, final_tracer_N + final_kelp_N; rtol) + + rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_C + initial_kelp_C), √eps(final_tracer_C + final_kelp_C)), 7e-7) + @test isapprox(initial_tracer_C + initial_kelp_C, final_tracer_C + final_kelp_C; rtol) +end \ No newline at end of file From 7fb7d783b3aba1f3c0244da42ca3174dc323f348 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 20:50:19 +0100 Subject: [PATCH 05/33] missed an export --- src/Models/Individuals/SugarKelp/SugarKelp.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index b2f1a8bc6..942dddb6a 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -1,5 +1,7 @@ module SugarKelpModel +export SugarKelp + using Roots using Oceananigans.Units From b3a391190d28fe8622c8f2640a609d1973e0690d Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 20:50:49 +0100 Subject: [PATCH 06/33] start making the docs at least --- docs/make.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 4a84eb77b..e3cc9c0a3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,7 @@ using Documenter, DocumenterCitations, Literate using OceanBioME -using OceanBioME: SLatissima, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus +using OceanBioME: SugarKelp, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation using OceanBioME: CarbonChemistry, GasExchange @@ -53,7 +53,7 @@ if !isdir(OUTPUT_DIR) mkdir(OUTPUT_DIR) end model_parameters = (LOBSTER(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry, - SLatissima(), + SugarKelp().biogeochemistry, TwoBandPhotosyntheticallyActiveRadiation(; grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))), SimpleMultiG(; grid = BoxModelGrid()), InstantRemineralisation(; grid = BoxModelGrid()), From 12b513cb546e8080fbe3574f884e4380298fce3e Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 22:55:56 +0100 Subject: [PATCH 07/33] fixed some bugs and updated example --- examples/kelp.jl | 61 +- src/Models/Individuals/SLatissima.jl | 570 ------------------ src/Models/Individuals/SugarKelp/SugarKelp.jl | 6 + src/Models/Individuals/SugarKelp/coupling.jl | 8 +- src/Particles/Particles.jl | 9 +- src/Particles/set.jl | 24 + test/test_particles.jl | 4 + test/test_slatissima.jl | 98 --- test/test_sugar_kelp.jl | 2 +- 9 files changed, 74 insertions(+), 708 deletions(-) delete mode 100644 src/Models/Individuals/SLatissima.jl create mode 100644 src/Particles/set.jl delete mode 100644 test/test_slatissima.jl diff --git a/examples/kelp.jl b/examples/kelp.jl index 2d1921824..efef43f9e 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -57,22 +57,21 @@ S = ConstantField(35) # ## Kelp Particle setup @info "Setting up kelp particles" -n = 5 # number of kelp fronds +n = 5 # number of kelp bundles z₀ = [-21:5:-1;] * 1.0 # depth of kelp fronds -particles = SLatissima(; architecture, - x = on_architecture(architecture, ones(n) * Lx / 2), - y = on_architecture(architecture, ones(n) * Ly / 2), - z = on_architecture(architecture, z₀), - A = on_architecture(architecture, ones(n) * 10.0), - latitude = 57.5, - scalefactor = 500.0) +particles = SugarKelp(n; grid, + advection = nothing, # we don't want them to move around + scalefactors = fill(2000, n)) # and we want them to look like there are 500 in each bundle + +set!(particles, A = 10, N = 0.01, C = 0.1, z = z₀, x = Lx / 2, y = Ly / 2) # ## Setup BGC model biogeochemistry = LOBSTER(; grid, surface_photosynthetically_active_radiation = PAR⁰, carbonates = true, variable_redfield = true, + oxygen = true, scale_negatives = true, particles) @@ -90,7 +89,7 @@ set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2 # - Store the model and particles output # - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implementation page) -simulation = Simulation(model, Δt = 3minutes, stop_time = 100days) +simulation = Simulation(model, Δt = 4minutes, stop_time = 150days) progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", iteration(sim), @@ -119,16 +118,10 @@ run!(simulation) # Now we can visualise the results with some post processing to diagnose the air-sea CO₂ flux - hopefully this looks different to the example without kelp! - P = FieldTimeSeries("$filename.jld2", "P") - NO₃ = FieldTimeSeries("$filename.jld2", "NO₃") - Z = FieldTimeSeries("$filename.jld2", "Z") -sPOC = FieldTimeSeries("$filename.jld2", "sPOC") -bPOC = FieldTimeSeries("$filename.jld2", "bPOC") - DIC = FieldTimeSeries("$filename.jld2", "DIC") - Alk = FieldTimeSeries("$filename.jld2", "Alk") +tracers = FieldDataset("$filename.jld2") -x, y, z = nodes(P) -times = P.times +x, y, z = nodes(tracers["P"]) +times = tracers["P"].times nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and @@ -142,34 +135,34 @@ using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity for (n, t) in enumerate(times) clock.time = t - air_sea_CO₂_flux[n] = CO₂_flux.condition.func(1, 1, grid, clock, (; DIC = DIC[n], Alk = Alk[n], T, S)) + air_sea_CO₂_flux[n] = CO₂_flux.condition.func(1, 1, grid, clock, (; DIC = tracers["DIC"][n], Alk = tracers["Alk"][n], T, S)) - carbon_export[n] = sPOC[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + - bPOC[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOC)).w[1, 1, grid.Nz-20] + carbon_export[n] = tracers["sPOC"][n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + + tracers["bPOC"][n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOC)).w[1, 1, grid.Nz-20] end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. using CairoMakie -fig = Figure(size = (1000, 1500), fontsize = 20) +fig = Figure(size = (1000, 1200), fontsize = 20) axis_kwargs = (xlabel = "Time (days)", ylabel = "z (m)", limits = ((0, times[end] / days), (-85meters, 0))) axP = Axis(fig[1, 1]; title = "Phytoplankton concentration (mmol N/m³)", axis_kwargs...) -hmP = heatmap!(times / days, z, interior(P, 1, 1, :, :)', colormap = :batlow) +hmP = heatmap!(times / days, z, interior(tracers["P"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[1, 2], hmP) axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N/m³)", axis_kwargs...) -hmNO₃ = heatmap!(times / days, z, interior(NO₃, 1, 1, :, :)', colormap = :batlow) +hmNO₃ = heatmap!(times / days, z, interior(tracers["NO₃"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[2, 2], hmNO₃) axZ = Axis(fig[3, 1]; title = "Zooplankton concentration (mmol N/m³)", axis_kwargs...) -hmZ = heatmap!(times / days, z, interior(Z, 1, 1, :, :)', colormap = :batlow) +hmZ = heatmap!(times / days, z, interior(tracers["Z"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[3, 2], hmZ) axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol C/m³)", axis_kwargs...) -hmD = heatmap!(times / days, z, interior(sPOC, 1, 1, :, :)' .+ interior(bPOC, 1, 1, :, :)', colormap = :batlow) +hmD = heatmap!(times / days, z, interior(tracers["sPOC"], 1, 1, :, :)' .+ interior(tracers["bPOC"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[4, 2], hmD) CO₂_molar_mass = (12 + 2 * 16) * 1e-3 # kg / mol @@ -194,14 +187,18 @@ A, N, C = ntuple(n -> zeros(5, length(iterations)), 3) times = zeros(length(iterations)) for (i, iter) in enumerate(iterations) - particles = file["timeseries/particles/$iter"] - A[:, i] = particles.A - N[:, i] = particles.N - C[:, i] = particles.C + particles_values = file["timeseries/particles/$iter"] + A[:, i] = particles_values.A + N[:, i] = particles_values.N + C[:, i] = particles_values.C times[i] = file["timeseries/t/$iter"] end +Nₛ = particles.biogeochemistry.structural_nitrogen +Cₛ = particles.biogeochemistry.structural_carbon +kₐ = particles.biogeochemistry.structural_dry_weight_per_area + fig = Figure(size = (1000, 800), fontsize = 20) axis_kwargs = (xlabel = "Time (days)", limits = ((0, times[end] / days), nothing)) @@ -210,9 +207,9 @@ ax1 = Axis(fig[1, 1]; ylabel = "Frond area (dm²)", axis_kwargs...) [lines!(ax1, times / day, A[n, :], linewidth = 3) for n in 1:5] ax2 = Axis(fig[2, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) -[lines!(ax2, times / day, (@. A * (N + particles.structural_nitrogen) * particles.structural_dry_weight_per_area)[n, :], linewidth = 3) for n in 1:5] +[lines!(ax2, times / day, (@. A * (N + Nₛ) * kₐ)[n, :], linewidth = 3) for n in 1:5] ax3 = Axis(fig[3, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) -[lines!(ax3, times / day, (@. A * (C + particles.structural_carbon) * particles.structural_dry_weight_per_area)[n, :], linewidth = 3) for n in 1:5] +[lines!(ax3, times / day, (@. A * (C + Cₛ) * kₐ)[n, :], linewidth = 3) for n in 1:5] fig diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl deleted file mode 100644 index 72f389ab1..000000000 --- a/src/Models/Individuals/SLatissima.jl +++ /dev/null @@ -1,570 +0,0 @@ -""" -Sugar kelp model of [Broch2012](@citet) and updated by [Broch2013](@citet), [Fossberg2018](@citet), and [Broch2019](@citet). - -Prognostic properties -===================== -* Area: A (dm²) -* Nitrogen reserve: N (gN/gSW) -* Carbon reserve: C (gC/gSW) - -Tracer dependencies -=================== -* Nitrates: NO₃ (mmol N/m³) -* Photosynthetically available radiation: PAR (einstein/m²/day) - -Optional -======== -* Ammonia: NH₄ (mmol N/m³) -""" -module SLatissimaModel - -export SLatissima - -using Roots, KernelAbstractions -using OceanBioME.Particles: BiogeochemicalParticles, get_node -using Oceananigans.Units -using Oceananigans: Center, CPU -using Oceananigans.Architectures: on_architecture, device, architecture -using Oceananigans.Biogeochemistry: required_biogeochemical_tracers, biogeochemical_auxiliary_fields -using Oceananigans.Models.LagrangianParticleTracking: _advect_particles! - -using Oceananigans.Operators: volume -using Oceananigans.Grids: AbstractGrid, Flat -using Oceananigans.Fields: fractional_indices, _interpolate, interpolator -using Oceananigans.Models: total_velocities - -import Adapt: adapt_structure, adapt -import Oceananigans.Biogeochemistry: update_tendencies! -import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties! - -@inline no_dynamics(args...) = nothing - - -""" - SLatissima(; architecture :: AR = CPU(), - growth_rate_adjustment :: FT = 4.5, - photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60), - minimum_carbon_reserve :: FT = 0.01, - structural_carbon :: FT = 0.2, - exudation :: FT = 0.5, - erosion :: FT = 0.22, - saturation_irradiance :: FT = 90 * day/ (10 ^ 6), - structural_dry_weight_per_area :: FT = 0.5, - structural_dry_to_wet_weight :: FT = 0.0785, - carbon_reserve_per_carbon :: FT = 2.1213, - nitrogen_reserve_per_nitrogen :: FT = 2.72, - minimum_nitrogen_reserve :: FT = 0.0126, - maximum_nitrogen_reserve :: FT = 0.0216, - growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)), - growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2, - maximum_specific_growth_rate :: FT = 0.18, - structural_nitrogen :: FT = 0.0146, - photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24, - photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24, - photosynthesis_ref_temp_1 :: FT = 285.0, - photosynthesis_ref_temp_2 :: FT = 288.0, - photoperiod_1 :: FT = 0.85, - photoperiod_2 :: FT = 0.3, - respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24, - respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24, - respiration_ref_temp_1 :: FT = 285.0, - respiration_ref_temp_2 :: FT = 290.0, - photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1), - photosynthesis_low_temp :: FT = 271.0, - photosynthesis_high_temp :: FT = 296.0, - photosynthesis_high_arrhenius_temp :: FT = 1414.87, - photosynthesis_low_arrhenius_temp :: FT = 4547.89, - respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1), - current_speed_for_0p65_uptake :: FT = 0.03, - nitrate_half_saturation :: FT = 4.0, - ammonia_half_saturation :: FT = 1.3, - maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6), - maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6), - current_1 :: FT = 0.72, - current_2 :: FT = 0.28, - current_3 :: FT = 0.045, - respiration_reference_A :: FT = 1.11e-4 * 24, - respiration_reference_B :: FT = 5.57e-5 * 24, - exudation_redfield_ratio :: FT = Inf, - - prescribed_velocity :: U = 0.1, - - #position - x :: P = on_architecture(architecture, [0.0]) - y :: P = on_architecture(architecture, zeros(Float64, length(x))), - z :: P = on_architecture(architecture, zeros(Float64, length(x))), - - #properties - A :: P = on_architecture(architecture, ones(Float64, length(x)) * 30), - N :: P = on_architecture(architecture, ones(Float64, length(x)) * 0.01), - C :: P = on_architecture(architecture, ones(Float64, length(x)) * 0.1), - - #feedback - nitrate_uptake :: P = on_architecture(architecture, zeros(Float64, length(x))), - ammonia_uptake :: P = on_architecture(architecture, zeros(Float64, length(x))), - primary_production :: P = on_architecture(architecture, zeros(Float64, length(x))), - frond_exudation :: P = on_architecture(architecture, zeros(Float64, length(x))), - nitrogen_erosion :: P = on_architecture(architecture, zeros(Float64, length(x))), - carbon_erosion :: P = on_architecture(architecture, zeros(Float64, length(x))), - - custom_dynamics :: F = no_dynamics, - - scalefactor :: FT = 1.0, - latitude :: FT = 57.5) - -Keyword Arguments -================= - -- `architecture`: the architecture to adapt arrays to -- `growth_rate_adjustment`, ..., `exudation_redfield_ratio`: parameter values -- `prescribed_velocity`: functions for the relative velocity `f(x, y, z, t)` -- `x`,`y` and `z`: positions of the particles -- `A`, `N`, and `C`: area, nitrogen, and carbon reserves -- `nitrate_uptake` ... `carbon_erosion`: diagnostic values coupled to tracer fields -- `custom_dynamics`: place to add any function of form `f!(particles, model, bgc, Δt)` -- `scalefactor`: scalar scaling for tracer coupling -- `latitude`: model latitude for seasonal growth modulation -""" -Base.@kwdef struct SLatissima{AR, FT, U, P, F} <: BiogeochemicalParticles - architecture :: AR = CPU() - growth_rate_adjustment :: FT = 4.5 - photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60) - minimum_carbon_reserve :: FT = 0.01 - structural_carbon :: FT = 0.2 - exudation :: FT = 0.5 - erosion :: FT = 0.22 - saturation_irradiance :: FT = 90 * day/ (10 ^ 6) - structural_dry_weight_per_area :: FT = 0.5 - structural_dry_to_wet_weight :: FT = 0.0785 - carbon_reserve_per_carbon :: FT = 2.1213 - nitrogen_reserve_per_nitrogen :: FT = 2.72 - minimum_nitrogen_reserve :: FT = 0.0126 - maximum_nitrogen_reserve :: FT = 0.0216 - growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2 - maximum_specific_growth_rate :: FT = 0.18 - structural_nitrogen :: FT = 0.0146 - photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24 - photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24 - photosynthesis_ref_temp_1 :: FT = 285.0 - photosynthesis_ref_temp_2 :: FT = 288.0 - photoperiod_1 :: FT = 0.85 - photoperiod_2 :: FT = 0.3 - respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24 - respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24 - respiration_ref_temp_1 :: FT = 285.0 - respiration_ref_temp_2 :: FT = 290.0 - photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1) - photosynthesis_low_temp :: FT = 271.0 - photosynthesis_high_temp :: FT = 296.0 - photosynthesis_high_arrhenius_temp :: FT = 1414.87 - photosynthesis_low_arrhenius_temp :: FT = 4547.89 - respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1) - current_speed_for_0p65_uptake :: FT = 0.03 - nitrate_half_saturation :: FT = 4.0 - ammonia_half_saturation :: FT = 1.3 - maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6) - maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6) - current_1 :: FT = 0.72 - current_2 :: FT = 0.28 - current_3 :: FT = 0.045 - respiration_reference_A :: FT = 1.11e-4 * 24 - respiration_reference_B :: FT = 5.57e-5 * 24 - exudation_redfield_ratio :: FT = Inf - - prescribed_velocity :: U = 0.1 - - #position - x :: P = on_architecture(architecture, [0.0]) - y :: P = on_architecture(architecture, zeros(Float64, length(x))) - z :: P = on_architecture(architecture, zeros(Float64, length(x))) - - #properties - A :: P = on_architecture(architecture, ones(Float64, length(x)) * 30) - N :: P = on_architecture(architecture, ones(Float64, length(x)) * 0.01) - C :: P = on_architecture(architecture, ones(Float64, length(x)) * 0.1) - - #feedback - nitrate_uptake :: P = on_architecture(architecture, zeros(Float64, length(x))) - ammonia_uptake :: P = on_architecture(architecture, zeros(Float64, length(x))) - primary_production :: P = on_architecture(architecture, zeros(Float64, length(x))) - frond_exudation :: P = on_architecture(architecture, zeros(Float64, length(x))) - nitrogen_erosion :: P = on_architecture(architecture, zeros(Float64, length(x))) - carbon_erosion :: P = on_architecture(architecture, zeros(Float64, length(x))) - - custom_dynamics :: F = no_dynamics - - scalefactor :: FT = 1.0 - latitude :: FT = 57.5 -end - -adapt_structure(to, kelp::SLatissima) = SLatissima(adapt(to, kelp.architecture), - adapt(to, kelp.growth_rate_adjustment), - adapt(to, kelp.photosynthetic_efficiency), - adapt(to, kelp.minimum_carbon_reserve), - adapt(to, kelp.structural_carbon), - adapt(to, kelp.exudation), - adapt(to, kelp.erosion), - adapt(to, kelp.saturation_irradiance), - adapt(to, kelp.structural_dry_weight_per_area), - adapt(to, kelp.structural_dry_to_wet_weight), - adapt(to, kelp.carbon_reserve_per_carbon), - adapt(to, kelp.nitrogen_reserve_per_nitrogen), - adapt(to, kelp.minimum_nitrogen_reserve), - adapt(to, kelp.maximum_nitrogen_reserve), - adapt(to, kelp.growth_adjustment_2), - adapt(to, kelp.growth_adjustment_1), - adapt(to, kelp.maximum_specific_growth_rate), - adapt(to, kelp.structural_nitrogen), - adapt(to, kelp.photosynthesis_at_ref_temp_1), - adapt(to, kelp.photosynthesis_at_ref_temp_2), - adapt(to, kelp.photosynthesis_ref_temp_1), - adapt(to, kelp.photosynthesis_ref_temp_2), - adapt(to, kelp.photoperiod_1), - adapt(to, kelp.photoperiod_2), - adapt(to, kelp.respiration_at_ref_temp_1), - adapt(to, kelp.respiration_at_ref_temp_2), - adapt(to, kelp.respiration_ref_temp_1), - adapt(to, kelp.respiration_ref_temp_2), - adapt(to, kelp.photosynthesis_arrhenius_temp), - adapt(to, kelp.photosynthesis_low_temp), - adapt(to, kelp.photosynthesis_high_temp), - adapt(to, kelp.photosynthesis_high_arrhenius_temp), - adapt(to, kelp.photosynthesis_low_arrhenius_temp), - adapt(to, kelp.respiration_arrhenius_temp), - adapt(to, kelp.current_speed_for_0p65_uptake), - adapt(to, kelp.nitrate_half_saturation), - adapt(to, kelp.ammonia_half_saturation), - adapt(to, kelp.maximum_nitrate_uptake), - adapt(to, kelp.maximum_ammonia_uptake), - adapt(to, kelp.current_1), - adapt(to, kelp.current_2), - adapt(to, kelp.current_3), - adapt(to, kelp.respiration_reference_A), - adapt(to, kelp.respiration_reference_B), - adapt(to, kelp.exudation_redfield_ratio), - adapt(to, kelp.prescribed_velocity), - adapt(to, kelp.x), - adapt(to, kelp.y), - adapt(to, kelp.z), - adapt(to, kelp.A), - adapt(to, kelp.N), - adapt(to, kelp.C), - adapt(to, kelp.nitrate_uptake), - adapt(to, kelp.ammonia_uptake), - adapt(to, kelp.primary_production), - adapt(to, kelp.frond_exudation), - adapt(to, kelp.nitrogen_erosion), - adapt(to, kelp.carbon_erosion), - adapt(to, kelp.custom_dynamics), - adapt(to, kelp.scalefactor), - adapt(to, kelp.latitude)) - -function update_tendencies!(bgc, particles::SLatissima, model) - num_particles = length(particles) - workgroup = min(num_particles, 256) - worksize = num_particles - - update_tracer_tendencies_kernel! = update_tracer_tendencies!(device(model.architecture), workgroup, worksize) - update_tracer_tendencies_kernel!(bgc, particles, model.timestepper.Gⁿ, model.grid) - - return nothing -end - -@kernel function update_tracer_tendencies!(bgc, p, tendencies, grid::AbstractGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} - idx = @index(Global) - - x = p.x[idx] - y = p.y[idx] - z = p.z[idx] - - bgc_tracers = required_biogeochemical_tracers(bgc) - - ii, jj, kk = fractional_indices((x, y, z), grid, ifelse(isa(TX(), Flat), nothing, Center()), - ifelse(isa(TY(), Flat), nothing, Center()), - ifelse(isa(TZ(), Flat), nothing, Center())) - - ix = interpolator(ii) - iy = interpolator(jj) - iz = interpolator(kk) - - i, j, k = (get_node(TX(), Int(ifelse(ix[3] < 0.5, ix[1], ix[2])), grid.Nx), - get_node(TY(), Int(ifelse(iy[3] < 0.5, iy[1], iy[2])), grid.Ny), - get_node(TZ(), Int(ifelse(iz[3] < 0.5, iz[1], iz[2])), grid.Nz)) - - node_volume = volume(i, j, k, grid, Center(), Center(), Center()) - - node_scalefactor = p.scalefactor / node_volume #* normfactor / (weight * node_volume) - - @inbounds begin - tendencies.NO₃[i, j, k] -= node_scalefactor * p.nitrate_uptake[idx] - - if :NH₄ in bgc_tracers - tendencies.NH₄[i, j, k] -= node_scalefactor * p.ammonia_uptake[idx] - end - - if :DIC in bgc_tracers - tendencies.DIC[i, j, k] -= node_scalefactor * p.primary_production[idx] - end - - if :O₂ in bgc_tracers - tendencies.O₂[i, j, k] += node_scalefactor * p.primary_production[idx] - end - - if :DOM in bgc_tracers - tendencies.DOM[i, j, k] += node_scalefactor * p.frond_exudation[idx] / p.exudation_redfield_ratio - elseif :DON in bgc_tracers - tendencies.DON[i, j, k] += node_scalefactor * p.frond_exudation[idx] / p.exudation_redfield_ratio - tendencies.DOC[i, j, k] += node_scalefactor * p.frond_exudation[idx] - end - - if :bPOM in bgc_tracers - tendencies.bPOM[i, j, k] += node_scalefactor * p.nitrogen_erosion[idx] - elseif :bPON in bgc_tracers - tendencies.bPON[i, j, k] += node_scalefactor * p.nitrogen_erosion[idx] - tendencies.bPOC[i, j, k] += node_scalefactor * p.carbon_erosion[idx] - end - end -end - -function update_lagrangian_particle_properties!(particles::SLatissima, model, bgc, Δt) - workgroup = min(length(particles), 256) - worksize = length(particles) - - arch = architecture(model) - - # Advect particles - advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) - - advect_particles_kernel!(particles, 1.0, model.grid, Δt, total_velocities(model)) - - update_particle_properties_kernel! = _update_lagrangian_particle_properties!(device(arch), workgroup, worksize) - - tracer_fields = merge(model.tracers, model.auxiliary_fields) - - update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying_biogeochemistry, model.grid, - total_velocities(model), tracer_fields, model.clock, Δt) - - particles.custom_dynamics(particles, model, bgc, Δt) -end - -@kernel function _update_lagrangian_particle_properties!(p, light_attenuation, bgc, grid, velocities, tracers, clock, Δt) - idx = @index(Global) - - @inbounds begin - x = p.x[idx] - y = p.y[idx] - z = p.z[idx] - - A = p.A[idx] - N = p.N[idx] - C = p.C[idx] - end - - t = clock.time - - @inbounds if p.A[idx] > 0 - NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(light_attenuation).PAR) - - photo = photosynthesis(T, PAR, p) - e = exudation(C, p) - ν = erosion(A, p) - - fᶜ = f_curr(u, p) - - j_NO₃ = max(0, p.maximum_ammonia_uptake * fᶜ * ((p.maximum_nitrogen_reserve - N) / - (p.maximum_nitrogen_reserve - p.minimum_nitrogen_reserve)) * NO₃ / - (p.nitrate_half_saturation + NO₃)) - - j̃_NH₄ = max(0, p.maximum_ammonia_uptake * fᶜ * NH₄ / (p.ammonia_half_saturation + NH₄)) - - μ_NH₄ = j̃_NH₄ / (p.structural_dry_weight_per_area * (N + p.structural_nitrogen)) - μ_NO₃ = 1 - p.minimum_nitrogen_reserve / N - μ_C = 1 - p.minimum_carbon_reserve / C - - n = floor(Int, mod(t, 364days) / day) - λ = normed_day_length_change(n, p.latitude) - - μ = @inbounds f_area(A, p) * - f_temp(T, p) * - f_photo(λ, p) * - min(μ_C, max(μ_NO₃, μ_NH₄)) - - j_NH₄ = min(j̃_NH₄, μ * p.structural_dry_weight_per_area * (N + p.structural_nitrogen)) - - r = respiration(T, μ, j_NO₃ + j_NH₄, p) - - dA = (μ - ν) * A / day - - dN = ((j_NO₃ + j_NH₄ - photo * e * 14 / (12 * p.exudation_redfield_ratio)) / p.structural_dry_weight_per_area - - μ * (N + p.structural_nitrogen)) / day - - dC = ((photo * (1 - e) - r) / p.structural_dry_weight_per_area - - μ * (C + p.structural_carbon)) / day - - A_new = A + dA * Δt - N_new = N + dN * Δt - C_new = C + dC * Δt - - if C_new < p.minimum_carbon_reserve - A_new *= (1 - (p.minimum_carbon_reserve - C) / p.structural_carbon) - C_new = p.minimum_carbon_reserve - N_new += p.structural_nitrogen * (p.minimum_carbon_reserve - C) / p.structural_carbon - end - - if N_new < p.minimum_nitrogen_reserve - A_new *= (1 - (p.minimum_nitrogen_reserve - N) / p.structural_nitrogen) - N_new = p.minimum_nitrogen_reserve - C_new += p.structural_carbon * (p.minimum_nitrogen_reserve - N) / p.structural_nitrogen - end - - @inbounds begin - p.primary_production[idx] = (photo - r) * A / (day * 12 * 0.001) #gC/dm^2/hr to mmol C/s - - p.frond_exudation[idx] = e * photo * A / (day * 12 * 0.001)#mmol C/s - - p.nitrogen_erosion[idx] = ν * p.structural_dry_weight_per_area * A * (p.structural_nitrogen + N) / (day * 14 * 0.001)#1/day to mmol N/s - - p.carbon_erosion[idx] = ν * p.structural_dry_weight_per_area * A * (p.structural_carbon + C) / (day * 12 * 0.001)#1/hr to mmol C/s - - p.nitrate_uptake[idx] = j_NO₃ * A / (day * 14 * 0.001)#gN/dm^2/hr to mmol N/s - - p.ammonia_uptake[idx] = j_NH₄ * A / (day * 14 * 0.001)#gN/dm^2/hr to mmol N/s - - p.A[idx] = A_new - p.N[idx] = N_new - p.C[idx] = C_new - end - end -end - -@inline function photosynthesis(T, PAR, p) - Tk = T + 273.15 - - pₘₐₓ = p.photosynthesis_at_ref_temp_1 * - exp(p.photosynthesis_arrhenius_temp / p.photosynthesis_ref_temp_1 - - p.photosynthesis_arrhenius_temp / Tk) / - (1 + exp(p.photosynthesis_low_arrhenius_temp / Tk - - p.photosynthesis_low_arrhenius_temp / p.photosynthesis_low_temp) - + exp(p.photosynthesis_high_arrhenius_temp / p.photosynthesis_high_temp - - p.photosynthesis_high_arrhenius_temp / Tk)) - - β = find_zero(β_func, (0, 0.1), Bisection(); p = (; pₘₐₓ, α = p.photosynthetic_efficiency, I_sat = p.saturation_irradiance)) - - pₛ = p.photosynthetic_efficiency * p.saturation_irradiance / log(1 + p.photosynthetic_efficiency / β) - - return pₛ * (1 - exp(- p.photosynthetic_efficiency * PAR / pₛ)) * exp(-β * PAR / pₛ) -end - -@inline β_func(x, p) = p.pₘₐₓ - (p.α * p.I_sat / log(1 + p.α / x)) * - (p.α / (p.α + x)) * (x / (p.α + x))^(x / p.α) - -@inline exudation(C, p) = 1 - exp(p.exudation * (p.minimum_carbon_reserve - C)) - -@inline erosion(A, p) = 1e-6 * exp(p.erosion * A) / (1 + 1e-6 * (exp(p.erosion * A) - 1)) - -@inline respiration(T, μ, j, p) = (p.respiration_reference_A * (μ / p.maximum_specific_growth_rate + - j / (p.maximum_nitrate_uptake + p.maximum_ammonia_uptake)) + - p.respiration_reference_B) * - exp(p.respiration_arrhenius_temp / p.respiration_ref_temp_1 - p.respiration_arrhenius_temp / (T + 273.15)) - -##### -##### Growth limitation -##### - -@inline f_curr(u, p) = p.current_1 * (1 - exp(-u / p.current_3)) + p.current_2 - -@inline f_area(a, p) = p.growth_adjustment_1 * exp(-(a / p.growth_rate_adjustment)^2) + p.growth_adjustment_2 - -@inline function f_temp(T, p) - # should probably parameterise these limits - if -1.8 <= T < 10 - return 0.08 * T + 0.2 - elseif 10 <= T <= 15 - return 1 - elseif 15 < T <= 19 - return 19 / 4 - T / 4 - elseif T > 19 - return 0 - else - return 0 - end -end - -@inline f_photo(λ, p) = p.photoperiod_1 * (1 + sign(λ) * abs(λ) ^ 0.5) + p.photoperiod_2 - -@inline function day_length(n, ϕ) - n -= 171 - M = mod((356.5291 + 0.98560028 * n), 360) - C = 1.9148 * sin(M * π / 180) + 0.02 * sin(2 * M * π / 180) + 0.0003 * sin(3 * M * π / 180) - λ = mod(M + C + 180 + 102.9372, 360) - δ = asin(sin(λ * π / 180) * sin(23.44 * π / 180)) - ω = (sin(-0.83 * π / 180) * sin(ϕ * π / 180) * sin(δ)) / (cos(ϕ * π / 180) * cos(δ)) - return ω / 180 -end - -@inline normed_day_length_change(n, ϕ) = (day_length(n, ϕ) - day_length(n - 1, ϕ)) / (day_length(76, ϕ) - day_length(75, ϕ)) - -@inline function get_arguments(x, y, z, t, particles, bgc, grid::AbstractGrid{FT, TX, TY, TZ}, velocities, tracers, PAR_field) where {FT, TX, TY, TZ} - bgc_tracers = required_biogeochemical_tracers(bgc) - - - ii, jj, kk = fractional_indices((x, y, z), grid, ifelse(isa(TX(), Flat), nothing, Center()), - ifelse(isa(TY(), Flat), nothing, Center()), - ifelse(isa(TZ(), Flat), nothing, Center())) - - ix = interpolator(ii) - iy = interpolator(jj) - iz = interpolator(kk) - - # TODO: ADD ALIASING/RENAMING OF TRACERS SO WE CAN USE E.G. N IN STEAD OF NO3 - - NO₃ = _interpolate(tracers.NO₃, ix, iy, iz) - PAR = _interpolate(PAR_field, ix, iy, iz) * day / (3.99e-10 * 545e12) # W / m² / s to einstein / m² / day - - if :NH₄ in bgc_tracers - NH₄ = _interpolate(tracers.NH₄, ix, iy, iz) - else - NH₄ = 0.0 - end - - T = _interpolate(tracers.T, ix, iy, iz) - - S = _interpolate(tracers.S, ix, iy, iz) - - if isnothing(particles.prescribed_velocity) - ii, jj, kk = fractional_indices((x, y, z), grid, Face(), Center(), Center()) - - ix = interpolator(ii) - iy = interpolator(jj) - iz = interpolator(kk) - - u = _interpolate(velocities.u, ix, iy, iz) - - ii, jj, kk = fractional_indices((x, y, z), grid, Center(), Face(), Center()) - - ix = interpolator(ii) - iy = interpolator(jj) - iz = interpolator(kk) - - v = _interpolate(velocities.v, ix, iy, iz) - - ii, jj, kk = fractional_indices((x, y, z), grid, Center(), Center(), Face()) - - ix = interpolator(ii) - iy = interpolator(jj) - iz = interpolator(kk) - - w = _interpolate(velocities.w, ix, iy, iz) - - speed = sqrt(u^2 + v^2 + w^2) - elseif isa(particles.prescribed_velocity, Number) - speed = particles.prescribed_velocity - else - speed = particles.prescribed_velocity(x, y, z, t) - end - - return NO₃, NH₄, PAR, speed, T, S -end - -end #module diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index 942dddb6a..01c1f95c6 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -6,6 +6,8 @@ using Roots using Oceananigans.Units +using OceanBioME.Particles: BiogeochemicalParticles + import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers # disgsuting number of parameters @@ -59,6 +61,10 @@ import OceanBioME.Particles: required_particle_fields, required_tracers, coupled adapted_latitude :: FT = 57.5 end +# convenience constructor that sets up with default parameters +SugarKelp(number; grid, kelp_parameters = NamedTuple(), kwargs...) = + BiogeochemicalParticles(number; grid, biogeochemistry = SugarKelp(; kelp_parameters...), kwargs...) + @inline required_particle_fields(::SugarKelp) = (:A, :N, :C) @inline required_tracers(::SugarKelp) = (:u, :v, :w, :T, :NO₃, :NH₄, :PAR) @inline coupled_tracers(::SugarKelp) = (:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON) diff --git a/src/Models/Individuals/SugarKelp/coupling.jl b/src/Models/Individuals/SugarKelp/coupling.jl index 9794c22e5..d25220048 100644 --- a/src/Models/Individuals/SugarKelp/coupling.jl +++ b/src/Models/Individuals/SugarKelp/coupling.jl @@ -3,13 +3,13 @@ @inline function (kelp::SugarKelp)(::Val{:NO₃}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) J = nitrate_uptake(kelp, N, NO₃, u, v, w) - return J * A / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s + return -J * A / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s end @inline function (kelp::SugarKelp)(::Val{:NH₄}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) J = ammonia_uptake(kelp, t, A, N, C, T, NH₄, u, v, w) - return J * A / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s + return -J * A / (day * 14 * 0.001) # gN/dm^2/hr to mmol N/s end @inline function (kelp::SugarKelp)(::Val{:DIC}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) @@ -19,7 +19,7 @@ end R = respiration(kelp, t, A, N, C, T, NO₃, NH₄, u, v, w, μ) - return (P - R) * A / (day * 12 * 0.001) # gC/dm^2/hr to mmol C/s + return -(P - R) * A / (day * 12 * 0.001) # gC/dm^2/hr to mmol C/s end # I now know that this may not be correct..., it probably should vary by where the nitrogen comes from @@ -53,5 +53,5 @@ end ν = erosion(kelp, t, A, N, C, T) - return ν * kₐ * A * (C + Cₛ) / (day * 14 * 0.001) # gC/dm^2/hr to mmol C/s + return ν * kₐ * A * (C + Cₛ) / (day * 12 * 0.001) # gC/dm^2/hr to mmol C/s end diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 12cb1e78a..5bb3ff7d0 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -12,6 +12,7 @@ using OceanBioME: Biogeochemistry using Oceananigans.Architectures: architecture, on_architecture +import Oceananigans.Architectures: architecture import Oceananigans.Biogeochemistry: update_tendencies! import Oceananigans.Fields: set! import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, step_lagrangian_particles! @@ -31,6 +32,8 @@ struct BiogeochemicalParticles{N, B, A, F, T, I, S, X, Y, Z} z :: Z end +architecture(p::BiogeochemicalParticles) = architecture(p.x) + function required_particle_fields end function required_tracers end function coupled_tracers end @@ -45,6 +48,7 @@ include("tracer_interpolation.jl") include("tendencies.jl") include("time_stepping.jl") include("update_tracer_tendencies.jl") +include("set.jl") function BiogeochemicalParticles(number; grid, @@ -133,9 +137,8 @@ end # User may want to overload this to not output parameters over and over again function fetch_output(particles::BiogeochemicalParticles, model) - names = propertynames(particles) - return NamedTuple{names}([getproperty(particles, name) for name in names]) + names = propertynames(particles.fields) + return NamedTuple{names}([getproperty(particles.fields, name) for name in names]) end -# TODO: implement set! end #module diff --git a/src/Particles/set.jl b/src/Particles/set.jl new file mode 100644 index 000000000..ea425408e --- /dev/null +++ b/src/Particles/set.jl @@ -0,0 +1,24 @@ +using Oceananigans.Architectures: architecture, on_architecture + +import Oceananigans.Fields: set! + +function set!(particles::BiogeochemicalParticles; kwargs...) + arch = architecture(particles) + + for (n, v) in pairs(kwargs) + if n == :x + particles.x .= on_architecture(arch, v) + elseif n == :y + particles.y .= on_architecture(arch, v) + elseif n == :z + particles.z .= on_architecture(arch, v) + elseif n == :scalefactors + particles.scalefactors .= on_architecture(arch, v) + else + particles.fields[n] .= on_architecture(arch, v) + end + end + + return nothing +end + diff --git a/test/test_particles.jl b/test/test_particles.jl index 176ca7213..746f80660 100644 --- a/test/test_particles.jl +++ b/test/test_particles.jl @@ -31,6 +31,10 @@ grid = RectilinearGrid(architecture; size = (3, 3, 3), extent = (3, 3, 3)) time_step!(model, 1) @test all(particles.fields.A .== 1) + + set!(particles, x = 1, A = fill(5, 3)) + + @test all(particles.fields.A .== 5) & all(particles.x .== 1) end @inline required_tracers(::SimpleParticleBiogeochemistry) = (:B, ) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl deleted file mode 100644 index 071fd47ee..000000000 --- a/test/test_slatissima.jl +++ /dev/null @@ -1,98 +0,0 @@ -include("dependencies_for_runtests.jl") - -using OceanBioME: SLatissima -using Oceananigans.Units -using Oceananigans.Fields: TracerFields -using Oceananigans.Architectures: on_architecture - -function intercept_tracer_tendencies!(model, intercepted_tendencies) - for (name, field) in enumerate(intercepted_tendencies) - field .= Array(interior(model.timestepper.Gⁿ[name + 3])) - end -end - -sum_tracer_nitrogen(tracers) = sum(Array(interior(tracers.NO₃))) + - sum(Array(interior(tracers.NH₄))) + - sum(Array(interior(tracers.P))) + - sum(Array(interior(tracers.Z))) + - sum(Array(interior(tracers.sPON))) + - sum(Array(interior(tracers.bPON))) + - sum(Array(interior(tracers.DON))) - -sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = - sum(Array(interior(tracers.sPOC))) + - sum(Array(interior(tracers.bPOC))) + - sum(Array(interior(tracers.DOC))) + - sum(Array(interior(tracers.DIC))) + - sum(Array(interior(tracers.P)) .* (1 + organic_carbon_calcate_ratio) .+ Array(interior(tracers.Z))) .* redfield - -@testset "SLatissima particle setup and conservations" begin - grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) - - # Initial properties - - particles = SLatissima(; architecture, - x = on_architecture(architecture, ones(Float64, 2) .* .5), - y = on_architecture(architecture, ones(Float64, 2) .* .5), - z = on_architecture(architecture, - ones(Float64, 2) .* .5), - A = on_architecture(architecture, ones(Float64, 2) .* 5), - N = on_architecture(architecture, ones(Float64, 2)), - C = on_architecture(architecture, ones(Float64, 2)), - latitude = 1.0) - - @test length(particles) == 2 - - # nitrogen and carbon conservation - - model = NonhydrostaticModel(; grid, - biogeochemistry = LOBSTER(; grid, carbonates = true, - variable_redfield = true, - sinking_speeds = NamedTuple(), - particles), - advection = nothing, - tracers = (:T, :S)) - - set!(model, NO₃ = 10.0, NH₄ = 1.0, DIC = 2000, Alk = 2000, T = 10, S = 35) - - initial_tracer_N = sum_tracer_nitrogen(model.tracers) - initial_kelp_N = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.N) .+ particles.structural_nitrogen)) ./ (14 * 0.001) - - initial_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield, model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) - initial_kelp_C = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.C) .+ particles.structural_carbon)) ./ (12 * 0.001) - - model.clock.time = 60days # get to a high growth phase - - for _ in 1:10 - time_step!(model, 1) - end - - final_tracer_N = sum_tracer_nitrogen(model.tracers) - final_kelp_N = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.N) .+ particles.structural_nitrogen)) ./ (14 * 0.001) - - final_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield, model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) - final_kelp_C = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.C) .+ particles.structural_carbon)) ./ (12 * 0.001) - - # kelp is being integrated - @test initial_kelp_N != final_kelp_N - @test initial_kelp_C != final_kelp_C - - # conservaitons - # (GPU eps is much larger (~10⁻⁷) than on CPU) - rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_N + initial_kelp_N), √eps(final_tracer_N + final_kelp_N)), 2e-7) - @test isapprox(initial_tracer_N + initial_kelp_N, final_tracer_N + final_kelp_N; rtol) - - rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_C + initial_kelp_C), √eps(final_tracer_C + final_kelp_C)), 7e-7) - @test isapprox(initial_tracer_C + initial_kelp_C, final_tracer_C + final_kelp_C; rtol) - - simulation = Simulation(model, Δt = 1.0, stop_iteration = 1) - - # slow but easiest to have this just done on CPU - intercepted_tendencies = Tuple(Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) - - simulation.callbacks[:intercept_tendencies] = Callback(intercept_tracer_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) - - run!(simulation) - - # the model is changing the tracer tendencies - not sure this test actually works as it didn't fail when it should have - @test any([any(intercepted_tendencies[idx] .!= Array(interior(model.timestepper.Gⁿ[tracer]))) for (idx, tracer) in enumerate((:NO₃, :NH₄, :DIC, :DOC, :bPON, :bPOC))]) -end diff --git a/test/test_sugar_kelp.jl b/test/test_sugar_kelp.jl index 113ab1516..1d95ac452 100644 --- a/test/test_sugar_kelp.jl +++ b/test/test_sugar_kelp.jl @@ -51,7 +51,6 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = Cₛ = particles.biogeochemistry.structural_carbon kₐ = particles.biogeochemistry.structural_dry_weight_per_area - initial_tracer_N = sum_tracer_nitrogen(model.tracers) initial_kelp_N = sum(A .* kₐ .* (N .+ Nₛ)) / (14 * 0.001) @@ -81,6 +80,7 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = # conservaitons # (GPU eps is much larger (~10⁻⁷) than on CPU), is this true??? And is this conservation good enough??? + # this is definitly too high since I didn't catch a sign error in the nitrogen uptakes rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_N + initial_kelp_N), √eps(final_tracer_N + final_kelp_N)), 2e-7) @test isapprox(initial_tracer_N + initial_kelp_N, final_tracer_N + final_kelp_N; rtol) From f0a927a6cd6a25a6c9fa2eeb2b3b977c615fb3ef Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 3 Oct 2024 23:29:18 +0100 Subject: [PATCH 08/33] labels were wrong way around --- examples/kelp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/kelp.jl b/examples/kelp.jl index efef43f9e..617b406ab 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -206,10 +206,10 @@ axis_kwargs = (xlabel = "Time (days)", limits = ((0, times[end] / days), nothing ax1 = Axis(fig[1, 1]; ylabel = "Frond area (dm²)", axis_kwargs...) [lines!(ax1, times / day, A[n, :], linewidth = 3) for n in 1:5] -ax2 = Axis(fig[2, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) +ax2 = Axis(fig[2, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) [lines!(ax2, times / day, (@. A * (N + Nₛ) * kₐ)[n, :], linewidth = 3) for n in 1:5] -ax3 = Axis(fig[3, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) +ax3 = Axis(fig[3, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) [lines!(ax3, times / day, (@. A * (C + Cₛ) * kₐ)[n, :], linewidth = 3) for n in 1:5] fig From 23f3dae6a2de8e772e630d315f5fb6d1043649e5 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 14:24:31 +0100 Subject: [PATCH 09/33] fix (?) adapt --- src/Particles/Particles.jl | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 5bb3ff7d0..67bea3224 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -30,6 +30,13 @@ struct BiogeochemicalParticles{N, B, A, F, T, I, S, X, Y, Z} x :: X # to maintain compatability with lagrangian particle advection from Oceananigans y :: Y z :: Z + + BiogeochemicalParticles(number, biogeochemistry::B, advection::A, + fields::F, timestepper::T, + field_interpolation::I, scalefactors::S, + x::X, y::Y, z::Z) where {B, A, F, T, I, S, X, Y, Z} = + new{number, B, A, F, T, I, S, X, Y, Z}(biogeochemistry, advection, fields, timestepper, + field_interpolation, scalefactors, x, y, z) end architecture(p::BiogeochemicalParticles) = architecture(p.x) @@ -88,16 +95,16 @@ function BiogeochemicalParticles(number; x, y, z) end -Adapt.adapt_structure(to, p::BiogeochemicalParticles) = - BiogeochemicalParticles(adapt(to, p.biogeochemistry), - nothing, - nothing, - nothing, - adapt(to, p.field_interpolation), - nothing, - adapt(to, p.x), - adapt(to, p.y), - adapt(to, p.z)) +Adapt.adapt_structure(to, p::BiogeochemicalParticles{N}) where N = + BiogeochemicalParticles(N, adapt(to, p.biogeochemistry), + nothing, + adapt(to, p.fields), + nothing, + adapt(to, p.field_interpolation), + adapt(to, p.scalefactors), + adapt(to, p.x), + adapt(to, p.y), + adapt(to, p.z)) # Type piracy...oops @inline step_lagrangian_particles!(::Nothing, From 38ac6d8bfe14b409eb083c06c1a390982a480ddc Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 14:30:08 +0100 Subject: [PATCH 10/33] adapts broke some stuff --- src/Particles/Particles.jl | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 67bea3224..d56fddb2b 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -59,11 +59,11 @@ include("set.jl") function BiogeochemicalParticles(number; grid, - biogeochemistry::B, - advection::A = LagrangianAdvection(), + biogeochemistry, + advection = LagrangianAdvection(), timestepper = ForwardEuler, - field_interpolation::I = NearestPoint(), - scalefactors = ones(number)) where {B, A, I} + field_interpolation = NearestPoint(), + scalefactors = ones(number)) arch = architecture(grid) @@ -79,20 +79,14 @@ function BiogeochemicalParticles(number; scalefactors = on_architecture(arch, scalefactors) - F = typeof(fields) - T = typeof(timestepper) - S = typeof(scalefactors) - X = typeof(x) - Y = typeof(y) - Z = typeof(z) - - return BiogeochemicalParticles{number, B, A, F, T, I, S, X, Y, Z}(biogeochemistry, - advection, - fields, - timestepper, - field_interpolation, - scalefactors, - x, y, z) + return BiogeochemicalParticles(number, + biogeochemistry, + advection, + fields, + timestepper, + field_interpolation, + scalefactors, + x, y, z) end Adapt.adapt_structure(to, p::BiogeochemicalParticles{N}) where N = From ab332aa375cc3aee1a48bce6f6ae747da6c4950c Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 14:41:19 +0100 Subject: [PATCH 11/33] some tidying --- test/test_particles.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/test/test_particles.jl b/test/test_particles.jl index 746f80660..2fbc3e2b7 100644 --- a/test/test_particles.jl +++ b/test/test_particles.jl @@ -1,4 +1,4 @@ -#include("dependencies_for_runtests.jl") +include("dependencies_for_runtests.jl") using OceanBioME.Particles: BiogeochemicalParticles @@ -94,9 +94,8 @@ coupled_tracers(::SimpleParticleBiogeochemistry) = (:B, ) set!(model, B = 1) # since the 0, 0, 0 point is ambigously closest to both 3, 3, 3 and 1, 1, 3 (and the logic makes it go to 3, 3, 3) - particles.x .= 0.5 - particles.y .= 0.5 - + set!(particles, x = 0.5, y = 0.5) + time_step!(model, 1) @test all(particles.fields.A .== 0.1) @@ -120,8 +119,7 @@ end set!(model, B = 1) # since the 0, 0, 0 point is ambigously closest to both 3, 3, 3 and 1, 1, 3 (and the logic makes it go to 3, 3, 3) - particles.x .= 0.5 - particles.y .= 0.5 + set!(particles, x = 0.5, y = 0.5) time_step!(model, 1) From cf55efe9a595ded2f928fb1808db9a52655b2039 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 14:54:47 +0100 Subject: [PATCH 12/33] strange bug in particle tests, seems to be related to order functions are defined in --- test/test_particles.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_particles.jl b/test/test_particles.jl index 2fbc3e2b7..6614a744d 100644 --- a/test/test_particles.jl +++ b/test/test_particles.jl @@ -6,7 +6,7 @@ import OceanBioME.Particles: required_particle_fields, required_tracers, coupled struct SimpleParticleBiogeochemistry end -@inline (::SimpleParticleBiogeochemistry)(val_name, t, A) = one(t) +@inline (::SimpleParticleBiogeochemistry)(::Val{:A}, t, A) = one(t) @inline required_particle_fields(::SimpleParticleBiogeochemistry) = (:A, ) @@ -39,7 +39,7 @@ end @inline required_tracers(::SimpleParticleBiogeochemistry) = (:B, ) -@inline (::SimpleParticleBiogeochemistry)(val_name, t, A, B) = 0.1 * B +@inline (::SimpleParticleBiogeochemistry)(::Val{:A}, t, A, B) = 0.1 * B @testset "Testing particle tracer detection" begin particle_biogeochemistry = SimpleParticleBiogeochemistry() From 00002adde84086041bf758088014f24c20872102 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 15:02:01 +0100 Subject: [PATCH 13/33] some clean up --- examples/kelp.jl | 5 +++-- src/Models/Individuals/SugarKelp/SugarKelp.jl | 4 ++-- src/Models/Models.jl | 2 +- src/OceanBioME.jl | 2 +- test/test_sugar_kelp.jl | 12 +++++++----- 5 files changed, 14 insertions(+), 11 deletions(-) diff --git a/examples/kelp.jl b/examples/kelp.jl index 617b406ab..c8ad85f9d 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -198,6 +198,7 @@ end Nₛ = particles.biogeochemistry.structural_nitrogen Cₛ = particles.biogeochemistry.structural_carbon kₐ = particles.biogeochemistry.structural_dry_weight_per_area +sf = particles.scalefactors[1] fig = Figure(size = (1000, 800), fontsize = 20) @@ -207,9 +208,9 @@ ax1 = Axis(fig[1, 1]; ylabel = "Frond area (dm²)", axis_kwargs...) [lines!(ax1, times / day, A[n, :], linewidth = 3) for n in 1:5] ax2 = Axis(fig[2, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) -[lines!(ax2, times / day, (@. A * (N + Nₛ) * kₐ)[n, :], linewidth = 3) for n in 1:5] +[lines!(ax2, times / day, (@. A * (N + Nₛ) * kₐ * sf)[n, :], linewidth = 3) for n in 1:5] ax3 = Axis(fig[3, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) -[lines!(ax3, times / day, (@. A * (C + Cₛ) * kₐ)[n, :], linewidth = 3) for n in 1:5] +[lines!(ax3, times / day, (@. A * (C + Cₛ) * kₐ * sf)[n, :], linewidth = 3) for n in 1:5] fig diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index 01c1f95c6..394fce487 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -1,6 +1,6 @@ module SugarKelpModel -export SugarKelp +export SugarKelp, SugarKelpParticles using Roots @@ -62,7 +62,7 @@ import OceanBioME.Particles: required_particle_fields, required_tracers, coupled end # convenience constructor that sets up with default parameters -SugarKelp(number; grid, kelp_parameters = NamedTuple(), kwargs...) = +SugarKelpParticles(number; grid, kelp_parameters = NamedTuple(), kwargs...) = BiogeochemicalParticles(number; grid, biogeochemistry = SugarKelp(; kelp_parameters...), kwargs...) @inline required_particle_fields(::SugarKelp) = (:A, :N, :C) diff --git a/src/Models/Models.jl b/src/Models/Models.jl index 20cb4a628..6184e17a9 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -6,7 +6,7 @@ export NPZD, NutrientPhytoplanktonZooplanktonDetritus, LOBSTER -export SugarKelp +export SugarKelp, SugarKelpParticles export CarbonChemistry diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 38ff25b1f..d1b758d32 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -8,7 +8,7 @@ module OceanBioME export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, NPZD, redfield # Macroalgae models -export BiogeochemicalParticles, SugarKelp +export BiogeochemicalParticles, SugarKelp, SugarKelpParticles # Box model export BoxModel, BoxModelGrid, SpeedyOutput, load_output diff --git a/test/test_sugar_kelp.jl b/test/test_sugar_kelp.jl index 1d95ac452..ee2d9b4e5 100644 --- a/test/test_sugar_kelp.jl +++ b/test/test_sugar_kelp.jl @@ -1,3 +1,5 @@ +include("dependencies_for_runtests.jl") + using Oceananigans.Architectures: on_architecture sum_tracer_nitrogen(tracers) = sum(on_architecture(CPU(), interior(tracers.NO₃))) + @@ -18,11 +20,11 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = @testset "SLatissima particle setup and conservations" begin grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) - particle_biogeochemistry = OceanBioME.Models.SugarKelpModel.SugarKelp() + particle = SugarKelpParticles(2; grid, advection = nothing) - particles = BiogeochemicalParticles(2; grid, - biogeochemistry = particle_biogeochemistry, - advection = nothing) + @test particles isa BiogeochemicalParticles + @test particles.biogeochemistry isa SugarKelp + @test length(particles) == 2 biogeochemistry = LOBSTER(; grid, particles, @@ -80,7 +82,7 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = # conservaitons # (GPU eps is much larger (~10⁻⁷) than on CPU), is this true??? And is this conservation good enough??? - # this is definitly too high since I didn't catch a sign error in the nitrogen uptakes + # this is definitly too high since I didn't catch a sign error of an uptake was the wrong way around rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_N + initial_kelp_N), √eps(final_tracer_N + final_kelp_N)), 2e-7) @test isapprox(initial_tracer_N + initial_kelp_N, final_tracer_N + final_kelp_N; rtol) From 71d8614c9a5af3d561d202eca37d87df96a1990b Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 15:03:02 +0100 Subject: [PATCH 14/33] docs fix --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index e3cc9c0a3..606aebc97 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -53,7 +53,7 @@ if !isdir(OUTPUT_DIR) mkdir(OUTPUT_DIR) end model_parameters = (LOBSTER(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry, - SugarKelp().biogeochemistry, + SugarKelp(), TwoBandPhotosyntheticallyActiveRadiation(; grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))), SimpleMultiG(; grid = BoxModelGrid()), InstantRemineralisation(; grid = BoxModelGrid()), From bd2f13678449eca940bc4e40a157fc51e05fa79e Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 15:44:42 +0100 Subject: [PATCH 15/33] forgot to update example --- examples/kelp.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/kelp.jl b/examples/kelp.jl index c8ad85f9d..8ea20ac5e 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -60,9 +60,9 @@ S = ConstantField(35) n = 5 # number of kelp bundles z₀ = [-21:5:-1;] * 1.0 # depth of kelp fronds -particles = SugarKelp(n; grid, - advection = nothing, # we don't want them to move around - scalefactors = fill(2000, n)) # and we want them to look like there are 500 in each bundle +particles = SugarKelpParticles(n; grid, + advection = nothing, # we don't want them to move around + scalefactors = fill(2000, n)) # and we want them to look like there are 500 in each bundle set!(particles, A = 10, N = 0.01, C = 0.1, z = z₀, x = Lx / 2, y = Ly / 2) From 080ea9c45c164f18fc8d931592e7cd3b3af34441 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 17:11:31 +0100 Subject: [PATCH 16/33] updated api library --- docs/src/appendix/library.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/appendix/library.md b/docs/src/appendix/library.md index 64b5ac576..29b51c35a 100644 --- a/docs/src/appendix/library.md +++ b/docs/src/appendix/library.md @@ -24,7 +24,7 @@ Modules = [OceanBioME.Models.LOBSTERModel] ### Sugar kelp (Saccharina latissima) ```@autodocs -Modules = [OceanBioME.Models.SLatissimaModel] +Modules = [OceanBioME.Models.SugarKelpModel] ``` ### Carbon Chemistry From d8b001688a62cfd4a99a5129e8e28c67ccc43e00 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 18:02:56 +0100 Subject: [PATCH 17/33] updated particle and sugar kelp docs --- .../src/model_components/individuals/index.md | 88 +++++-------------- .../individuals/slatissima.md | 25 +++++- 2 files changed, 47 insertions(+), 66 deletions(-) diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index a42988ccf..f06745486 100644 --- a/docs/src/model_components/individuals/index.md +++ b/docs/src/model_components/individuals/index.md @@ -1,92 +1,50 @@ # [Individuals](@id individuals) -The effects of individuals can be modelled in OceanBioME. We have implemented this through custom dynamics in the [Lagrangian Particle tracking feature of Oceananigans](https://clima.github.io/OceananigansDocumentation/stable/model_setup/lagrangian_particles/). We have extended these functionalities to make it easier to implement "active" particles which interact with the tracers. We have then implemented a model of [sugar kelp](@ref SLatissima) which can be followed as an example of using this functionality. +The effects of individuals can be modelled in OceanBioME. We have implemented this through custom dynamics in the [Lagrangian Particle tracking feature of Oceananigans](https://clima.github.io/OceananigansDocumentation/stable/model_setup/lagrangian_particles/). We have extended these functionalities to make it easier to implement "active" particles which interact with the tracers. We have then implemented a model of [sugar kelp](@ref sugar-kelp) which can be followed as an example of using this functionality. -To setup particles first create a particle type with the desired properties, e.g.: +To setup particles first create a particle biogeochemistry, e.g.: ```@example particles -using OceanBioME.Particles: BiogeochemicalParticles -struct GrowingParticles{FT, VT} <: BiogeochemicalParticles +struct GrowingParticles{FT} nutrients_half_saturation :: FT - - size :: VT - nitrate_uptake :: VT - - x :: VT - y :: VT - z :: VT end ``` -You then need to overload particular functions to integrate the growth, so they need to first be `import`ed: +We then need to add some methods to tell `OceanBioME` what properties this particle has, and what tracers it interacts with: ```@example particles -import Oceananigans.Biogeochemistry: update_tendencies! -import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties! -``` -First, to integrate the particles properties we overload `update_lagrangian_particle_properties!`; -in this fictitious case we will have a Mondo-quota nutrient uptake and growth: - -```@example particles -using Oceananigans.Fields: interpolate +import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers -function update_lagrangian_particle_properties!(particles::GrowingParticles, model, bgc, Δt) - @inbounds for p in 1:length(particles) - nutrients = @inbounds interpolate(model.tracers.NO₃, particle.x[p], particle.y[p], particle.z[p]) +required_particle_fields(::GrowingParticles) = (:S, ) +required_tracers(::GrowingParticles) = (:N, ) +coupled_tracers(::GrowingParticles) = (:N, ) - uptake = nutrients / (particle.nutrients_half_saturation + nutrients) - - particles.size[p] += uptake * Δt - particles.nitrate_uptake[p] = uptake - end - return nothing -end - -nothing #hide ``` -In this example the particles will not move around, and are only integrated on a single thread. For a more comprehensive example see the [Sugar Kelp](@ref SLatissima) implementation. We then need to update the tracer tendencies to match the nutrients' uptake: - +So our model is going to track the `S`ize of the particles and take up `N`utrients. +Now we need to how this growth happens. +The forcing functions should be of the form `(particles::ParticleBiogeochemistry)(::Val{:PROPERTY}, t, required_particle_fields..., required_tracers...)`, so in this example: ```@example particles -using OceanBioME.Particles: get_node - -function update_tendencies!(bgc, particles::GrowingParticles, model) - @inbounds for p in 1:length(particles) - i, j, k = fractional_indices((x, y, z), grid, Center(), Center(), Center()) - - # Convert fractional indices to unit cell coordinates 0 ≤ (ξ, η, ζ) ≤ 1 - # and integer indices (with 0-based indexing). - ξ, i = modf(i) - η, j = modf(j) - ζ, k = modf(k) +(p::GrowingParticles)(::Val{:S}, t, S, N) = N / (N + p.nutrient_half_saturation) +(p::GrowingParticles)(::Val{:N}, t, S, N) = - N / (N + p.nutrient_half_saturation) +``` - # Round to nearest node and enforce boundary conditions - i, j, k = (get_node(TX(), Int(ifelse(ξ < 0.5, i + 1, i + 2)), grid.Nx), - get_node(TY(), Int(ifelse(η < 0.5, j + 1, j + 2)), grid.Ny), - get_node(TZ(), Int(ifelse(ζ < 0.5, k + 1, k + 2)), grid.Nz)) +We can then create an instance of this particle model using `BiogeochemicalParticles`, and set their initial position and size: +```@example particles +using OceanBioME - node_volume = volume(i, j, k, grid, Center(), Center(), Center()) +Lx, Ly, Lz = 100, 100, 100 +grid = RectilinearGrid(; size = (8, 8, 8), extent = (Lx, Ly, Lz)) - model.timestepper.Gⁿ.NO₃[i, j, k] += particles.nitrate_uptake[p] / (d * node_volume) - end - return nothing -end +particles = BiogeochemicalParticles(10; grid, biogeochemistry = GrowingParticles()) -nothing #hide +set!(particles, S = 0.1, x = rand(10) * Lx, y = rand(10) * Ly, z = rand(10) * Lz) ``` -Now we can just plug this into any biogeochemical model setup to have particles (currently [NPZD](@ref NPZD) and [LOBSTER](@ref LOBSTER)): - +We can then put these into a compatible biogeochemical model, for example: ```@example particles -using OceanBioME, Oceananigans - -Lx, Ly, Lz = 1000, 1000, 100 -grid = RectilinearGrid(; size = (64, 64, 16), extent = (Lx, Ly, Lz)) - -# Start the particles randomly distributed, floating on the surface -particles = GrowingParticles(0.5, zeros(3), zeros(3), rand(3) * Lx, rand(3) * Ly, zeros(3)) -biogeochemistry = LOBSTER(; grid, particles) +biogeochemistry = NPZD(; grid, particles) ``` diff --git a/docs/src/model_components/individuals/slatissima.md b/docs/src/model_components/individuals/slatissima.md index eb18809f2..fe95bbff9 100644 --- a/docs/src/model_components/individuals/slatissima.md +++ b/docs/src/model_components/individuals/slatissima.md @@ -1,4 +1,4 @@ -# [Sugar kelp (Saccharina latissima) individuals](@id SLatissima) +# [Sugar kelp (Saccharina latissima) individuals](@id sugar-kelp) We have implemented a model of sugar kelp growth within this spatially infinitesimal Lagrangian particles framework originally based on the model of [Broch2012](@citet) and updated by [Broch2013](@citet), [Fossberg2018](@citet), and [Broch2019](@citet). This is the same model passively forced by [StrongWright2022](@citet). @@ -7,6 +7,29 @@ The model tracks three variables, the frond area, A (dm²), carbon reserve, C (g Results could look something like this (from [StrongWright2022](@citet)): ![Example A, N, and C profiles from [StrongWright2022](@citet)](https://www.frontiersin.org/files/Articles/793977/fmars-08-793977-HTML/image_m/fmars-08-793977-g002.jpg) +You can access the model biogeochemistry by setting up `SugarKelp`, i.e.: +```jldoctest +using OceanBioME + +kelp_bgc = SugarKelp() + +# output +SugarKelp{FT} biogeochemistry (Broch & Slagstad, 2012) tracking the `N`itrogen and `C`arbon in a frond of `A`rea +``` +which can be put into `BiogeochemicalParticles`, or you can directly manifest particles: +```jldoctest +using OceanBioME, Oceananigans + +grid = RectilinearGrid(size = (1, 1, 1), extent = (1, 1, 1)); +particles = SugarKelpParticles(10; grid) + +# output +10 BiogeochemicalParticles with SugarKelp{FT} biogeochemistry: +├── fields: (:A, :N, :C) +└── coupled tracers: (:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON) + +``` + ## Model equations As per [Broch2012](@citet) this model variables evolve as: From 7af67ba445200e7c87ec406a25564d47237ee59d Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 18:11:36 +0100 Subject: [PATCH 18/33] docstrings --- src/Models/Individuals/SugarKelp/SugarKelp.jl | 33 +++++++++++++++++-- src/Particles/Particles.jl | 18 ++++++++++ src/Particles/advection.jl | 5 +++ src/Particles/time_stepping.jl | 6 ++++ src/Particles/tracer_interpolation.jl | 5 +++ 5 files changed, 64 insertions(+), 3 deletions(-) diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index 394fce487..4d6ae9741 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -1,3 +1,20 @@ +""" +Sugar kelp model of [Broch2012](@citet) and updated by [Broch2013](@citet), [Fossberg2018](@citet), and [Broch2019](@citet). + +Prognostic properties +===================== +* Area: A (dm²) +* Nitrogen reserve: N (gN/gSW) +* Carbon reserve: C (gC/gSW) + +Tracer dependencies +=================== +* Nitrates: NO₃ (mmol N/m³) +* Ammonia: NH₄ (mmol N/m³) +* Photosynthetically available radiation: PAR (einstein/m²/day) +* Temperature: T (°C) + +""" module SugarKelpModel export SugarKelp, SugarKelpParticles @@ -11,6 +28,11 @@ using OceanBioME.Particles: BiogeochemicalParticles import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers # disgsuting number of parameters +""" + SugarKelp + +Defines the parameters for `SugarKelp` biogeochemistry. +""" @kwdef struct SugarKelp{FT, TL} temperature_limit :: TL = LinearOptimalTemperatureRange() # TODO: split up more parameterisations like this growth_rate_adjustment :: FT = 4.5 @@ -61,9 +83,14 @@ import OceanBioME.Particles: required_particle_fields, required_tracers, coupled adapted_latitude :: FT = 57.5 end -# convenience constructor that sets up with default parameters -SugarKelpParticles(number; grid, kelp_parameters = NamedTuple(), kwargs...) = - BiogeochemicalParticles(number; grid, biogeochemistry = SugarKelp(; kelp_parameters...), kwargs...) +""" + SugarKelpParticles(n; grid, kelp_parameters = NamedTuple(), kwargs...) + +Sets up `n` sugar kelp `BiogeochemicalParticles` with default parameters except +those specified in `kelp_parameters`. `kwagrs` are passed onto `BiogeochemicalParticles`. +""" +SugarKelpParticles(n; grid, kelp_parameters = NamedTuple(), kwargs...) = + BiogeochemicalParticles(n; grid, biogeochemistry = SugarKelp(; kelp_parameters...), kwargs...) @inline required_particle_fields(::SugarKelp) = (:A, :N, :C) @inline required_tracers(::SugarKelp) = (:u, :v, :w, :T, :NO₃, :NH₄, :PAR) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index d56fddb2b..145e38795 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -57,6 +57,24 @@ include("time_stepping.jl") include("update_tracer_tendencies.jl") include("set.jl") +""" + BiogeochemicalParticles(number; + grid, + biogeochemistry, + advection = LagrangianAdvection(), + timestepper = ForwardEuler, + field_interpolation = NearestPoint(), + scalefactors = ones(number)) + +Creates `number` particles with `biogeochemistry` on `grid`, advected by +`advection` which defaults to `LagrangianAdvection` (i.e. they comove with +the water). The biogeochemistry is stepped by `timestepper` and tracer fields +are interpolated by `field_interpolation`, which defaults to directly reading +the nearest center point and taking up/depositing in the same. + +Particles can also have a `scalefactor` which scales their tracer interaction +(e.g. to mimic the particle representing multiple particles). +""" function BiogeochemicalParticles(number; grid, biogeochemistry, diff --git a/src/Particles/advection.jl b/src/Particles/advection.jl index 1ac64208f..b940def40 100644 --- a/src/Particles/advection.jl +++ b/src/Particles/advection.jl @@ -5,6 +5,11 @@ using Oceananigans.Models.LagrangianParticleTracking: _advect_particles!, total_ # put nothing to do nothing advect_particles!(advection, particles, model, Δt) = nothing +""" + LagrangianAdvection + +Specifies that particles should move in a purley lagrangian mannor. +""" struct LagrangianAdvection end function advect_particles!(::LagrangianAdvection, particles, model, Δt) diff --git a/src/Particles/time_stepping.jl b/src/Particles/time_stepping.jl index 0217b26fb..a5b81ea6d 100644 --- a/src/Particles/time_stepping.jl +++ b/src/Particles/time_stepping.jl @@ -1,3 +1,9 @@ +""" + ForwardEuler + +Step particle biogeochemistry with a `ForwardEuler` methods with `Δt` from +the physical model substep. +""" struct ForwardEuler{T} tendencies :: T diff --git a/src/Particles/tracer_interpolation.jl b/src/Particles/tracer_interpolation.jl index aefaab7f9..9c9489c80 100644 --- a/src/Particles/tracer_interpolation.jl +++ b/src/Particles/tracer_interpolation.jl @@ -5,6 +5,11 @@ using Oceananigans.Grids: AbstractGrid, Flat, Bounded, Periodic @inline get_node(::Bounded, i, N) = min(max(i, 1), N) @inline get_node(::Periodic, i, N) = ifelse(i < 1, N, ifelse(i > N, 1, i)) +""" + NearestPoint + +Specifies that tracer values should be taken from the nearst center point. +""" struct NearestPoint end @inline function extract_tracer_values(::NearestPoint, particles, grid, fields, n) From bb2fd2817ac9a7e874ea0d9eb192ed483928c0b2 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 20:51:17 +0100 Subject: [PATCH 19/33] docs fixes --- docs/src/appendix/library.md | 6 ++++++ docs/src/model_components/biogeochemical/LOBSTER.md | 2 +- docs/src/model_components/individuals/index.md | 2 +- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/src/appendix/library.md b/docs/src/appendix/library.md index 29b51c35a..8efdcde71 100644 --- a/docs/src/appendix/library.md +++ b/docs/src/appendix/library.md @@ -56,3 +56,9 @@ Modules = [OceanBioME.Models.GasExchangeModel, OceanBioME.Models.GasExchangeMode ```@autodocs Modules = [OceanBioME.BoxModels] ``` + +## Particles + +```@autodocs +Modules = [OceanBioME.Particles] +``` \ No newline at end of file diff --git a/docs/src/model_components/biogeochemical/LOBSTER.md b/docs/src/model_components/biogeochemical/LOBSTER.md index dd562b901..b04deec94 100644 --- a/docs/src/model_components/biogeochemical/LOBSTER.md +++ b/docs/src/model_components/biogeochemical/LOBSTER.md @@ -103,7 +103,7 @@ When the oxygen chemistry is activated additional tracer ``O_2`` evolve like: ### Variable Redfield -When the variable Redfield modification is activated the organic components are modified to evolve their nitrogen and carbon content separately. This means that the waste from non-Redfield models (e.g. loss from the [kelp](@ref SLatissima)) can be accounted for. +When the variable Redfield modification is activated the organic components are modified to evolve their nitrogen and carbon content separately. This means that the waste from non-Redfield models (e.g. loss from the [kelp](@ref sugar-kelp)) can be accounted for. In this case the organic components are split into nitrogen and carbon compartments, so the tracers ``sPOM``, ``bPOM``, and ``DOM`` are replaced with ``sPON``, ``sPOC``, ``bPON``, ``bPOC``, ``DON``, and ``DOC``. The nitrogen compartments evolve as per the organic matter equations above (i.e. replacing each ``XOM`` with ``XON``), while the carbon compartments evolve like: diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index f06745486..fea676a63 100644 --- a/docs/src/model_components/individuals/index.md +++ b/docs/src/model_components/individuals/index.md @@ -33,7 +33,7 @@ The forcing functions should be of the form `(particles::ParticleBiogeochemistry We can then create an instance of this particle model using `BiogeochemicalParticles`, and set their initial position and size: ```@example particles -using OceanBioME +using OceanBioME, Oceananigans Lx, Ly, Lz = 100, 100, 100 grid = RectilinearGrid(; size = (8, 8, 8), extent = (Lx, Ly, Lz)) From d852f1ebf580be173298c6ac2a97421c3f6fa700 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 4 Oct 2024 21:40:17 +0100 Subject: [PATCH 20/33] final oops --- docs/src/model_components/individuals/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index fea676a63..eed2eb2b4 100644 --- a/docs/src/model_components/individuals/index.md +++ b/docs/src/model_components/individuals/index.md @@ -38,7 +38,7 @@ using OceanBioME, Oceananigans Lx, Ly, Lz = 100, 100, 100 grid = RectilinearGrid(; size = (8, 8, 8), extent = (Lx, Ly, Lz)) -particles = BiogeochemicalParticles(10; grid, biogeochemistry = GrowingParticles()) +particles = BiogeochemicalParticles(10; grid, biogeochemistry = GrowingParticles(0.5)) set!(particles, S = 0.1, x = rand(10) * Lx, y = rand(10) * Ly, z = rand(10) * Lz) ``` From fbf28deae093c3e5180f4da603ca4bbb63a49d93 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sat, 5 Oct 2024 11:41:30 +0100 Subject: [PATCH 21/33] updated index and attempted to make kelp example figure show gain --- docs/src/index.md | 2 +- examples/kelp.jl | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 026553160..1c6423541 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,7 +4,7 @@ OceanBioME.jl is a fast and flexible ocean biogeochemical modelling environment. OceanBioME.jl currently provides a core of several biogeochemical models Nutrient--Phytoplankton--Zooplankton--Detritus (NPZD) and [LOBSTER](https://doi.org/10.1029/2004JC002588), a medium complexity model, air-sea gas exchange models to provide appropriate top boundary conditions, and sediment models to for the benthic boundary. [PISCES](https://doi.org/10.5194/gmd-8-2465-2015) and other higher complexity models are in our future development plans. -OceanBioME.jl includes a framework for integrating the growth of biological/active Lagrangian particles which move around and can interact with the (Eulerian) tracer fields - for example, consuming nutrients and carbon dioxide while releasing dissolved organic material. A growth model for sugar kelp is currently implemented using active particles, and this model can be used in a variety of dynamical scenarios including free-floating or bottom-attached particles. +OceanBioME.jl includes a framework for integrating the growth of [biological/active particles](@ref individuals) which move around and can interact with the (Eulerian) tracer fields - for example, consuming nutrients and carbon dioxide while releasing dissolved organic material. A growth model for sugar kelp is currently implemented using active particles, and this model can be used in a variety of dynamical scenarios including free-floating or bottom-attached particles. ## Quick install diff --git a/examples/kelp.jl b/examples/kelp.jl index 8ea20ac5e..e12a51a78 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -145,7 +145,7 @@ end using CairoMakie -fig = Figure(size = (1000, 1200), fontsize = 20) +fig = Figure(size = (1000, 900), fontsize = 20) axis_kwargs = (xlabel = "Time (days)", ylabel = "z (m)", limits = ((0, times[end] / days), (-85meters, 0))) @@ -157,10 +157,6 @@ axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N/m³)", axis_kwa hmNO₃ = heatmap!(times / days, z, interior(tracers["NO₃"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[2, 2], hmNO₃) -axZ = Axis(fig[3, 1]; title = "Zooplankton concentration (mmol N/m³)", axis_kwargs...) -hmZ = heatmap!(times / days, z, interior(tracers["Z"], 1, 1, :, :)', colormap = :batlow) -Colorbar(fig[3, 2], hmZ) - axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol C/m³)", axis_kwargs...) hmD = heatmap!(times / days, z, interior(tracers["sPOC"], 1, 1, :, :)' .+ interior(tracers["bPOC"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[4, 2], hmD) From b849e93f87cbbb7e8b6adf2a6e5597a989467aa0 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sat, 5 Oct 2024 11:42:10 +0100 Subject: [PATCH 22/33] better attempt --- examples/kelp.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/kelp.jl b/examples/kelp.jl index e12a51a78..dbefd4d9a 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -169,7 +169,9 @@ lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year, lines!(axfDIC, times / days, carbon_export / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export") Legend(fig[5, 2], axfDIC, framevisible = false) -fig +save("kelp.png", fig) + +# ![](kelp.png) # We can also have a look at how the kelp particles evolve using JLD2 From bf6603e04fc4fbf1083cb6f5e53fc66be8998a65 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sat, 5 Oct 2024 14:11:46 +0100 Subject: [PATCH 23/33] fix figure --- examples/kelp.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/kelp.jl b/examples/kelp.jl index dbefd4d9a..e877ec495 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -157,17 +157,17 @@ axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N/m³)", axis_kwa hmNO₃ = heatmap!(times / days, z, interior(tracers["NO₃"], 1, 1, :, :)', colormap = :batlow) Colorbar(fig[2, 2], hmNO₃) -axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol C/m³)", axis_kwargs...) +axD = Axis(fig[3, 1]; title = "Detritus concentration (mmol C/m³)", axis_kwargs...) hmD = heatmap!(times / days, z, interior(tracers["sPOC"], 1, 1, :, :)' .+ interior(tracers["bPOC"], 1, 1, :, :)', colormap = :batlow) -Colorbar(fig[4, 2], hmD) +Colorbar(fig[3, 2], hmD) CO₂_molar_mass = (12 + 2 * 16) * 1e-3 # kg / mol -axfDIC = Axis(fig[5, 1], xlabel = "Time (days)", ylabel = "Flux (kgCO₂/m²/year)", +axfDIC = Axis(fig[4, 1], xlabel = "Time (days)", ylabel = "Flux (kgCO₂/m²/year)", title = "Air-sea CO₂ flux and Sinking", limits = ((0, times[end] / days), nothing)) lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Air-sea flux") lines!(axfDIC, times / days, carbon_export / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export") -Legend(fig[5, 2], axfDIC, framevisible = false) +Legend(fig[4, 2], axfDIC, framevisible = false) save("kelp.png", fig) From 23e5ea721807b1deb42ad8e733a980f808871141 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sat, 5 Oct 2024 14:13:28 +0100 Subject: [PATCH 24/33] changed units --- examples/kelp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/kelp.jl b/examples/kelp.jl index e877ec495..e68a03a93 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -208,7 +208,7 @@ ax1 = Axis(fig[1, 1]; ylabel = "Frond area (dm²)", axis_kwargs...) ax2 = Axis(fig[2, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) [lines!(ax2, times / day, (@. A * (N + Nₛ) * kₐ * sf)[n, :], linewidth = 3) for n in 1:5] -ax3 = Axis(fig[3, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) -[lines!(ax3, times / day, (@. A * (C + Cₛ) * kₐ * sf)[n, :], linewidth = 3) for n in 1:5] +ax3 = Axis(fig[3, 1]; ylabel = "Total Carbon (kgCO₂(eq))", axis_kwargs...) +[lines!(ax3, times / day, (@. A * (C + Cₛ) * kₐ * sf)[n, :] / 1000 * 44 / 12, linewidth = 3) for n in 1:5] fig From 84b927a3c32c011f23cb9ac293b82e0d83b8ba0f Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 7 Oct 2024 15:27:50 +0100 Subject: [PATCH 25/33] formatting --- docs/src/model_components/biogeochemical/LOBSTER.md | 2 +- src/Models/Models.jl | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/docs/src/model_components/biogeochemical/LOBSTER.md b/docs/src/model_components/biogeochemical/LOBSTER.md index b04deec94..3fb3f584b 100644 --- a/docs/src/model_components/biogeochemical/LOBSTER.md +++ b/docs/src/model_components/biogeochemical/LOBSTER.md @@ -95,7 +95,7 @@ When the carbonate chemistry is activated additional tracers ``DIC`` and ``Alk`` ### Oxygen chemistry -When the oxygen chemistry is activated additional tracer ``O_2`` evolve like: +When the oxygen chemistry is activated, additional tracer ``O_2`` evolve like: ```math \frac{\partial O_2}{\partial t} = \mu_P L_{PAR}\left(L_{NO_3} + L_{NH_4}\right)R_{O_2}P - (R_{O_2} - R_{nit})\frac{\partial NH_4}{\partial t} - R_{O_2}\mu_nNH_4. diff --git a/src/Models/Models.jl b/src/Models/Models.jl index 6184e17a9..1c0d28907 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -6,7 +6,7 @@ export NPZD, NutrientPhytoplanktonZooplanktonDetritus, LOBSTER -export SugarKelp, SugarKelpParticles +export SugarKelp, SugarKelpParticles, GiantKelp export CarbonChemistry @@ -22,7 +22,6 @@ export GasExchange, include("Sediments/Sediments.jl") include("AdvectedPopulations/LOBSTER/LOBSTER.jl") include("AdvectedPopulations/NPZD.jl") -#include("Individuals/SLatissima.jl") include("Individuals/SugarKelp/SugarKelp.jl") include("seawater_density.jl") include("CarbonChemistry/CarbonChemistry.jl") From 563c9ccf2aaaffbee470c2810a524486d28f8642 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 7 Oct 2024 15:39:44 +0100 Subject: [PATCH 26/33] oops --- src/Models/Individuals/SugarKelp/equations.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Models/Individuals/SugarKelp/equations.jl b/src/Models/Individuals/SugarKelp/equations.jl index b064f5a2a..1c8d8e25d 100644 --- a/src/Models/Individuals/SugarKelp/equations.jl +++ b/src/Models/Individuals/SugarKelp/equations.jl @@ -1,4 +1,3 @@ - @inline function (kelp::SugarKelp)(::Val{:A}, t, A, N, C, u, v, w, T, NO₃, NH₄, PAR) μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) @@ -26,7 +25,7 @@ end Cₛ = kelp.structural_carbon P = photosynthesis(kelp, T, PAR) - + μ = growth(kelp, t, A, N, C, T, NH₄, u, v, w) R = respiration(kelp, t, A, N, C, T, NO₃, NH₄, u, v, w, μ) @@ -185,6 +184,8 @@ end fₜ = kelp.temperature_limit(T) fₐ = area_limitation(kelp, A) fₚ = seasonal_limitation(kelp, t) + + return fₜ * fₐ * fₚ end @inline function area_limitation(kelp, A) @@ -200,7 +201,7 @@ end upper_optimal :: FT = 15.0 lower_gradient :: FT = 1/(lower_optimal + 1.8) #0.08 - I think its important that the minimum cut off (-1.8) is observed, # because the literature doesn't say anything about it growing okay and then suddenly dying off - #at low temperature as the origional form suggests, although for higher temperature... + # at low temperature as the origional form suggests, although for higher temperature... upper_gradient :: FT = -0.25 end From f6ceb8d5d78acd044b373d6bb5ae9290bb966819 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 7 Oct 2024 16:01:40 +0100 Subject: [PATCH 27/33] units for PAR oops --- src/Models/Individuals/SugarKelp/equations.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Models/Individuals/SugarKelp/equations.jl b/src/Models/Individuals/SugarKelp/equations.jl index 1c8d8e25d..a59d1c32e 100644 --- a/src/Models/Individuals/SugarKelp/equations.jl +++ b/src/Models/Individuals/SugarKelp/equations.jl @@ -90,6 +90,8 @@ end end @inline function photosynthesis(kelp, T, PAR) + PAR *= day / (3.99e-10 * 545e12) # W / m² / s to einstein / m² / day + Tk = T + 273.15 P₁ = kelp.photosynthesis_at_ref_temp_1 From 708a856b8d99fa8b3ab5896f7693dd76ae77006d Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 8 Oct 2024 14:05:59 +0100 Subject: [PATCH 28/33] factor of 4 issue --- src/Models/Individuals/SugarKelp/SugarKelp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index 4d6ae9741..dfced410f 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -72,8 +72,8 @@ Defines the parameters for `SugarKelp` biogeochemistry. current_speed_for_0p65_uptake :: FT = 0.03 nitrate_half_saturation :: FT = 4.0 ammonia_half_saturation :: FT = 1.3 - maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6) - maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6) + maximum_nitrate_uptake :: FT = 10 / structural_dry_weight_per_area * 24 * 14 / (10^6) # / dm^2 / h to + maximum_ammonia_uptake :: FT = 12 / structural_dry_weight_per_area * 24 * 14 / (10^6) current_1 :: FT = 0.72 current_2 :: FT = 0.28 current_3 :: FT = 0.045 From 200c2fa41bbaac6e6053cdb29e9b96c1f23ae798 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 9 Oct 2024 15:04:08 +0100 Subject: [PATCH 29/33] turns out its just the test being werid --- src/Particles/Particles.jl | 8 +------- test/test_particles.jl | 2 +- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index e706b0155..d71ea5587 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -121,13 +121,7 @@ Adapt.adapt_structure(to, p::BiogeochemicalParticles{N}) where N = const BGC_WITH_PARTICLES = Union{<:DiscreteBiogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, <:ContinuousBiogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}} -@inline step_lagrangian_particles!(::Nothing, - model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}, - Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) - -@inline step_lagrangian_particles!(::Nothing, - model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}, - Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) +@inline step_lagrangian_particles!(::Nothing, model, Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) @inline update_lagrangian_particle_properties!(model, bgc::BGC_WITH_PARTICLES, Δt) = update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) diff --git a/test/test_particles.jl b/test/test_particles.jl index 6614a744d..5581292ac 100644 --- a/test/test_particles.jl +++ b/test/test_particles.jl @@ -74,7 +74,7 @@ end time_step!(model, 1) @test all(particles.fields.A .== 0.1) - @test all(particles.x .== 0.1) && all(particles.y .== 0.2) && all(particles.z .== 0) + @test all(particles.x .== 0.1) && all(particles.y .≈ 0.2) && all(particles.z .== 0) end coupled_tracers(::SimpleParticleBiogeochemistry) = (:B, ) From a4666f1a1f0c555b125580b9db78165dcc4c2539 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 9 Oct 2024 15:46:49 +0100 Subject: [PATCH 30/33] warnings go away please --- src/Particles/Particles.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index d71ea5587..5673c7453 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -121,7 +121,11 @@ Adapt.adapt_structure(to, p::BiogeochemicalParticles{N}) where N = const BGC_WITH_PARTICLES = Union{<:DiscreteBiogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, <:ContinuousBiogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}} -@inline step_lagrangian_particles!(::Nothing, model, Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) +const MODEL_WITH_BGC_PARTICLES = Union{<:NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}, + <:HydrostaticFreeSurfaceModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}} + +@inline step_lagrangian_particles!(::Nothing, model::MODEL_WITH_BGC_PARTICLES, Δt) = + update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) @inline update_lagrangian_particle_properties!(model, bgc::BGC_WITH_PARTICLES, Δt) = update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) From cb1d6c16d1799afe4702784da1e946a11ab24316 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 9 Oct 2024 15:47:13 +0100 Subject: [PATCH 31/33] formatting --- src/Particles/Particles.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 5673c7453..e1d733c8d 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -121,8 +121,9 @@ Adapt.adapt_structure(to, p::BiogeochemicalParticles{N}) where N = const BGC_WITH_PARTICLES = Union{<:DiscreteBiogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, <:ContinuousBiogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}} -const MODEL_WITH_BGC_PARTICLES = Union{<:NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}, - <:HydrostaticFreeSurfaceModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}} +const MODEL_WITH_BGC_PARTICLES = + Union{<:NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}, + <:HydrostaticFreeSurfaceModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:BGC_WITH_PARTICLES}} @inline step_lagrangian_particles!(::Nothing, model::MODEL_WITH_BGC_PARTICLES, Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) From 19a634988ad5bf86b1d195895daa06f97dd992d6 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 11 Oct 2024 13:19:52 +0100 Subject: [PATCH 32/33] numerical errors... --- test/test_particles.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_particles.jl b/test/test_particles.jl index 5581292ac..3e0027d30 100644 --- a/test/test_particles.jl +++ b/test/test_particles.jl @@ -34,7 +34,7 @@ grid = RectilinearGrid(architecture; size = (3, 3, 3), extent = (3, 3, 3)) set!(particles, x = 1, A = fill(5, 3)) - @test all(particles.fields.A .== 5) & all(particles.x .== 1) + @test all(particles.fields.A .≈ 5) & all(particles.x .== 1) end @inline required_tracers(::SimpleParticleBiogeochemistry) = (:B, ) @@ -54,7 +54,7 @@ end time_step!(model, 1) - @test all(particles.fields.A .== 0.1) + @test all(particles.fields.A .≈ 0.1) # also these shouldn't have moved @test all(particles.x .== 0) @@ -73,8 +73,8 @@ end time_step!(model, 1) - @test all(particles.fields.A .== 0.1) - @test all(particles.x .== 0.1) && all(particles.y .≈ 0.2) && all(particles.z .== 0) + @test all(particles.fields.A .≈ 0.1) + @test all(particles.x .≈ 0.1) && all(particles.y .≈ 0.2) && all(particles.z .≈ 0) end coupled_tracers(::SimpleParticleBiogeochemistry) = (:B, ) @@ -98,9 +98,9 @@ coupled_tracers(::SimpleParticleBiogeochemistry) = (:B, ) time_step!(model, 1) - @test all(particles.fields.A .== 0.1) + @test all(particles.fields.A .≈ 0.1) - @test interior(model.tracers.B, 1, 1, 3) .== 0.7 # 1 - 3 * 0.1 + @test interior(model.tracers.B, 1, 1, 3) .≈ 0.7 # 1 - 3 * 0.1 end @@ -123,7 +123,7 @@ end time_step!(model, 1) - @test all(particles.fields.A .== 0.1) + @test all(particles.fields.A .≈ 0.1) - @test interior(model.tracers.B, 1, 1, 3) .== 0.8 + @test interior(model.tracers.B, 1, 1, 3) .≈ 0.8 end \ No newline at end of file From 60182109627f779d0f81292b05fd07bc9ca095aa Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 14 Oct 2024 13:34:42 +0100 Subject: [PATCH 33/33] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 45edcade4..e73d523e8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright and contributors"] -version = "0.12.0" +version = "0.13.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"