From 158936c7f5429695b6ab6439d326f6bb6f1f4cb0 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 14:37:07 +0200 Subject: [PATCH 01/68] refactored models to abstract out light attenuation/sediments/particles more --- src/Boundaries/Sediments/Sediments.jl | 7 +- .../Sediments/coupled_timesteppers.jl | 7 +- src/Boundaries/gasexchange.jl | 2 +- src/Light/2band.jl | 4 +- src/Light/Light.jl | 5 +- src/Light/surface_PAR.jl | 0 .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 211 ++++++++---------- .../AdvectedPopulations/LOBSTER/fallbacks.jl | 24 +- src/Models/AdvectedPopulations/NPZD.jl | 110 ++++----- src/Models/Individuals/SLatissima.jl | 4 +- src/Models/Models.jl | 13 +- src/OceanBioME.jl | 74 +++++- src/Particles/Particles.jl | 10 +- 13 files changed, 231 insertions(+), 240 deletions(-) create mode 100644 src/Light/surface_PAR.jl diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 3487b8cad..ebe8268eb 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -3,7 +3,6 @@ module Sediments export SimpleMultiG, InstantRemineralisation using KernelAbstractions -using OceanBioME: ContinuousFormBiogeochemistry using Oceananigans using Oceananigans.Architectures: device using Oceananigans.Utils: launch! @@ -14,17 +13,17 @@ using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity using Oceananigans.Grids: zspacing using Oceananigans.Operators: volume using Oceananigans.Fields: Center +using OceanBioME: Biogeochemistry import Adapt: adapt_structure, adapt +import Oceananigans.Biogeochemistry: update_tendencies! abstract type AbstractSediment end abstract type FlatSediment <: AbstractSediment end sediment_fields(::AbstractSediment) = () -@inline update_sediment_tendencies!(bgc, sediment, model) = nothing - -function update_sediment_tendencies!(bgc, sediment::FlatSediment, model) +function update_tendencies!(bgc, sediment::FlatSediment, model) arch = model.grid.architecture for (i, tracer) in enumerate(sediment_tracers(sediment)) diff --git a/src/Boundaries/Sediments/coupled_timesteppers.jl b/src/Boundaries/Sediments/coupled_timesteppers.jl index 3942a7031..867785861 100644 --- a/src/Boundaries/Sediments/coupled_timesteppers.jl +++ b/src/Boundaries/Sediments/coupled_timesteppers.jl @@ -1,5 +1,4 @@ using Oceananigans: NonhydrostaticModel, prognostic_fields, HydrostaticFreeSurfaceModel -using OceanBioME: ContinuousFormBiogeochemistry using OceanBioME.Boundaries.Sediments: AbstractSediment using Oceananigans.TimeSteppers: ab2_step_field!, rk3_substep_field!, stage_Δt using Oceananigans.Utils: work_layout, launch! @@ -9,7 +8,7 @@ using Oceananigans.Architectures: AbstractArchitecture import Oceananigans.TimeSteppers: ab2_step!, rk3_substep! -@inline function ab2_step!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ) +@inline function ab2_step!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:FlatSediment}}, Δt, χ) workgroup, worksize = work_layout(model.grid, :xyz) arch = model.architecture step_field_kernel! = ab2_step_field!(device(arch), workgroup, worksize) @@ -46,7 +45,7 @@ import Oceananigans.TimeSteppers: ab2_step!, rk3_substep! return nothing end -@inline function ab2_step!(model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ) +@inline function ab2_step!(model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:FlatSediment}}, Δt, χ) # Step locally velocity and tracers @apply_regionally local_ab2_step!(model, Δt, χ) @@ -75,7 +74,7 @@ end @inbounds u[i, j, 1] += Δt * ((one_point_five + χ) * Gⁿ[i, j, 1] - (oh_point_five + χ) * G⁻[i, j, 1]) end -function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, γⁿ, ζⁿ) +function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:FlatSediment}}, Δt, γⁿ, ζⁿ) workgroup, worksize = work_layout(model.grid, :xyz) arch = model.architecture substep_field_kernel! = rk3_substep_field!(device(arch), workgroup, worksize) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index 8f18ab88b..badb7610a 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -236,7 +236,7 @@ function GasExchange(; gas, return FluxBoundaryCondition(gasexchange_function, field_dependencies = field_dependencies, parameters = gasexchange) end -# Hack this nicer format into Oceananigans boundary conditions because `typeof(gasexhcange) = DataType` and `BoundaryCondition` has a special use for `DataType` in the first argument so doesn't work +# Hack this nicer format into Oceananigans boundary conditions because `typeof(gasexchange) = DataType` and `BoundaryCondition` has a special use for `DataType` in the first argument so doesn't work @inline gasexchange_function(x, y, t, args...) = @inbounds args[end](x, y, t, args[1:end-1]...) @inline (gasexchange::GasExchange)(x, y, t, conc) = gasexchange(x, y, t, conc, gasexchange.temperature(x, y, 0.0, t), gasexchange.salinity(x, y, 0.0, t)) diff --git a/src/Light/2band.jl b/src/Light/2band.jl index f4b1f4a75..b966540c9 100644 --- a/src/Light/2band.jl +++ b/src/Light/2band.jl @@ -120,15 +120,13 @@ function TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR) end -function update_PAR!(model, PAR::TwoBandPhotosyntheticallyActiveRadiation) +function update_biogeochemical_state!(model, PAR::TwoBandPhotosyntheticallyActiveRadiation) arch = architecture(model.grid) launch!(arch, model.grid, :xy, update_TwoBandPhotosyntheticallyActiveRadiation!, PAR.field, model.grid, model.tracers.P, PAR.surface_PAR, model.clock.time, PAR) fill_halo_regions!(PAR.field, model.clock, fields(model)) end -required_PAR_fields(::TwoBandPhotosyntheticallyActiveRadiation) = (:PAR, ) - summary(::TwoBandPhotosyntheticallyActiveRadiation{FT}) where {FT} = string("Two-band light attenuation model ($FT)") show(io::IO, model::TwoBandPhotosyntheticallyActiveRadiation{FT}) where {FT} = print(io, summary(model)) diff --git a/src/Light/Light.jl b/src/Light/Light.jl index 2a561d3b3..fc871fc1a 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -22,13 +22,10 @@ using Oceananigans.Units import Adapt: adapt_structure, adapt import Base: show, summary -import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields +import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields, update_biogeochemical_state!, required_biogeochemical_auxiliary_fields import Oceananigans.BoundaryConditions: _fill_top_halo! # Fallback -update_PAR!(model, PAR) = nothing -required_PAR_fields(PAR) = () - include("2band.jl") include("morel.jl") diff --git a/src/Light/surface_PAR.jl b/src/Light/surface_PAR.jl new file mode 100644 index 000000000..e69de29bb diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 36cd9202e..08882f335 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -43,13 +43,11 @@ module LOBSTERModel export LOBSTER -using OceanBioME: ContinuousFormBiogeochemistry - using Oceananigans.Units using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField -using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, update_PAR!, required_PAR_fields -using OceanBioME: setup_velocity_fields, show_sinking_velocities +using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation +using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux @@ -57,9 +55,7 @@ import OceanBioME.BoxModels: update_boxmodel_state! import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - biogeochemical_drift_velocity, - update_biogeochemical_state!, - biogeochemical_auxiliary_fields + biogeochemical_drift_velocity import OceanBioME: maximum_sinking_velocity @@ -68,7 +64,7 @@ import Base: show, summary import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver -struct LOBSTER{FT, LA, S, B, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} +struct LOBSTER{FT, B, W} phytoplankton_preference :: FT maximum_grazing_rate :: FT grazing_half_saturation :: FT @@ -99,15 +95,10 @@ struct LOBSTER{FT, LA, S, B, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} disolved_organic_breakdown_rate :: FT zooplankton_calcite_dissolution :: FT - light_attenuation_model :: LA - sediment_model :: S - optionals :: B sinking_velocities :: W - particles :: P - function LOBSTER(phytoplankton_preference::FT, maximum_grazing_rate::FT, grazing_half_saturation::FT, @@ -138,53 +129,43 @@ struct LOBSTER{FT, LA, S, B, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} disolved_organic_breakdown_rate::FT, zooplankton_calcite_dissolution::FT, - light_attenuation_model::LA, - sediment_model::S, - optionals::B, - sinking_velocities::W, - - particles::P) where {FT, LA, S, B, W, P} - - return new{FT, LA, S, B, W, P}(phytoplankton_preference, - maximum_grazing_rate, - grazing_half_saturation, - light_half_saturation, - nitrate_ammonia_inhibition, - nitrate_half_saturation, - ammonia_half_saturation, - maximum_phytoplankton_growthrate, - zooplankton_assimilation_fraction, - zooplankton_mortality, - zooplankton_excretion_rate, - phytoplankton_mortality, - small_detritus_remineralisation_rate, - large_detritus_remineralisation_rate, - phytoplankton_exudation_fraction, - nitrifcaiton_rate, - ammonia_fraction_of_exudate, - ammonia_fraction_of_excriment, - ammonia_fraction_of_detritus, - phytoplankton_redfield, - organic_redfield, - phytoplankton_chlorophyll_ratio, - organic_carbon_calcate_ratio, - respiraiton_oxygen_nitrogen_ratio, - nitrifcation_oxygen_nitrogen_ratio, - slow_sinking_mortality_fraction, - fast_sinking_mortality_fraction, - disolved_organic_breakdown_rate, - zooplankton_calcite_dissolution, - - light_attenuation_model, - sediment_model, - - optionals, - - sinking_velocities, - - particles) + sinking_velocities::W) where {FT, B, W} + + return new{FT, B, W}(phytoplankton_preference, + maximum_grazing_rate, + grazing_half_saturation, + light_half_saturation, + nitrate_ammonia_inhibition, + nitrate_half_saturation, + ammonia_half_saturation, + maximum_phytoplankton_growthrate, + zooplankton_assimilation_fraction, + zooplankton_mortality, + zooplankton_excretion_rate, + phytoplankton_mortality, + small_detritus_remineralisation_rate, + large_detritus_remineralisation_rate, + phytoplankton_exudation_fraction, + nitrifcaiton_rate, + ammonia_fraction_of_exudate, + ammonia_fraction_of_excriment, + ammonia_fraction_of_detritus, + phytoplankton_redfield, + organic_redfield, + phytoplankton_chlorophyll_ratio, + organic_carbon_calcate_ratio, + respiraiton_oxygen_nitrogen_ratio, + nitrifcation_oxygen_nitrogen_ratio, + slow_sinking_mortality_fraction, + fast_sinking_mortality_fraction, + disolved_organic_breakdown_rate, + zooplankton_calcite_dissolution, + + optionals, + + sinking_velocities) end end @@ -325,57 +306,55 @@ function LOBSTER(; grid, optionals = Val((carbonates, oxygen, variable_redfield)) - return LOBSTER(phytoplankton_preference, - maximum_grazing_rate, - grazing_half_saturation, - light_half_saturation, - nitrate_ammonia_inhibition, - nitrate_half_saturation, - ammonia_half_saturation, - maximum_phytoplankton_growthrate, - zooplankton_assimilation_fraction, - zooplankton_mortality, - zooplankton_excretion_rate, - phytoplankton_mortality, - small_detritus_remineralisation_rate, - large_detritus_remineralisation_rate, - phytoplankton_exudation_fraction, - nitrifcaiton_rate, - ammonia_fraction_of_exudate, - ammonia_fraction_of_excriment, - ammonia_fraction_of_detritus, - phytoplankton_redfield, - organic_redfield, - phytoplankton_chlorophyll_ratio, - organic_carbon_calcate_ratio, - respiraiton_oxygen_nitrogen_ratio, - nitrifcation_oxygen_nitrogen_ratio, - slow_sinking_mortality_fraction, - fast_sinking_mortality_fraction, - disolved_organic_breakdown_rate, - zooplankton_calcite_dissolution, - - light_attenuation_model, - sediment_model, - - optionals, - - sinking_velocities, - - particles) + underlying_biogeochemistry = LOBSTER(phytoplankton_preference, + maximum_grazing_rate, + grazing_half_saturation, + light_half_saturation, + nitrate_ammonia_inhibition, + nitrate_half_saturation, + ammonia_half_saturation, + maximum_phytoplankton_growthrate, + zooplankton_assimilation_fraction, + zooplankton_mortality, + zooplankton_excretion_rate, + phytoplankton_mortality, + small_detritus_remineralisation_rate, + large_detritus_remineralisation_rate, + phytoplankton_exudation_fraction, + nitrifcaiton_rate, + ammonia_fraction_of_exudate, + ammonia_fraction_of_excriment, + ammonia_fraction_of_detritus, + phytoplankton_redfield, + organic_redfield, + phytoplankton_chlorophyll_ratio, + organic_carbon_calcate_ratio, + respiraiton_oxygen_nitrogen_ratio, + nitrifcation_oxygen_nitrogen_ratio, + slow_sinking_mortality_fraction, + fast_sinking_mortality_fraction, + disolved_organic_breakdown_rate, + zooplankton_calcite_dissolution, + + optionals, + + sinking_velocities) + + return Biogeochemistry(; underlying_biogeochemistry, + light_attenuation = light_attenuation_model, + sediment = sediment_model, + particles) end # wrote this functionally and it took 2.5x longer so even though this is long going to use this way instead -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, false, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :O₂) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, false, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk, :O₂) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :O₂, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :O₂, :sPOC, :bPOC, :DOC) - -@inline required_biogeochemical_auxiliary_fields(model::LOBSTER) = required_PAR_fields(model.light_attenuation_model) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, false, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, false, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, true, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :O₂) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, false, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk, :O₂) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :O₂, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :O₂, :sPOC, :bPOC, :DOC) const small_detritus = Union{Val{:sPON}, Val{:sPOC}} const large_detritus = Union{Val{:bPON}, Val{:bPOC}} @@ -397,10 +376,6 @@ const DOM = Union{Val{:DOM}, Val{:DON}} end end -function update_biogeochemical_state!(bgc::LOBSTER, model) - update_PAR!(model, bgc.light_attenuation_model) -end - function update_boxmodel_state!(model::BoxModel{<:LOBSTER, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) end @@ -441,11 +416,9 @@ adapt_structure(to, lobster::LOBSTER) = adapt(to, lobster.sinking_velocities), adapt(to, lobster.particles)) -summary(::LOBSTER{FT, LA, S, B, W, P}) where {FT, LA, S, B, W, P} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") -show(io::IO, model::LOBSTER{FT, LA, S, Val{B}, W, P}) where {FT, LA, S, B, W, P} = +summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") +show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = print(io, summary(model), " \n", - " Light Attenuation Model: ", "\n", - " └── ", summary(model.light_attenuation_model), "\n", " Optional components:", "\n", " ├── Carbonates $(B[1] ? :✅ : :❌) \n", " ├── Oxygen $(B[2] ? :✅ : :❌) \n", @@ -454,8 +427,6 @@ show(io::IO, model::LOBSTER{FT, LA, S, Val{B}, W, P}) where {FT, LA, S, B, W, P} @inline maximum_sinking_velocity(bgc::LOBSTER) = maximum(abs, bgc.sinking_velocities.bPOM.w) -@inline biogeochemical_auxiliary_fields(bgc::LOBSTER) = biogeochemical_auxiliary_fields(bgc.light_attenuation_model) - include("fallbacks.jl") include("core.jl") @@ -463,10 +434,10 @@ include("carbonate_chemistry.jl") include("oxygen_chemistry.jl") include("variable_redfield.jl") -const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, false, true)}, <:Any, <:Any}, - LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any}, - LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any}, - LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, true)}, <:Any, <:Any}} +const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true)}, <:Any}, + LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any}, + LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}, + LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}} @inline nitrogen_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = diff --git a/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl b/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl index bf9cc5d28..fd8902c26 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl @@ -12,22 +12,22 @@ # Carbonates and oxygen # We can't tell the difference between carbonates and oxygen, and variable redfields without specifying more # We're lucky that the optional groups have 2, 1, 3 extra variables each so this is the only non-unique conbination -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, O₂, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, O₂, PAR) # Carbonates and variable redfield -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) # Oxygen and variable redfield -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) # Carbonates, oxygen, and variable redfield @inline (bgc::LOBSTER)(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 3d4981b0a..7ab12fd1a 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -16,12 +16,12 @@ module NPZDModel export NutrientPhytoplanktonZooplanktonDetritus, NPZD -using OceanBioME: ContinuousFormBiogeochemistry +using OceanBioME: Biogeochemistry using Oceananigans.Units using Oceananigans.Fields: ZeroField -using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, update_PAR!, required_PAR_fields +using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation using OceanBioME: setup_velocity_fields, show_sinking_velocities using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux @@ -31,7 +31,6 @@ import Base: show, summary import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - update_biogeochemical_state!, biogeochemical_drift_velocity, biogeochemical_auxiliary_fields @@ -41,7 +40,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa import Adapt: adapt_structure, adapt -struct NutrientPhytoplanktonZooplanktonDetritus{FT, LA, S, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} +struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} # phytoplankton initial_photosynthetic_slope :: FT # α, 1/(W/m²)/s base_maximum_growth :: FT # μ₀, 1/s @@ -59,17 +58,9 @@ struct NutrientPhytoplanktonZooplanktonDetritus{FT, LA, S, W, P} <: ContinuousFo # detritus remineralization_rate :: FT # rᵈⁿ, 1/s - # light attenuation - light_attenuation_model :: LA - - # sediment - sediment_model :: S - # sinking sinking_velocities :: W - particles :: P - function NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope::FT, base_maximum_growth::FT, nutrient_half_saturation::FT, @@ -84,34 +75,22 @@ struct NutrientPhytoplanktonZooplanktonDetritus{FT, LA, S, W, P} <: ContinuousFo remineralization_rate::FT, - light_attenuation_model::LA, - - sediment_model::S, - - sinking_velocities::W, - - particles::P) where {FT, LA, S, W, P} - return new{FT, LA, S, W, P}(initial_photosynthetic_slope, - base_maximum_growth, - nutrient_half_saturation, - base_respiration_rate, - phyto_base_mortality_rate, - - maximum_grazing_rate, - grazing_half_saturation, - assimulation_efficiency, - base_excretion_rate, - zoo_base_mortality_rate, - - remineralization_rate, - - light_attenuation_model, - - sediment_model, - - sinking_velocities, - - particles) + sinking_velocities::W) where {FT, W} + return new{FT, W}(initial_photosynthetic_slope, + base_maximum_growth, + nutrient_half_saturation, + base_respiration_rate, + phyto_base_mortality_rate, + + maximum_grazing_rate, + grazing_half_saturation, + assimulation_efficiency, + base_excretion_rate, + zoo_base_mortality_rate, + + remineralization_rate, + + sinking_velocities) end end @@ -199,21 +178,24 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom) - return NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope, - base_maximum_growth, - nutrient_half_saturation, - base_respiration_rate, - phyto_base_mortality_rate, - maximum_grazing_rate, - grazing_half_saturation, - assimulation_efficiency, - base_excretion_rate, - zoo_base_mortality_rate, - remineralization_rate, - light_attenuation_model, - sediment_model, - sinking_velocities, - particles) + underlying_biogeochemistry = + NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope, + base_maximum_growth, + nutrient_half_saturation, + base_respiration_rate, + phyto_base_mortality_rate, + maximum_grazing_rate, + grazing_half_saturation, + assimulation_efficiency, + base_excretion_rate, + zoo_base_mortality_rate, + remineralization_rate, + sinking_velocities) + + return Biogeochemistry(; underlying_biogeochemistry, + light_attenuation = light_attenuation_model, + sediment = sediment_model, + particles) end const NPZD = NutrientPhytoplanktonZooplanktonDetritus @@ -298,24 +280,17 @@ end end end -function update_biogeochemical_state!(bgc::NPZD, model) - update_PAR!(model, bgc.light_attenuation_model) -end - function update_boxmodel_state!(model::BoxModel{<:NPZD, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) end -summary(::NPZD{FT, LA, W, P}) where {FT, LA, W, P} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") -show(io::IO, model::NPZD{FT, LA, W, P}) where {FT, LA, W, P} = +summary(::NPZD{FT, W}) where {FT, W} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") +show(io::IO, model::NPZD{FT, W}) where {FT, W} = print(io, summary(model), " \n", - " Light Attenuation Model: ", "\n", - " └── ", summary(model.light_attenuation_model), "\n", " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w) -@inline biogeochemical_auxiliary_fields(bgc::NPZD) = biogeochemical_auxiliary_fields(bgc.light_attenuation_model) adapt_structure(to, npzd::NPZD) = NutrientPhytoplanktonZooplanktonDetritus(npzd.initial_photosynthetic_slope, @@ -332,12 +307,7 @@ adapt_structure(to, npzd::NPZD) = npzd.remineralization_rate, - adapt(to, npzd.light_attenuation_model), - adapt(to, npzd.sediment_model), - - adapt(to, npzd.sinking_velocities), - - adapt(to, npzd.particles)) + adapt(to, npzd.sinking_velocities)) @inline nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 5b2a68ae2..1d8ead7c3 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -32,7 +32,7 @@ using Oceananigans.Grids: AbstractGrid using Oceananigans.Fields: fractional_indices, _interpolate, datatuple import Adapt: adapt_structure, adapt -import OceanBioME.Particles: update_particles_tendencies! +import Oceananigans.Biogeochemistry: update_tendencies! import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, _advect_particles! @inline no_dynamics(args...) = nothing @@ -264,7 +264,7 @@ adapt_structure(to, kelp::SLatissima) = SLatissima(kelp.architecture, kelp.scalefactor, kelp.latitude) -function update_particles_tendencies!(bgc, particles::SLatissima, model) +function update_tendencies!(bgc, particles::SLatissima, model) num_particles = length(particles) workgroup = min(num_particles, 256) worksize = num_particles diff --git a/src/Models/Models.jl b/src/Models/Models.jl index af590ffc1..c50245e38 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -1,14 +1,3 @@ -using OceanBioME.Boundaries.Sediments: update_sediment_tendencies! -using OceanBioME.Particles: update_particles_tendencies! - -import Oceananigans.Biogeochemistry: update_tendencies! - include("AdvectedPopulations/LOBSTER/LOBSTER.jl") include("AdvectedPopulations/NPZD.jl") -include("Individuals/SLatissima.jl") - -@inline function update_tendencies!(bgc::ContinuousFormBiogeochemistry, model) - update_sediment_tendencies!(bgc, bgc.sediment_model, model) - update_particles_tendencies!(bgc, bgc.particles, model) - return nothing -end +include("Individuals/SLatissima.jl") \ No newline at end of file diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 432dcf2a3..cf0529fec 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -17,7 +17,7 @@ export BoxModel, BoxModelGrid, SaveBoxModel, run!, set! export Particles # Light models -export TwoBandPhotosyntheticallyActiveRadiation, update_PAR! +export TwoBandPhotosyntheticallyActiveRadiation # Boundaries export Boundaries, Sediments, GasExchange, FlatSediment @@ -32,8 +32,76 @@ export zero_negative_tracers!, error_on_neg!, warn_on_neg!, ScaleNegativeTracers export ColumnField, isacolumn using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry - -abstract type ContinuousFormBiogeochemistry{LA, S, P} <: AbstractContinuousFormBiogeochemistry end +using Adapt + +import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, + required_biogeochemical_auxiliary_fields, + update_biogeochemical_state!, + biogeochemical_drift_velocity, + biogeochemical_auxiliary_fields, + update_tendencies! + +import Adapt: adapt_structure + +struct Biogeochemistry{B, L, S, P, I} <: AbstractContinuousFormBiogeochemistry + underlying_biogeochemistry :: B + light_attenuation :: L + sediment :: S + particles :: P + inputs :: I + + Biogeochemistry(underlying_biogeochemistry::B, + light_attenuation::L, + sediment::S, + particles::P, + inputs::I) where {B, L, S, P, I} = + new{B, L, S, P, I}(underlying_biogeochemistry, + light_attenuation, + sediment, + particles, + inputs) +end + +Biogeochemistry(; underlying_biogeochemistry, + light_attenuation = nothing, + sediment = nothing, + particles = nothing, + inputs = nothing) = + Biogeochemistry(underlying_biogeochemistry, + light_attenuation, + sediment, + particles, + inputs) + +required_biogeochemical_tracers(bgc::Biogeochemistry) = required_biogeochemical_tracers(bgc.underlying_biogeochemistry) + +required_biogeochemical_auxiliary_fields(bgc::Biogeochemistry) = required_biogeochemical_auxiliary_fields(bgc.underlying_biogeochemistry) + +biogeochemical_drift_velocity(bgc::Biogeochemistry, val_name) = biogeochemical_drift_velocity(bgc.underlying_biogeochemistry, val_name) + +biogeochemical_auxiliary_fields(bgc::Biogeochemistry) = merge(biogeochemical_auxiliary_fields(bgc.underlying_biogeochemistry), + biogeochemical_auxiliary_fields(bgc.light_attenuation)) + +adapt_structure(to, bgc::Biogeochemistry) = Biogeochemistry(adapt(to, bgc.underlying_biogeochemistry), + adapt(to, bgc.light_attenuation), + adapt(to, bgc.sediment), + adapt(to, bgc.particles), + adapt(to, bgc.inputs)) + +@inline function update_tendencies!(bgc::Biogeochemistry, model) + update_tendencies!(bgc, bgc.sediment, model) + update_tendencies!(bgc, bgc.particles, model) + #update_input_tendencies!(bgc, bgc.inputs, model) + return nothing +end + +@inline (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) + +function update_biogeochemical_state!(bgc::Biogeochemistry, model) + update_biogeochemical_state!(model, bgc.light_attenuation) +end + +update_tendencies!(bgc, ::Nothing, model) = nothing struct BoxModelGrid end diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 8972bdfa7..a9f28472b 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -1,8 +1,9 @@ module Particles -using OceanBioME: ContinuousFormBiogeochemistry using Oceananigans: NonhydrostaticModel, HydrostaticFreeSurfaceModel +using OceanBioME: Biogeochemistry +import Oceananigans.Biogeochemistry: update_tendencies! import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, step_lagrangian_particles! import Oceananigans.OutputWriters: fetch_output import Base: length, size, show, summary @@ -13,20 +14,19 @@ abstract type BiogeochemicalParticles end @inline step_lagrangian_particles!(::Nothing, model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, - <:ContinuousFormBiogeochemistry{<:Any, <:Any, <:BiogeochemicalParticles}}, + <:Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}}, Δ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, - <:ContinuousFormBiogeochemistry{<:Any, <:Any, <:BiogeochemicalParticles}, + <:Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, <:Any, <:Any, <:Any, <:Any, <:Any,}, Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) -@inline update_lagrangian_particle_properties!(model, bgc::ContinuousFormBiogeochemistry{<:Any, <:Any, <:BiogeochemicalParticles}, Δt) = +@inline update_lagrangian_particle_properties!(model, bgc::Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, Δt) = update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) update_lagrangian_particle_properties!(::BiogeochemicalParticles, model, bgc, Δt) = nothing -update_particles_tendencies!(bgc, particles, model) = nothing size(particles::BiogeochemicalParticles) = size(particles.x) length(particles::BiogeochemicalParticles) = length(particles.x) From b4d7e9c65379fedcefc1b7b7af9b67ef337a3a62 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 14:42:22 +0200 Subject: [PATCH 02/68] fixed --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 4 ++-- src/Models/AdvectedPopulations/NPZD.jl | 4 ++-- src/OceanBioME.jl | 4 ++++ 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 08882f335..07b77a39d 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -47,7 +47,7 @@ using Oceananigans.Units using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation -using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry +using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, UnderlyingBiogeochemicalModel using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux @@ -64,7 +64,7 @@ import Base: show, summary import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver -struct LOBSTER{FT, B, W} +struct LOBSTER{FT, B, W} <: UnderlyingBiogeochemicalModel phytoplankton_preference :: FT maximum_grazing_rate :: FT grazing_half_saturation :: FT diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 7ab12fd1a..160d83168 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -16,7 +16,7 @@ module NPZDModel export NutrientPhytoplanktonZooplanktonDetritus, NPZD -using OceanBioME: Biogeochemistry +using OceanBioME: Biogeochemistry, UnderlyingBiogeochemicalModel using Oceananigans.Units using Oceananigans.Fields: ZeroField @@ -40,7 +40,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa import Adapt: adapt_structure, adapt -struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} +struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} <: UnderlyingBiogeochemicalModel # phytoplankton initial_photosynthetic_slope :: FT # α, 1/(W/m²)/s base_maximum_growth :: FT # μ₀, 1/s diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index cf0529fec..fc4489365 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -103,6 +103,10 @@ end update_tendencies!(bgc, ::Nothing, model) = nothing +abstract type UnderlyingBiogeochemicalModel end + +@inline (::UnderlyingBiogeochemicalModel)(val_tracer_name, x, y, z, t, fields...) = zero(x) + struct BoxModelGrid end @inline maximum_sinking_velocity(bgc) = 0.0 From f353bcc526b55c0222bba7fa618e093b0fc7de5c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 14:44:34 +0200 Subject: [PATCH 03/68] didn't mean to include that --- src/Light/surface_PAR.jl | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/Light/surface_PAR.jl diff --git a/src/Light/surface_PAR.jl b/src/Light/surface_PAR.jl deleted file mode 100644 index e69de29bb..000000000 From d8c520020552f464f8808cab7fbbc07ad3c231a8 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 14:54:36 +0200 Subject: [PATCH 04/68] fixed LOBSTER --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 07b77a39d..671705e2a 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -356,6 +356,8 @@ end @inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :O₂, :sPOC, :bPOC, :DOC) @inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :O₂, :sPOC, :bPOC, :DOC) +required_biogeochemical_auxiliary_fields(::LOBSTER) = (:PAR, ) + const small_detritus = Union{Val{:sPON}, Val{:sPOC}} const large_detritus = Union{Val{:bPON}, Val{:bPOC}} const disolved_organic_matter = Union{Val{:DON}, Val{:DOC}} From cb6cdfaf822427ba7972d2a1a2bf784c37fcc586 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 15:09:33 +0200 Subject: [PATCH 05/68] updated light attenuation test --- test/test_light.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_light.jl b/test/test_light.jl index 7be4c12af..38a500bbc 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -1,6 +1,6 @@ using Oceananigans, Test using OceanBioME: TwoBandPhotosyntheticallyActiveRadiation, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus -using Oceananigans.Biogeochemistry: update_biogeochemical_state!, required_biogeochemical_tracers +using Oceananigans.Biogeochemistry: update_biogeochemical_state!, required_biogeochemical_tracers, biogeochemical_auxiliary_fields Pᵢ(x,y,z) = 2.5 + z @@ -37,7 +37,7 @@ function test_two_band(grid, bgc, model_type) expected_PAR = 100.0 .* [exp(- 0.5 * kʳ - ∫Chlʳ[1] * χʳ) + exp(- 0.5 * kᵇ - ∫Chlᵇ[1] * χᵇ), exp(- 1.5 * kʳ - ∫Chlʳ[2] * χʳ) + exp(- 1.5 * kᵇ - ∫Chlᵇ[2] * χᵇ)] ./ 2 - results_PAR = convert(Array, model.biogeochemistry.light_attenuation_model.field)[1, 1, 1:2] + results_PAR = convert(Array, biogeochemical_auxiliary_fields(biogeochemistry).PAR)[1, 1, 1:2] return all(results_PAR .≈ reverse(expected_PAR)) end From 9a8fc8e54ef8d8614969261ecffa485b0e28c170 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 15:55:23 +0200 Subject: [PATCH 06/68] updated light attenuation test --- test/test_light.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_light.jl b/test/test_light.jl index 38a500bbc..14d3e582f 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -17,7 +17,7 @@ function test_two_band(grid, bgc, model_type) update_biogeochemical_state!(model.biogeochemistry, model) - PAR_model = model.biogeochemistry.light_attenuation_model + PAR_model = model.biogeochemistry.light_attenuation kʳ = PAR_model.water_red_attenuation kᵇ = PAR_model.water_blue_attenuation From f8feb4e120c16935a53bcbe21250bfd5c25a0c03 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:08:54 +0200 Subject: [PATCH 07/68] fixed test lobster --- test/test_LOBSTER.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index c3ef4eeca..1f8bcd792 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -11,7 +11,7 @@ function ΣC(model, carbonates, variable_redfield) if variable_redfield OC = sum(model.tracers.sPOC .+ model.tracers.bPOC .+ model.tracers.DOC) else - OC = sum(model.tracers.sPOM .+ model.tracers.bPOM .+ model.tracers.DOM) * model.biogeochemistry. organic_redfield + OC = sum(model.tracers.sPOM .+ model.tracers.bPOM .+ model.tracers.DOM) * model.biogeochemistry.underlying_model.organic_redfield end if carbonates @@ -20,7 +20,7 @@ function ΣC(model, carbonates, variable_redfield) IC = 0.0 end - LC = sum(model.tracers.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.phytoplankton_redfield + LC = sum(model.tracers.P * (1 + model.biogeochemistry.underlying_model.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_model.phytoplankton_redfield return OC + IC + LC end @@ -90,9 +90,9 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open model.tracers.bPON .= rand() model.tracers.DON .= rand() - model.tracers.sPOC .= model.tracers.sPON * model.biogeochemistry.phytoplankton_redfield - model.tracers.bPOC .= model.tracers.bPON * model.biogeochemistry.phytoplankton_redfield - model.tracers.DOC .= model.tracers.DON * model.biogeochemistry.phytoplankton_redfield + model.tracers.sPOC .= model.tracers.sPON * model.biogeochemistry.underlying_model.phytoplankton_redfield + model.tracers.bPOC .= model.tracers.bPON * model.biogeochemistry.underlying_model.phytoplankton_redfield + model.tracers.DOC .= model.tracers.DON * model.biogeochemistry.underlying_model.phytoplankton_redfield else model.tracers.sPOM .= rand() model.tracers.bPOM .= rand() From c4999288a0c32b8497f2d4397b8837949c3ce9b8 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:14:00 +0200 Subject: [PATCH 08/68] fixed LOBSTER adapt --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 671705e2a..132c0b3f6 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -412,11 +412,8 @@ adapt_structure(to, lobster::LOBSTER) = lobster.fast_sinking_mortality_fraction, lobster.disolved_organic_breakdown_rate, lobster.zooplankton_calcite_dissolution, - adapt(to, lobster.light_attenuation_model), - adapt(to, lobster.sediment_model), lobster.optionals, - adapt(to, lobster.sinking_velocities), - adapt(to, lobster.particles)) + adapt(to, lobster.sinking_velocities)) summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = From c19bb631df2a41fbc02834e3880a686e87c17a08 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:15:04 +0200 Subject: [PATCH 09/68] fixed examples --- examples/column.jl | 8 +++++--- examples/data_forced.jl | 6 ++++-- examples/kelp.jl | 4 ++-- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/examples/column.jl b/examples/column.jl index b28f72f62..4ffdd344c 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -102,10 +102,12 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +redfield = model.biogeochemistry.underlying_biogeochemistry.organic_redfield + for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20]) * model.biogeochemistry.organic_redfield + air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) + carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 0e8774a43..14116a5d1 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -135,10 +135,12 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +redfield = model.biogeochemistry.underlying_biogeochemistry.organic_redfield + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20]) * model.biogeochemistry.organic_redfield + carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/kelp.jl b/examples/kelp.jl index e8e3eea2b..e38a4c478 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -132,8 +132,8 @@ carbon_export = zeros(length(times)) for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + - bPOC[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20] + carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + + bPOC[1, 1, grid.Nz-20, i] * 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)`. From 5d7d0193cc5becc6b19db87a5c9a9db3dc9fa58f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:34:07 +0200 Subject: [PATCH 10/68] added redfield function --- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 14 +++++++++++++- src/Models/AdvectedPopulations/NPZD.jl | 9 +++++++-- src/OceanBioME.jl | 19 ++++++++++++++++++- 3 files changed, 38 insertions(+), 4 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 132c0b3f6..dae6ea723 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -51,6 +51,7 @@ using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistr using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux +import OceanBioME: redfield import OceanBioME.BoxModels: update_boxmodel_state! import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, @@ -438,12 +439,23 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}, LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}} +@inline redfield(i, j, k, val_tracer_name, bgc::LOBSTER, tracers) = redfield(val_tracer_name, bgc) + +@inline redfield(::Val{:P}, bgc::LOBSTER) = (1 + bgc.organic_carbon_calcate_ratio) * bgc.phytoplankton_redfield +@inline redfield(::Val{:Z}, bgc::LOBSTER) = bgc.phytoplankton_redfield +@inline redfield(::Union{Val{:NO₃}, Val{:NH₄}, Val{:Alk}, Val{:O₂}}, bgc::LOBSTER) = 0 +@inline redfield(::Union{Val{:sPOM}, Val{:bPOM}, Val{:DOM}}, bgc::LOBSTER) = bgc.organic_redfield +@inline redfield(::Union{Val{:sPOC}, Val{:bPOC}, Val{:DOC}, Val{:DIC}}, bgc::LOBSTER) = 1 + +@inline redfield(i, j, k, ::Val{:sPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.sPOC[i, j, k] / tracers.sPON[i, j, k] +@inline redfield(i, j, k, ::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.bPOC[i, j, k] / tracers.bPON[i, j, k] +@inline redfield(i, j, k, ::Val{:DON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.DOC[i, j, k] / tracers.DON[i, j, k] @inline nitrogen_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:sPOM), bgc, tracers) + sinking_flux(i, j, grid, advection, Val(:bPOM), bgc, tracers) -@inline carbon_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = nitrogen_flux(grid, advection, bgc, tracers, i, j) * bgc.organic_redfield +@inline carbon_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = nitrogen_flux(grid, advection, bgc, tracers, i, j) * redfield(Val(:sPOM), bgc) @inline nitrogen_flux(grid, advection, bgc::lobster_variable_redfield, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:sPON), bgc, tracers) + diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 160d83168..53eb04910 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -26,6 +26,7 @@ using OceanBioME: setup_velocity_fields, show_sinking_velocities using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux +import OceanBioME: redfield import OceanBioME.BoxModels: update_boxmodel_state! import Base: show, summary @@ -309,9 +310,13 @@ adapt_structure(to, npzd::NPZD) = adapt(to, npzd.sinking_velocities)) +@inline redfield(i, j, k, val_tracer_name, bgc::NPZD, tracers) = redfield(val_tracer_name, bgc) +@inline redfield(::Union{Val{:N}}, bgc::NPZD) = 0 +@inline redfield(::Union{Val{:P}, Val{:Z}, Val{:D}}, bgc::NPZD) = 6.56 + @inline nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) - -@inline carbon_flux(bgc::NPZD, tracers, i, j) = nitrogen_flux(bgc, tracers, i, j) * 6.56 + +@inline carbon_flux(bgc::NPZD, tracers, i, j) = nitrogen_flux(bgc, tracers, i, j) * redfield(Val(:P), bgc, tracers) @inline remineralisation_receiver(::NPZD) = :N end # module diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index fc4489365..8f2e3df5a 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -5,7 +5,7 @@ between ocean biogeochemistry, carbonate chemistry, and physics. module OceanBioME # Biogeochemistry models -export LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, PISCES, NPZD +export LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, PISCES, NPZD, redfield # Macroalgae models export SLatissima @@ -111,6 +111,23 @@ struct BoxModelGrid end @inline maximum_sinking_velocity(bgc) = 0.0 +""" + redfield(i, j, k, val_tracer_name, bgc, tracers) + +Returns the redfield ratio of `tracer_name` from `bgc` at `i`, `j`, `k`. +""" +@inline redfield(i, j, k, val_tracer_name, bgc, tracers) = NaN + +""" + redfield(val_tracer_name, bgc, tracers) + +Returns the redfield ratio of `tracer_name` from `bgc` when it is constant across the domain. +""" +@inline redfield(val_tracer_name, bgc) = NaN + +@inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) +@inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) + include("Utils/Utils.jl") include("Boundaries/Boundaries.jl") include("Light/Light.jl") From cb01e55c720a1371f953107ce44b33118d6fd6da Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:36:58 +0200 Subject: [PATCH 11/68] updated examples to use redfield function --- examples/column.jl | 4 +--- examples/data_forced.jl | 6 ++---- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/examples/column.jl b/examples/column.jl index 4ffdd344c..27349b7e3 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -102,12 +102,10 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) -redfield = model.biogeochemistry.underlying_biogeochemistry.organic_redfield - for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 14116a5d1..4b5cec8c0 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -17,7 +17,7 @@ # ## Model setup # First load the required packages -using OceanBioME +using OceanBioME using Oceananigans, Random, Printf, NetCDF, Interpolations, DataDeps using Oceananigans.Units @@ -135,12 +135,10 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) -redfield = model.biogeochemistry.underlying_biogeochemistry.organic_redfield - for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. From d8cb2059146c5e976205aea4280f1be670117497 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:39:59 +0200 Subject: [PATCH 12/68] fix LOBSTER test --- test/test_LOBSTER.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index 1f8bcd792..a1121b671 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -11,7 +11,7 @@ function ΣC(model, carbonates, variable_redfield) if variable_redfield OC = sum(model.tracers.sPOC .+ model.tracers.bPOC .+ model.tracers.DOC) else - OC = sum(model.tracers.sPOM .+ model.tracers.bPOM .+ model.tracers.DOM) * model.biogeochemistry.underlying_model.organic_redfield + OC = sum(model.tracers.sPOM .+ model.tracers.bPOM .+ model.tracers.DOM) * model.biogeochemistry.underlying_biogeochemistry.organic_redfield end if carbonates @@ -20,7 +20,7 @@ function ΣC(model, carbonates, variable_redfield) IC = 0.0 end - LC = sum(model.tracers.P * (1 + model.biogeochemistry.underlying_model.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_model.phytoplankton_redfield + LC = sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield return OC + IC + LC end @@ -90,9 +90,9 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open model.tracers.bPON .= rand() model.tracers.DON .= rand() - model.tracers.sPOC .= model.tracers.sPON * model.biogeochemistry.underlying_model.phytoplankton_redfield - model.tracers.bPOC .= model.tracers.bPON * model.biogeochemistry.underlying_model.phytoplankton_redfield - model.tracers.DOC .= model.tracers.DON * model.biogeochemistry.underlying_model.phytoplankton_redfield + model.tracers.sPOC .= model.tracers.sPON * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield + model.tracers.bPOC .= model.tracers.bPON * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield + model.tracers.DOC .= model.tracers.DON * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield else model.tracers.sPOM .= rand() model.tracers.bPOM .= rand() From daa7e9e5b0d9110700a6c9d0e4a13b12bf58f7ed Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:44:23 +0200 Subject: [PATCH 13/68] fixed column example --- examples/column.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/column.jl b/examples/column.jl index 27349b7e3..68e9774c5 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -103,7 +103,7 @@ air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) + air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), (args...) -> 35) carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) end From a15377a2692eecb614e05190e246556ffb5298de Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:53:03 +0200 Subject: [PATCH 14/68] added some docstrings and moved `underlying_model` to arg (rather than kwarg --- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 8 ++--- src/Models/AdvectedPopulations/NPZD.jl | 8 ++--- src/OceanBioME.jl | 33 ++++++++++++++----- 3 files changed, 33 insertions(+), 16 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index dae6ea723..b28be47ac 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -341,10 +341,10 @@ function LOBSTER(; grid, sinking_velocities) - return Biogeochemistry(; underlying_biogeochemistry, - light_attenuation = light_attenuation_model, - sediment = sediment_model, - particles) + return Biogeochemistry(underlying_biogeochemistry; + light_attenuation = light_attenuation_model, + sediment = sediment_model, + particles) end # wrote this functionally and it took 2.5x longer so even though this is long going to use this way instead diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 53eb04910..31ae9b050 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -193,10 +193,10 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, remineralization_rate, sinking_velocities) - return Biogeochemistry(; underlying_biogeochemistry, - light_attenuation = light_attenuation_model, - sediment = sediment_model, - particles) + return Biogeochemistry(underlying_biogeochemistry; + light_attenuation = light_attenuation_model, + sediment = sediment_model, + particles) end const NPZD = NutrientPhytoplanktonZooplanktonDetritus diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 8f2e3df5a..01c2e6b90 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -62,11 +62,29 @@ struct Biogeochemistry{B, L, S, P, I} <: AbstractContinuousFormBiogeochemistry inputs) end -Biogeochemistry(; underlying_biogeochemistry, - light_attenuation = nothing, - sediment = nothing, - particles = nothing, - inputs = nothing) = +""" + Biogeochemistry(underlying_biogeochemistry; + light_attenuation = nothing, + sediment = nothing, + particles = nothing, + inputs = nothing) + +Construct a biogeochemical model based on `underlying_biogeochemistry` which may have +a `light_attenuation` model, `sediment`, `particles`, and `inputs`. + +Keyword Arguments +================= + +- `light_attenuation_model`: light attenuation model which integrated the attenuation of available light +- `sediment_model`: slot for `AbstractSediment` +- `particles`: slot for `BiogeochemicalParticles` +- `inputs`: slot for nutrient inputs such as rivers (work in progress) +""" +Biogeochemistry(underlying_biogeochemistry, + light_attenuation = nothing, + sediment = nothing, + particles = nothing, + inputs = nothing) = Biogeochemistry(underlying_biogeochemistry, light_attenuation, sediment, @@ -88,11 +106,10 @@ adapt_structure(to, bgc::Biogeochemistry) = Biogeochemistry(adapt(to, bgc.underl adapt(to, bgc.particles), adapt(to, bgc.inputs)) -@inline function update_tendencies!(bgc::Biogeochemistry, model) +function update_tendencies!(bgc::Biogeochemistry, model) update_tendencies!(bgc, bgc.sediment, model) update_tendencies!(bgc, bgc.particles, model) - #update_input_tendencies!(bgc, bgc.inputs, model) - return nothing + update_tendencies!(bgc, bgc.inputs, model) end @inline (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) From da3307791a8785158071d33d89f2f2c3f2168d40 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 19:57:40 +0200 Subject: [PATCH 15/68] fixed stuff and added show/summary --- src/OceanBioME.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 01c2e6b90..f4ab1ce4f 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -42,6 +42,7 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, update_tendencies! import Adapt: adapt_structure +import Base: show, summary struct Biogeochemistry{B, L, S, P, I} <: AbstractContinuousFormBiogeochemistry underlying_biogeochemistry :: B @@ -80,7 +81,7 @@ Keyword Arguments - `particles`: slot for `BiogeochemicalParticles` - `inputs`: slot for nutrient inputs such as rivers (work in progress) """ -Biogeochemistry(underlying_biogeochemistry, +Biogeochemistry(underlying_biogeochemistry; light_attenuation = nothing, sediment = nothing, particles = nothing, @@ -145,6 +146,13 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros @inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) @inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) +summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") +show(io::IO, model::Biogeochemistry) = + print(io, show(model.underlying_biogeochemistry), " \n", + " Light attenuation: $(summary(model.light_attenuation))", "\n", + " Sediment: $(summary(model.sediment))", "\n", + " Particles: $(summary(model.particles))") + include("Utils/Utils.jl") include("Boundaries/Boundaries.jl") include("Light/Light.jl") From 0b2def66f6ed511b9069177aeef494b441805992 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 20:15:24 +0200 Subject: [PATCH 16/68] fixed doctests (hopefully again) --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 14 +++++++++----- src/Models/AdvectedPopulations/NPZD.jl | 14 ++++++++------ src/OceanBioME.jl | 6 +++--- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index b28be47ac..dd05994ff 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -245,15 +245,16 @@ julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); julia> model = LOBSTER(; grid) Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model (Float64) - Light Attenuation Model: - └── Two-band light attenuation model (Float64) Optional components: ├── Carbonates ❌ ├── Oxygen ❌ └── Variable Redfield Ratio ❌ Sinking Velocities: ├── sPOM: 0.0 to -3.47e-5 m/s - └── bPOM: 0.0 to -0.0023148148148148147 m/s + └── bPOM: 0.0 to -0.0023148148148148147 m/s + Light attenuation: Two-band light attenuation model (Float64) + Sediment: Nothing + Particles: Nothing ``` """ function LOBSTER(; grid, @@ -417,14 +418,17 @@ adapt_structure(to, lobster::LOBSTER) = adapt(to, lobster.sinking_velocities)) summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") -show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = - print(io, summary(model), " \n", + +show(model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = + string(summary(model), " \n", " Optional components:", "\n", " ├── Carbonates $(B[1] ? :✅ : :❌) \n", " ├── Oxygen $(B[2] ? :✅ : :❌) \n", " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) +show(io::IO, model::LOBSTER) = print(io, show(model)) + @inline maximum_sinking_velocity(bgc::LOBSTER) = maximum(abs, bgc.sinking_velocities.bPOM.w) include("fallbacks.jl") diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 31ae9b050..c3ae789c3 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -146,11 +146,13 @@ julia> grid = RectilinearGrid(size=(20, 30), extent=(200, 200), topology=(Bounde julia> model = NutrientPhytoplanktonZooplanktonDetritus(; grid) Nutrient Phytoplankton Zooplankton Detritus model (Float64) - Light Attenuation Model: - └── Two-band light attenuation model (Float64) Sinking Velocities: ├── P: 0.0 to -2.9525462962962963e-6 m/s - └── D: 0.0 to -3.181597222222222e-5 m/s + └── D: 0.0 to -3.181597222222222e-5 m/s + Light attenuation: Two-band light attenuation model (Float64) + Sediment: Nothing + Particles: Nothing + ``` """ function NutrientPhytoplanktonZooplanktonDetritus(; grid, @@ -287,9 +289,9 @@ function update_boxmodel_state!(model::BoxModel{<:NPZD, <:Any, <:Any, <:Any, <:A end summary(::NPZD{FT, W}) where {FT, W} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") -show(io::IO, model::NPZD{FT, W}) where {FT, W} = - print(io, summary(model), " \n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) +show(model::NPZD) = string(summary(model), " \n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) +show(io::IO, model::NPZD) = print(io, show(model)) @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index f4ab1ce4f..4abbe0feb 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -149,9 +149,9 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") show(io::IO, model::Biogeochemistry) = print(io, show(model.underlying_biogeochemistry), " \n", - " Light attenuation: $(summary(model.light_attenuation))", "\n", - " Sediment: $(summary(model.sediment))", "\n", - " Particles: $(summary(model.particles))") + " Light attenuation: ", summary(model.light_attenuation), "\n", + " Sediment: ", summary(model.sediment), "\n", + " Particles: ", summary(model.particles)) include("Utils/Utils.jl") include("Boundaries/Boundaries.jl") From b88d638113db06353518a0ba6873f1d19a861cec Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 20:35:14 +0200 Subject: [PATCH 17/68] Fixed other LOBSTER test function (oops) --- test/test_LOBSTER.jl | 4 ++-- test/test_light.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index a1121b671..ed667b850 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -35,7 +35,7 @@ function ΣGᶜ(model, carbonates, variable_redfield) if variable_redfield OC = sum(model.timestepper.Gⁿ.sPOC .+ model.timestepper.Gⁿ.bPOC .+ model.timestepper.Gⁿ.DOC) else - OC = sum(model.timestepper.Gⁿ.sPOM .+ model.timestepper.Gⁿ.bPOM .+ model.timestepper.Gⁿ.DOM) * model.biogeochemistry. organic_redfield + OC = sum(model.timestepper.Gⁿ.sPOM .+ model.timestepper.Gⁿ.bPOM .+ model.timestepper.Gⁿ.DOM) * model.biogeochemistry.underlying_biogeochemistry.organic_redfield end if carbonates @@ -44,7 +44,7 @@ function ΣGᶜ(model, carbonates, variable_redfield) IC = 0.0 end - LC = sum(model.timestepper.Gⁿ.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.timestepper.Gⁿ.Z) * model.biogeochemistry.phytoplankton_redfield + LC = sum(model.timestepper.Gⁿ.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.timestepper.Gⁿ.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield return OC + IC + LC end diff --git a/test/test_light.jl b/test/test_light.jl index 14d3e582f..fdcf99169 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -49,7 +49,7 @@ archs = (CPU(), ) arch in archs, grid in (RectilinearGrid(arch; size = (2, 2, 2), extent = (2, 2, 2)), LatitudeLongitudeGrid(arch; size = (5, 5, 2), longitude = (-180, 180), latitude = (-85, 85), z = (-2, 0))), - bgc in (LOBSTER, NutrientPhytoplanktonZooplanktonDetritus) + bgc in (LOBSTER, NutrientPhytoplanktonZooplanktonDetritus) # this is now redundant since each model doesn't deal with the light separatly if !((model == NonhydrostaticModel) && ((grid isa LatitudeLongitudeGrid) | (grid isa OrthogonalSphericalShellGrid))) @info "Testing $bgc in $model on $grid..." From e0c93f1f5950b4994a8306834eae82e6d88a6bb8 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 20:37:21 +0200 Subject: [PATCH 18/68] Fixed example --- examples/column.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/column.jl b/examples/column.jl index 68e9774c5..14571a6f0 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -103,7 +103,7 @@ air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), (args...) -> 35) + air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) end From def4adf99bf58322df45f5b52a9b2fb4f16395ff Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 20:43:05 +0200 Subject: [PATCH 19/68] fixed examples (again) --- examples/column.jl | 4 ++-- examples/data_forced.jl | 6 ++++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/column.jl b/examples/column.jl index 14571a6f0..1453467ca 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -104,8 +104,8 @@ carbon_export = zeros(length(times)) for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) + carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 4b5cec8c0..4370b591b 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -135,10 +135,12 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) + carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. From 749d35902fdb6e363e58a43ac93426e27474d7b1 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 20:46:52 +0200 Subject: [PATCH 20/68] bumped compats --- Manifest.toml | 155 ++++++++++++++++++++++++++++---------------------- Project.toml | 8 +-- 2 files changed, 90 insertions(+), 73 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 94e39d56e..bcc3f050c 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,18 +1,19 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.2" +julia_version = "1.9.1" manifest_format = "2.0" -project_hash = "1df4fb4dcc221fff4ae2f940474ce7322fc6228d" +project_hash = "b8336a615981d764881d88f0d8193dd7f4667beb" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] -git-tree-sha1 = "cad4c758c0038eea30394b1b671526921ca85b21" +git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef" uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" -version = "1.4.0" -weakdeps = ["ChainRulesCore"] +version = "1.5.0" +weakdeps = ["ChainRulesCore", "Test"] [deps.AbstractFFTs.extensions] AbstractFFTsChainRulesCoreExt = "ChainRulesCore" + AbstractFFTsTestExt = "Test" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] @@ -87,9 +88,9 @@ version = "0.1.2" [[deps.CUDA]] deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Preferences", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "35160ef0f03b14768abfd68b830f8e3940e8e0dc" +git-tree-sha1 = "968c1365e2992824c3e7a794e30907483f8469a9" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "4.4.0" +version = "4.4.1" [[deps.CUDA_Driver_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] @@ -128,9 +129,9 @@ version = "0.2.4" [[deps.Compat]] deps = ["UUIDs"] -git-tree-sha1 = "4e88377ae7ebeaf29a047aa1ee40826e0b708a5d" +git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.7.0" +version = "4.9.0" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -139,13 +140,13 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+0" +version = "1.0.2+0" [[deps.ConstructionBase]] deps = ["LinearAlgebra"] -git-tree-sha1 = "738fec4d684a9a6ee9598a8bfee305b26831f28c" +git-tree-sha1 = "fe2838a593b5f776e1597e086dcd47560d94e816" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.2" +version = "1.5.3" [deps.ConstructionBase.extensions] ConstructionBaseIntervalSetsExt = "IntervalSets" @@ -173,9 +174,9 @@ version = "1.15.0" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "cf25ccb972fec4e4817764d01c82386ae94f77b4" +git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.14" +version = "0.18.15" [[deps.DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" @@ -186,6 +187,16 @@ version = "1.0.0" deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[deps.Distances]] +deps = ["LinearAlgebra", "Statistics", "StatsAPI"] +git-tree-sha1 = "b6def76ffad15143924a2199f72a5cd883a2e8a9" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.9" +weakdeps = ["SparseArrays"] + + [deps.Distances.extensions] + DistancesSparseArraysExt = "SparseArrays" + [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -207,9 +218,9 @@ uuid = "b305315f-e792-5b7a-8f41-49f472929428" version = "1.0.1" [[deps.ExprTools]] -git-tree-sha1 = "c1d06d129da9f55715c6c212866f5b1bddc5fa00" +git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" -version = "0.1.9" +version = "0.1.10" [[deps.FFTW]] deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] @@ -261,9 +272,9 @@ version = "1.3.1" [[deps.HDF5_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"] -git-tree-sha1 = "3b20c3ce9c14aedd0adca2bc8c882927844bd53d" +git-tree-sha1 = "592e1c427983a465831fc73c5ae0ca5d0ac13a9e" uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" -version = "1.14.0+0" +version = "1.14.1+0" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -278,9 +289,9 @@ version = "0.2.1" [[deps.IntelOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "0cb9352ef2e01574eeebdb102948a58740dcaf83" +git-tree-sha1 = "ad37c091f7d7daf900963171600d7c1c5c3ede32" uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2023.1.0+0" +version = "2023.2.0+0" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -304,27 +315,27 @@ version = "1.0.0" [[deps.JLD2]] deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"] -git-tree-sha1 = "5df8278ad24772c0c6dbbeb97b162ccf29ced2a9" +git-tree-sha1 = "aa6ffef1fd85657f4999030c52eaeec22a279738" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.32" +version = "0.4.33" [[deps.JLLWrappers]] -deps = ["Preferences"] -git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.4.1" +version = "1.5.0" [[deps.JSON3]] deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] -git-tree-sha1 = "5b62d93f2582b09e469b3099d839c2d2ebf5066d" +git-tree-sha1 = "95220473901735a0f4df9d1ca5b171b568b2daa3" uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -version = "1.13.1" +version = "1.13.2" [[deps.KernelAbstractions]] deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "6d08ca80b621635fed9cdfeb9a4280a574020bf3" +git-tree-sha1 = "4c5875e4c228247e1c2b087669846941fb6e0118" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.7" +version = "0.9.8" [deps.KernelAbstractions.extensions] EnzymeExt = "EnzymeCore" @@ -388,9 +399,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "c3ce8e7420b3a6e071e0fe4745f5d4300e37b13f" +git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.24" +version = "0.3.26" [deps.LogExpFunctions.extensions] LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -407,15 +418,15 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] -git-tree-sha1 = "154d7aaa82d24db6d8f7e4ffcfe596f40bff214b" +git-tree-sha1 = "eb006abbd7041c28e0d16260e50a24f8f9104913" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2023.1.0+0" +version = "2023.2.0+0" [[deps.MPI]] deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"] -git-tree-sha1 = "98fc280a56d3b5cc218435f82df60e05771fefa6" +git-tree-sha1 = "32cafbe56c7f0b7160a1a6c492773af66c0b722f" uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195" -version = "0.20.12" +version = "0.20.14" [deps.MPI.extensions] AMDGPUExt = "AMDGPU" @@ -433,9 +444,9 @@ version = "4.1.2+0" [[deps.MPIPreferences]] deps = ["Libdl", "Preferences"] -git-tree-sha1 = "d86a788b336e8ae96429c0c42740ccd60ac0dfcc" +git-tree-sha1 = "781916a2ebf2841467cda03b6f1af43e23839d85" uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" -version = "0.1.8" +version = "0.1.9" [[deps.MPItrampoline_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] @@ -445,9 +456,9 @@ version = "5.3.1+0" [[deps.MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.10" +version = "0.5.11" [[deps.Markdown]] deps = ["Base64"] @@ -488,10 +499,10 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" [[deps.Oceananigans]] -deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] -git-tree-sha1 = "c5120ad3b0d67ab23e9401063ecb9c14999f2613" +deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] +git-tree-sha1 = "996eb159d0a1380d09beb7190ed0807d51ba0924" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.85.0" +version = "0.87.1" [[deps.OffsetArrays]] deps = ["Adapt"] @@ -517,9 +528,9 @@ version = "4.1.5+0" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "cae3153c7f6cf3f069a853883fd1919a6e5bab5b" +git-tree-sha1 = "e78db7bd5c26fc5a6911b50a47ee302219157ea8" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.9+0" +version = "3.0.10+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -528,21 +539,21 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" [[deps.OrderedCollections]] -git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" +git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.0" +version = "1.6.2" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "4b2e829ee66d4218e0cef22c0a64ee37cf258c29" +git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.7.1" +version = "2.7.2" [[deps.PencilArrays]] deps = ["Adapt", "JSON3", "LinearAlgebra", "MPI", "OffsetArrays", "Random", "Reexport", "StaticArrayInterface", "StaticArrays", "StaticPermutations", "Strided", "TimerOutputs", "VersionParsing"] -git-tree-sha1 = "c1d2b659ce423897de45e998f49f77e789e1859d" +git-tree-sha1 = "c30a7fb1e424ea572962dac493ad6c3f556695e0" uuid = "0e08944d-e94e-41b1-9406-dcf66b6a9d2e" -version = "0.18.1" +version = "0.19.0" [deps.PencilArrays.extensions] PencilArraysDiffEqExt = ["DiffEqBase"] @@ -554,14 +565,14 @@ version = "0.18.1" [[deps.PencilFFTs]] deps = ["AbstractFFTs", "FFTW", "LinearAlgebra", "MPI", "PencilArrays", "Reexport", "TimerOutputs"] -git-tree-sha1 = "af200cef52069f3428127b237919d7c69eb6b10d" +git-tree-sha1 = "bd69f3f0ee248cfb4241800aefb705b5ded592ff" uuid = "4a48f351-57a6-4416-9ec4-c37015456aae" -version = "0.15.0" +version = "0.15.1" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" +version = "1.9.0" [[deps.PkgVersion]] deps = ["Pkg"] @@ -571,9 +582,9 @@ version = "0.3.2" [[deps.PrecompileTools]] deps = ["Preferences"] -git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81" +git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f" uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.1.2" +version = "1.2.0" [[deps.Preferences]] deps = ["TOML"] @@ -587,9 +598,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[deps.ProgressBars]] deps = ["Printf"] -git-tree-sha1 = "9d84c8646109eb8bc7a006d59b157c64d5155c81" +git-tree-sha1 = "b437cdb0385ed38312d91d9c00c20f3798b30256" uuid = "49802e3a-d2f1-5c88-81d8-b72133a6f568" -version = "1.5.0" +version = "1.5.1" [[deps.Quaternions]] deps = ["LinearAlgebra", "Random", "RealDot"] @@ -642,9 +653,9 @@ version = "1.3.0" [[deps.Roots]] deps = ["ChainRulesCore", "CommonSolve", "Printf", "Setfield"] -git-tree-sha1 = "de432823e8aab4dd1a985be4be768f95acf152d4" +git-tree-sha1 = "ff42754a57bb0d6dcfe302fd0d4272853190421f" uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -version = "2.0.17" +version = "2.0.19" [deps.Roots.extensions] RootsForwardDiffExt = "ForwardDiff" @@ -701,9 +712,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "7beb031cf8145577fbccacd94b8a8f4ce78428d3" +git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.3.0" +version = "2.3.1" weakdeps = ["ChainRulesCore"] [deps.SpecialFunctions.extensions] @@ -711,9 +722,9 @@ weakdeps = ["ChainRulesCore"] [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "dbde6766fc677423598138a5951269432b0fcc90" +git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.7" +version = "0.8.8" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "Static", "SuiteSparse"] @@ -728,18 +739,18 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "Random", "StaticArraysCore"] -git-tree-sha1 = "fffc14c695c17bfdbfa92a2a01836cdc542a1e46" +git-tree-sha1 = "9cabadf6e7cd2349b6cf49f1915ad2028d65e881" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.6.1" +version = "1.6.2" weakdeps = ["Statistics"] [deps.StaticArrays.extensions] StaticArraysStatisticsExt = "Statistics" [[deps.StaticArraysCore]] -git-tree-sha1 = "1d5708d926c76a505052d0d24a846d5da08bc3a4" +git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.1" +version = "1.4.2" [[deps.StaticPermutations]] git-tree-sha1 = "193c3daa18ff3e55c1dae66acb6a762c4a3bdb0b" @@ -751,11 +762,17 @@ deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" version = "1.9.0" +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.6.0" + [[deps.Strided]] deps = ["LinearAlgebra", "StridedViews", "TupleTools"] -git-tree-sha1 = "b32eadf6ac726a790567fdc872b63117712e16a8" +git-tree-sha1 = "137303f5e0a39f966b462c53ae2c5c6e34c4828b" uuid = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" -version = "2.0.1" +version = "2.0.3" [[deps.StridedViews]] deps = ["LinearAlgebra"] diff --git a/Project.toml b/Project.toml index ebbdf412a..c2cb9b761 100644 --- a/Project.toml +++ b/Project.toml @@ -12,12 +12,12 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" [compat] -Adapt = "3.4" -JLD2 = "0.4.31" +Adapt = "3" +JLD2 = "^0.4" KernelAbstractions = "0.9" -Oceananigans = "0.84.1, 0.85" +Oceananigans = "0.84.1, 0.85, 0.86, 0.87" Roots = "2" -StructArrays = "0.6.15" +StructArrays = "0.4, 0.5, 0.6" julia = "1.9" [extras] From 8d2914b03aae7966318c629830104c134958b7ec Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 21:31:07 +0200 Subject: [PATCH 21/68] fixed kelp test --- test/test_slatissima.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index d9738887a..ada407cd9 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -40,7 +40,7 @@ end initial_kelp_N = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.N .+ particles.structural_nitrogen)) ./ (14 * 0.001) initial_tracer_C = sum(model.tracers.sPOC) + sum(model.tracers.bPOC) + sum(model.tracers.DOC) + sum(model.tracers.DIC) + - sum(model.tracers.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.phytoplankton_redfield + sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield initial_kelp_C = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.C .+ particles.structural_carbon)) ./ (12 * 0.001) @@ -54,7 +54,7 @@ end final_kelp_N = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.N .+ particles.structural_nitrogen)) ./ (14 * 0.001) final_tracer_C = sum(model.tracers.sPOC) + sum(model.tracers.bPOC) + sum(model.tracers.DOC) + sum(model.tracers.DIC) + - sum(model.tracers.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.phytoplankton_redfield + sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield final_kelp_C = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.C .+ particles.structural_carbon)) ./ (12 * 0.001) # kelp is being integrated From 285c82b213a1ac279975fea75278e88a98b57d9f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 21:31:52 +0200 Subject: [PATCH 22/68] fixed kelp example --- examples/kelp.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/kelp.jl b/examples/kelp.jl index e38a4c478..3a10a155c 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -130,6 +130,8 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + From b688c68b69f43d81b6accb1b401b3b201cdee98c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 21:49:03 +0200 Subject: [PATCH 23/68] oops hadn't fixed column example --- examples/column.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/column.jl b/examples/column.jl index 1453467ca..201490814 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -102,10 +102,12 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), biogeochemistry) + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. From 9ea977e0afa1bf385029cec5f1f9aa6d3f6542e5 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 22:35:45 +0200 Subject: [PATCH 24/68] fix sediment tests --- test/test_sediments.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 2d6ebede1..d43e84c03 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -103,8 +103,8 @@ display_name(::InstantRemineralisation) = "Instant remineralisation" sediment_model)) if !(isa(sediment_model, SimpleMultiG) && isa(biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)) - @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry))" - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry))" test_flat_sediment(grid, biogeochemistry; timestepper) + @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry; timestepper) end end end From 1b5693bc349c1e2db2ffe815262a03d83109e14c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 23:49:25 +0200 Subject: [PATCH 25/68] fix sediment tests (again) --- test/test_sediments.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index d43e84c03..47db2ee9e 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -51,7 +51,7 @@ function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBash closure = nothing, timestepper) - set_defaults!(model.biogeochemistry.sediment_model) + set_defaults!(model.biogeochemistry.sediment) set_defaults!(biogeochemistry, model) @@ -61,7 +61,7 @@ function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBash simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) - N₀ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment_model) * Azᶠᶜᶜ(1, 1, 1, grid) + N₀ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) run!(simulation) @@ -69,14 +69,14 @@ function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBash @test any([any(intercepted_tendencies[tracer] .!= model.timestepper.Gⁿ[tracer]) for tracer in keys(model.tracers)]) # the sediment tendencies are being updated - @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment_model.tendencies.Gⁿ]) - @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment_model.tendencies.G⁻]) + @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.Gⁿ]) + @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.G⁻]) # the sediment values are being integrated initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) - @test all([any(field .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment_model.fields)]) + @test all([any(field .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) - N₁ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment_model) * Azᶠᶜᶜ(1, 1, 1, grid) + N₁ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) # conservations @test N₁ ≈ N₀ From b6530fd7cfc5c077738178f8e78a7671954204cd Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 27 Aug 2023 23:58:11 +0200 Subject: [PATCH 26/68] not sure why column example is being more unstable now --- examples/data_forced.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 4370b591b..19e368283 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -105,7 +105,7 @@ scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄ simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, - max_change = 2.0, min_change = 0.5, + max_change = 1.5, min_change = 0.75, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) From 5dd1832f69e8409657c04919d1997cc3d37bd68d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 28 Aug 2023 00:42:19 +0200 Subject: [PATCH 27/68] fixed sediment tests again again --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 47db2ee9e..ce221687f 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -53,7 +53,7 @@ function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBash set_defaults!(model.biogeochemistry.sediment) - set_defaults!(biogeochemistry, model) + set_defaults!(biogeochemistry.underlying_biogeochemistry, model) simulation = Simulation(model, Δt = 50, stop_time = 1day) From 080488bbdffbf466b9cd81b2f534467ffa5fc279 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 28 Aug 2023 00:44:12 +0200 Subject: [PATCH 28/68] doctest fix --- docs/src/model_components/biogeochemical/LOBSTER.md | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/src/model_components/biogeochemical/LOBSTER.md b/docs/src/model_components/biogeochemical/LOBSTER.md index d5686452b..848a57c73 100644 --- a/docs/src/model_components/biogeochemical/LOBSTER.md +++ b/docs/src/model_components/biogeochemical/LOBSTER.md @@ -13,15 +13,16 @@ julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); julia> bgc_model = LOBSTER(; grid, carbonates = true) Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model (Float64) - Light Attenuation Model: - └── Two-band light attenuation model (Float64) Optional components: ├── Carbonates ✅ ├── Oxygen ❌ └── Variable Redfield Ratio ❌ Sinking Velocities: ├── sPOM: 0.0 to -3.47e-5 m/s - └── bPOM: 0.0 to -0.0023148148148148147 m/s + └── bPOM: 0.0 to -0.0023148148148148147 m/s + Light attenuation: Two-band light attenuation model (Float64) + Sediment: Nothing + Particles: Nothing ``` ## Model equations From b10a4ec68907bfb6f1f1e0bf6bc0556d1c169762 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 11:18:34 +0200 Subject: [PATCH 29/68] fixed (and locally tested) sediment --- src/Boundaries/Sediments/Sediments.jl | 8 ++++---- src/Boundaries/Sediments/coupled_timesteppers.jl | 2 +- src/Boundaries/Sediments/instant_remineralization.jl | 3 +++ test/test_sediments.jl | 4 ++-- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index ebe8268eb..b0cb521ba 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -32,7 +32,7 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper) + sediment, bgc, model.grid, model.advection, model.tracers, model.timestepper) return nothing end @@ -43,9 +43,9 @@ end end -@inline nitrogen_flux(grid, adveciton, bgc, tracers, i, j) = 0 -@inline carbon_flux(grid, adveciton, bgc, tracers, i, j) = 0 -@inline remineralisation_receiver(bgc) = nothing +@inline nitrogen_flux(grid, adveciton, bgc::Biogeochemistry, tracers, i, j) = nitrogen_flux(grid, adveciton, bgc.underlying_biogeochemistry, tracers, i, j) +@inline carbon_flux(grid, adveciton, bgc::Biogeochemistry, tracers, i, j) = carbon_flux(grid, adveciton, bgc.underlying_biogeochemistry, tracers, i, j) +@inline remineralisation_receiver(bgc::Biogeochemistry) = remineralisation_receiver(bgc.underlying_biogeochemistry) @inline sinking_flux(i, j, grid, advection, val_tracer::Val{T}, bgc, tracers) where T = - advective_tracer_flux_z(i, j, 1, grid, advection, biogeochemical_drift_velocity(bgc, val_tracer).w, tracers[T]) / diff --git a/src/Boundaries/Sediments/coupled_timesteppers.jl b/src/Boundaries/Sediments/coupled_timesteppers.jl index 867785861..e32c52dab 100644 --- a/src/Boundaries/Sediments/coupled_timesteppers.jl +++ b/src/Boundaries/Sediments/coupled_timesteppers.jl @@ -32,7 +32,7 @@ import Oceananigans.TimeSteppers: ab2_step!, rk3_substep! Δt) end - sediment = model.biogeochemistry.sediment_model + sediment = model.biogeochemistry.sediment for (i, field) in enumerate(sediment_fields(sediment)) launch!(arch, model.grid, :xy, ab2_step_flat_field!, diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 2c08caad3..d3ed03bcb 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -91,3 +91,6 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st @inbounds timestepper.Gⁿ[remineralisation_receiver(bgc)][i, j, 1] += flux * (1 - burial_efficiency) / Δz end + +summary(::InstantRemineralisation) = string("Single-layer instant remineralisaiton ($FT)") +show(io::IO, model::InstantRemineralisation) = print(io, summary(model)) \ No newline at end of file diff --git a/test/test_sediments.jl b/test/test_sediments.jl index ce221687f..e60eaf7f4 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -61,7 +61,7 @@ function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBash simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) - N₀ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + N₀ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) run!(simulation) @@ -76,7 +76,7 @@ function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBash initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) @test all([any(field .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) - N₁ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + N₁ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) # conservations @test N₁ ≈ N₀ From 3479ba85526e1aa5fac54992ea345a940d0866c9 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 12:21:38 +0200 Subject: [PATCH 30/68] fixed rk3 --- src/Boundaries/Sediments/coupled_timesteppers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Boundaries/Sediments/coupled_timesteppers.jl b/src/Boundaries/Sediments/coupled_timesteppers.jl index e32c52dab..61c2f085d 100644 --- a/src/Boundaries/Sediments/coupled_timesteppers.jl +++ b/src/Boundaries/Sediments/coupled_timesteppers.jl @@ -52,7 +52,7 @@ end # blocking step for implicit free surface, non blocking for explicit ab2_step_free_surface!(model.free_surface, model, Δt, χ) - sediment = model.biogeochemistry.sediment_model + sediment = model.biogeochemistry.sediment arch = model.architecture for (i, field) in enumerate(sediment_fields(sediment)) @@ -97,7 +97,7 @@ function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:A stage_Δt(Δt, γⁿ, ζⁿ)) end - sediment = model.biogeochemistry.sediment_model + sediment = model.biogeochemistry.sediment for (i, field) in enumerate(sediment_fields(sediment)) launch!(arch, model.grid, :xy, rk3_step_flat_field!, From 991d480f5d85bbd99fc398e53ee6deb339051206 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 14:49:29 +0200 Subject: [PATCH 31/68] fixed NPZD cabon flux --- src/Models/AdvectedPopulations/NPZD.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index c3ae789c3..7a7b50789 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -319,6 +319,6 @@ adapt_structure(to, npzd::NPZD) = @inline nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(bgc::NPZD, tracers, i, j) = nitrogen_flux(bgc, tracers, i, j) * redfield(Val(:P), bgc, tracers) +@inline carbon_flux(grid, advection, bgc::NPZD, tracers, i, j) = nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) * redfield(Val(:P), bgc, tracers) @inline remineralisation_receiver(::NPZD) = :N end # module From 660f6a7c7bb505256770c3aff169d20a1bab9a5e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 16:27:27 +0200 Subject: [PATCH 32/68] again --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index e60eaf7f4..b00bae378 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -102,7 +102,7 @@ display_name(::InstantRemineralisation) = "Instant remineralisation" open_bottom = true, sediment_model)) - if !(isa(sediment_model, SimpleMultiG) && isa(biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)) + if !(isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)) @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry; timestepper) end From 63d062ba0d632558d53d8f3d500e6ac580046c2c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 16:28:32 +0200 Subject: [PATCH 33/68] fixed NPZD redfield --- src/Models/AdvectedPopulations/NPZD.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 7a7b50789..d7d9d49dd 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -319,6 +319,6 @@ adapt_structure(to, npzd::NPZD) = @inline nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(grid, advection, bgc::NPZD, tracers, i, j) = nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) * redfield(Val(:P), bgc, tracers) +@inline carbon_flux(grid, advection, bgc::NPZD, tracers, i, j) = nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) * redfield(Val(:P), bgc) @inline remineralisation_receiver(::NPZD) = :N end # module From 172192335f682ca665f6533f39e754cb7ae36c56 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 17:42:18 +0200 Subject: [PATCH 34/68] fixed show stuff --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 15 ++++++--------- src/Models/AdvectedPopulations/NPZD.jl | 5 ++--- src/OceanBioME.jl | 2 +- 3 files changed, 9 insertions(+), 13 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index dd05994ff..e1003fb84 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -419,15 +419,12 @@ adapt_structure(to, lobster::LOBSTER) = summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") -show(model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = - string(summary(model), " \n", - " Optional components:", "\n", - " ├── Carbonates $(B[1] ? :✅ : :❌) \n", - " ├── Oxygen $(B[2] ? :✅ : :❌) \n", - " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) - -show(io::IO, model::LOBSTER) = print(io, show(model)) +show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = print(io, string(summary(model), " \n", + " Optional components:", "\n", + " ├── Carbonates $(B[1] ? :✅ : :❌) \n", + " ├── Oxygen $(B[2] ? :✅ : :❌) \n", + " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities))) @inline maximum_sinking_velocity(bgc::LOBSTER) = maximum(abs, bgc.sinking_velocities.bPOM.w) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index d7d9d49dd..e3b8733c6 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -289,9 +289,8 @@ function update_boxmodel_state!(model::BoxModel{<:NPZD, <:Any, <:Any, <:Any, <:A end summary(::NPZD{FT, W}) where {FT, W} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") -show(model::NPZD) = string(summary(model), " \n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) -show(io::IO, model::NPZD) = print(io, show(model)) +show(io::IO, model::NPZD) = print(io, string(summary(model), " \n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities))) @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 4abbe0feb..a1ef6a2a4 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -148,7 +148,7 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") show(io::IO, model::Biogeochemistry) = - print(io, show(model.underlying_biogeochemistry), " \n", + print(io, show(io, model.underlying_biogeochemistry), " \n", " Light attenuation: ", summary(model.light_attenuation), "\n", " Sediment: ", summary(model.sediment), "\n", " Particles: ", summary(model.particles)) From 7f064242da8734b99031eaa3b2977fa64009fe37 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 29 Aug 2023 20:31:46 +0200 Subject: [PATCH 35/68] fixed docstrings --- src/Boundaries/Sediments/instant_remineralization.jl | 2 +- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 12 ++++++------ src/Models/AdvectedPopulations/NPZD.jl | 4 ++-- src/OceanBioME.jl | 2 +- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index d3ed03bcb..e7881f089 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -92,5 +92,5 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st @inbounds timestepper.Gⁿ[remineralisation_receiver(bgc)][i, j, 1] += flux * (1 - burial_efficiency) / Δz end -summary(::InstantRemineralisation) = string("Single-layer instant remineralisaiton ($FT)") +summary(::InstantRemineralisation{FT}) where {FT} = string("Single-layer instant remineralisaiton ($FT)") show(io::IO, model::InstantRemineralisation) = print(io, summary(model)) \ No newline at end of file diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index e1003fb84..915fc9dd1 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -419,12 +419,12 @@ adapt_structure(to, lobster::LOBSTER) = summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") -show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = print(io, string(summary(model), " \n", - " Optional components:", "\n", - " ├── Carbonates $(B[1] ? :✅ : :❌) \n", - " ├── Oxygen $(B[2] ? :✅ : :❌) \n", - " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities))) +show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = string(summary(model), " \n", + " Optional components:", "\n", + " ├── Carbonates $(B[1] ? :✅ : :❌) \n", + " ├── Oxygen $(B[2] ? :✅ : :❌) \n", + " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) @inline maximum_sinking_velocity(bgc::LOBSTER) = maximum(abs, bgc.sinking_velocities.bPOM.w) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index e3b8733c6..1f18f00ed 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -289,8 +289,8 @@ function update_boxmodel_state!(model::BoxModel{<:NPZD, <:Any, <:Any, <:Any, <:A end summary(::NPZD{FT, W}) where {FT, W} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") -show(io::IO, model::NPZD) = print(io, string(summary(model), " \n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities))) +show(io::IO, model::NPZD) = string(summary(model), " \n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index a1ef6a2a4..4abbe0feb 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -148,7 +148,7 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") show(io::IO, model::Biogeochemistry) = - print(io, show(io, model.underlying_biogeochemistry), " \n", + print(io, show(model.underlying_biogeochemistry), " \n", " Light attenuation: ", summary(model.light_attenuation), "\n", " Sediment: ", summary(model.sediment), "\n", " Particles: ", summary(model.particles)) From 47610ede2a1613c637bd2cf235862d59359fbc80 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 09:45:52 +0200 Subject: [PATCH 36/68] attempt to update model implimentation page --- docs/src/model_implementation.md | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 1d3d864ac..47c588d58 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -18,20 +18,18 @@ For this example we are going to implement the simple Nutrient-Phytoplankton mod The first step is to import the abstract type from OceanBioME, some units from Oceananigans (for ease of parameter definition), and [`import`](https://stackoverflow.com/questions/27086159/what-is-the-difference-between-using-and-import-in-julia-when-building-a-mod) some functions from Oceananigans in order to add methods to: ```@example implementing -using OceanBioME: ContinuousFormBiogeochemistry +using OceanBioME: Biogeochemistry using Oceananigans.Units import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - update_biogeochemical_state!, - biogeochemical_drift_velocity, - biogeochemical_auxiliary_fields + biogeochemical_drift_velocity ``` We then define our `struct` with the model parameters, as well as slots for the particles, light attenuation, and sediment models: ```@example implementing -@kwdef struct NutrientPhytoplankton{FT, LA, S, W, P} <:ContinuousFormBiogeochemistry{LA, S, P} +@kwdef struct NutrientPhytoplankton{FT, W} base_growth_rate :: FT = 1.27 / day # 1 / seconds nutrient_half_saturation :: FT = 0.025 * 1000 / 14 # mmol N / m³ light_half_saturation :: FT = 300.0 # micro einstein / m² / s @@ -41,10 +39,7 @@ We then define our `struct` with the model parameters, as well as slots for the mortality_rate :: FT = 0.15 / day # 1 / seconds crowding_mortality_rate :: FT = 0.004 / day / 1000 * 14 # 1 / seconds / mmol N / m³ - light_attenuation_model :: LA = nothing - sediment_model :: S = nothing sinking_velocity :: W = 2 / day - particles :: P = nothing end ``` @@ -54,8 +49,6 @@ Here, we use descriptive names for the parameters. Below, each of these paramete required_biogeochemical_tracers(::NutrientPhytoplankton) = (:N, :P, :T) required_biogeochemical_auxiliary_fields(::NutrientPhytoplankton) = (:PAR, ) - -biogeochemical_auxiliary_fields(bgc::NutrientPhytoplankton) = biogeochemical_auxiliary_fields(bgc.light_attenuation_model) ``` Next, we define the functions that specify how the phytoplankton ``P`` evolve. In the absence of advection and diffusion (both of which are handled by Oceananigans), we want the phytoplankton to evolve at the rate given by: @@ -126,10 +119,6 @@ Finally, we need to define some functions to allow us to update the time-depende using OceanBioME: BoxModel import OceanBioME.BoxModels: update_boxmodel_state! -function update_biogeochemical_state!(bgc::NutrientPhytoplankton, model) - update_PAR!(model, bgc.light_attenuation_model) -end - function update_boxmodel_state!(model::BoxModel{<:NutrientPhytoplankton, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) @@ -228,7 +217,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa @inline nitrogen_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(bgc::NutrientPhytoplankton, tracers, i, j) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 +@inline carbon_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 @inline remineralisation_receiver(::NutrientPhytoplankton) = :N ``` From 5db1ac2bb3bbdb9bafd577c0e1857acf50722d9d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 09:53:53 +0200 Subject: [PATCH 37/68] fixed missed merge conflicts --- src/Boundaries/Sediments/Sediments.jl | 16 +++++----------- src/Models/AdvectedPopulations/NPZD.jl | 12 ++---------- 2 files changed, 7 insertions(+), 21 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index e941d2105..aff6af773 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -3,26 +3,20 @@ module Sediments export SimpleMultiG, InstantRemineralisation using KernelAbstractions -<<<<<<< HEAD -======= -using OceanBioME: ContinuousFormBiogeochemistry, BoxModelGrid ->>>>>>> origin + +using OceanBioME: Biogeochemistry, BoxModelGrid + using Oceananigans using Oceananigans.Architectures: device using Oceananigans.Utils: launch! using Oceananigans.Advection: advective_tracer_flux_z using Oceananigans.Units: day -using Oceananigans.Fields: CenterField, Face +using Oceananigans.Fields: CenterField, Face, Center, ConstantField using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity using Oceananigans.Grids: zspacing using Oceananigans.Operators: volume -<<<<<<< HEAD -using Oceananigans.Fields: Center -using OceanBioME: Biogeochemistry -======= -using Oceananigans.Fields: Center, ConstantField using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, immersed_cell ->>>>>>> origin + import Adapt: adapt_structure, adapt import Oceananigans.Biogeochemistry: update_tendencies! diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 1962e5128..5b05b3410 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -32,8 +32,7 @@ import Base: show, summary import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - biogeochemical_drift_velocity, - biogeochemical_auxiliary_fields + biogeochemical_drift_velocity import OceanBioME: maximum_sinking_velocity @@ -315,17 +314,10 @@ adapt_structure(to, npzd::NPZD) = @inline redfield(::Union{Val{:N}}, bgc::NPZD) = 0 @inline redfield(::Union{Val{:P}, Val{:Z}, Val{:D}}, bgc::NPZD) = 6.56 -<<<<<<< HEAD -@inline nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:D), bgc, tracers) + - sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) - -@inline carbon_flux(grid, advection, bgc::NPZD, tracers, i, j) = nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) * redfield(Val(:P), bgc) -======= @inline nitrogen_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = sinking_flux(i, j, k, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 +@inline carbon_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * redfield(Val(:P), bgc) ->>>>>>> origin @inline remineralisation_receiver(::NPZD) = :N end # module From 6a1110aa15daad981e07559acc75aaa747e991bf Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 30 Aug 2023 12:44:11 +0200 Subject: [PATCH 38/68] Update test_sediments.jl --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 65b6f497c..73c7f940a 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -135,7 +135,7 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) if run @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry; timestepper) + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry, model; timestepper) end end end From ea4926855b5e6b95e68b0e68aaa874c91af312b3 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 30 Aug 2023 12:48:09 +0200 Subject: [PATCH 39/68] Update model_implementation.md --- docs/src/model_implementation.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 47c588d58..a46d0a390 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -243,9 +243,9 @@ grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = (32, ), extent = # setup the biogeochemical model -light_attenuation_model = TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR) +light_attenuation = TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR) -sediment_model = InstantRemineralisation(; grid) +sediment = InstantRemineralisation(; grid) sinking_velocity = ZFaceField(grid) @@ -253,7 +253,9 @@ w_sink(x, y, z) = 2 / day * tanh(z / 5) set!(sinking_velocity, w_sink) -biogeochemistry = NutrientPhytoplankton(; light_attenuation_model, sinking_velocity, sediment_model) +biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; light_attenuation_model, sinking_velocity); + light_attenuation, + sediment) # put the model together From 4b7593a2758339873066847e81f07a0ff100a1a8 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 16:57:23 +0200 Subject: [PATCH 40/68] fixed implimentation docs --- docs/src/model_implementation.md | 8 ++++---- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 2 +- src/Models/AdvectedPopulations/NPZD.jl | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index a46d0a390..03b4968e8 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -18,7 +18,7 @@ For this example we are going to implement the simple Nutrient-Phytoplankton mod The first step is to import the abstract type from OceanBioME, some units from Oceananigans (for ease of parameter definition), and [`import`](https://stackoverflow.com/questions/27086159/what-is-the-difference-between-using-and-import-in-julia-when-building-a-mod) some functions from Oceananigans in order to add methods to: ```@example implementing -using OceanBioME: Biogeochemistry +using OceanBioME: Biogeochemistry, UnderlyingBiogeochemicalModel using Oceananigans.Units import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, @@ -29,7 +29,7 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, We then define our `struct` with the model parameters, as well as slots for the particles, light attenuation, and sediment models: ```@example implementing -@kwdef struct NutrientPhytoplankton{FT, W} +@kwdef struct NutrientPhytoplankton{FT, W} <: UnderlyingBiogeochemicalModel base_growth_rate :: FT = 1.27 / day # 1 / seconds nutrient_half_saturation :: FT = 0.025 * 1000 / 14 # mmol N / m³ light_half_saturation :: FT = 300.0 # micro einstein / m² / s @@ -253,7 +253,7 @@ w_sink(x, y, z) = 2 / day * tanh(z / 5) set!(sinking_velocity, w_sink) -biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; light_attenuation_model, sinking_velocity); +biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; sinking_velocity); light_attenuation, sediment) @@ -275,7 +275,7 @@ simulation.output_writers[:tracers] = JLD2OutputWriter(model, model.tracers, schedule = TimeInterval(1day), overwrite_existing = true) -simulation.output_writers[:sediment] = JLD2OutputWriter(model, model.biogeochemistry.sediment_model.fields, +simulation.output_writers[:sediment] = JLD2OutputWriter(model, model.biogeochemistry.sediment.fields, indices = (:, :, 1), filename = "column_np_sediment.jld2", schedule = TimeInterval(1day), diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 9d9a051fc..a7af7c5da 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -384,7 +384,7 @@ const DOM = Union{Val{:DOM}, Val{:DON}} end end -function update_boxmodel_state!(model::BoxModel{<:LOBSTER, <:Any, <:Any, <:Any, <:Any, <:Any}) +function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:LOBSTER}, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) end diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 5b05b3410..4ac9dd4ff 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -282,7 +282,7 @@ end end end -function update_boxmodel_state!(model::BoxModel{<:NPZD, <:Any, <:Any, <:Any, <:Any, <:Any}) +function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:NPZD}, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) end From eb90278ffacd1f02295250e2cd5db96bf6db2962 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 18:15:16 +0200 Subject: [PATCH 41/68] finished 'modifier' framework --- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 6 ++-- src/Models/AdvectedPopulations/NPZD.jl | 6 ++-- src/OceanBioME.jl | 28 ++++++++++--------- 3 files changed, 23 insertions(+), 17 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index a7af7c5da..1e0ea9881 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -302,7 +302,8 @@ function LOBSTER(; grid, sinking_speeds = (sPOM = 3.47e-5, bPOM = 200/day), open_bottom::Bool = true, - particles::P = nothing) where {FT, LA, S, P} + particles::P = nothing, + modifiers::M = nothing) where {FT, LA, S, P, M} if !isnothing(sediment_model) && !open_bottom @warn "You have specified a sediment model but not `open_bottom` which will not work as the tracer will settle in the bottom cell" @@ -349,7 +350,8 @@ function LOBSTER(; grid, return Biogeochemistry(underlying_biogeochemistry; light_attenuation = light_attenuation_model, sediment = sediment_model, - particles) + particles, + modifiers) end # wrote this functionally and it took 2.5x longer so even though this is long going to use this way instead diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 4ac9dd4ff..acf8cb0ee 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -176,7 +176,8 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, sinking_speeds = (P = 0.2551/day, D = 2.7489/day), open_bottom::Bool = true, - particles::P = nothing) where {FT, LA, S, P} + particles::P = nothing, + modifiers::M = nothing) where {FT, LA, S, P, M} sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom) @@ -197,7 +198,8 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, return Biogeochemistry(underlying_biogeochemistry; light_attenuation = light_attenuation_model, sediment = sediment_model, - particles) + particles, + modifiers) end const NPZD = NutrientPhytoplanktonZooplanktonDetritus diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 4abbe0feb..b36436a8f 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -5,7 +5,7 @@ between ocean biogeochemistry, carbonate chemistry, and physics. module OceanBioME # Biogeochemistry models -export LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, PISCES, NPZD, redfield +export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, PISCES, NPZD, redfield # Macroalgae models export SLatissima @@ -44,23 +44,23 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, import Adapt: adapt_structure import Base: show, summary -struct Biogeochemistry{B, L, S, P, I} <: AbstractContinuousFormBiogeochemistry +struct Biogeochemistry{B, L, S, P, M} <: AbstractContinuousFormBiogeochemistry underlying_biogeochemistry :: B light_attenuation :: L sediment :: S particles :: P - inputs :: I + modifiers :: M Biogeochemistry(underlying_biogeochemistry::B, light_attenuation::L, sediment::S, particles::P, - inputs::I) where {B, L, S, P, I} = - new{B, L, S, P, I}(underlying_biogeochemistry, + modifiers::M) where {B, L, S, P, M} = + new{B, L, S, P, M}(underlying_biogeochemistry, light_attenuation, sediment, particles, - inputs) + modifiers) end """ @@ -68,10 +68,10 @@ end light_attenuation = nothing, sediment = nothing, particles = nothing, - inputs = nothing) + modifiers = nothing) Construct a biogeochemical model based on `underlying_biogeochemistry` which may have -a `light_attenuation` model, `sediment`, `particles`, and `inputs`. +a `light_attenuation` model, `sediment`, `particles`, and `modifiers`. Keyword Arguments ================= @@ -79,18 +79,18 @@ Keyword Arguments - `light_attenuation_model`: light attenuation model which integrated the attenuation of available light - `sediment_model`: slot for `AbstractSediment` - `particles`: slot for `BiogeochemicalParticles` -- `inputs`: slot for nutrient inputs such as rivers (work in progress) +- `modifiers`: slot for components which modfiy the biogeochemistry when the tendencies have been calculated """ Biogeochemistry(underlying_biogeochemistry; light_attenuation = nothing, sediment = nothing, particles = nothing, - inputs = nothing) = + modifiers = nothing) = Biogeochemistry(underlying_biogeochemistry, light_attenuation, sediment, particles, - inputs) + modifiers) required_biogeochemical_tracers(bgc::Biogeochemistry) = required_biogeochemical_tracers(bgc.underlying_biogeochemistry) @@ -105,14 +105,16 @@ adapt_structure(to, bgc::Biogeochemistry) = Biogeochemistry(adapt(to, bgc.underl adapt(to, bgc.light_attenuation), adapt(to, bgc.sediment), adapt(to, bgc.particles), - adapt(to, bgc.inputs)) + adapt(to, bgc.modifiers)) function update_tendencies!(bgc::Biogeochemistry, model) update_tendencies!(bgc, bgc.sediment, model) update_tendencies!(bgc, bgc.particles, model) - update_tendencies!(bgc, bgc.inputs, model) + update_tendencies!(bgc, bgc.modifiers, model) end +update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modifier, model) for modifier in modifiers] + @inline (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) function update_biogeochemical_state!(bgc::Biogeochemistry, model) From 7e0b694d1d3fa802388422a58c9ee23f6159fce4 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 18:32:38 +0200 Subject: [PATCH 42/68] updated zero negative tracers to be more reliable --- src/OceanBioME.jl | 10 +-- src/Utils/negative_tracers.jl | 132 +++++++++------------------------- 2 files changed, 40 insertions(+), 102 deletions(-) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index b36436a8f..acf94c150 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -5,7 +5,7 @@ between ocean biogeochemistry, carbonate chemistry, and physics. module OceanBioME # Biogeochemistry models -export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, PISCES, NPZD, redfield +export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, NPZD, redfield # Macroalgae models export SLatissima @@ -26,7 +26,7 @@ export Boundaries, Sediments, GasExchange, FlatSediment export column_advection_timescale, column_diffusion_timescale, sinking_advection_timescale, Budget # Positivity preservaiton utilities -export zero_negative_tracers!, error_on_neg!, warn_on_neg!, ScaleNegativeTracers, remove_NaN_tendencies! +export ScaleNegativeTracers, ZeroNegativeTracers # Oceananigans extensions export ColumnField, isacolumn @@ -79,7 +79,7 @@ Keyword Arguments - `light_attenuation_model`: light attenuation model which integrated the attenuation of available light - `sediment_model`: slot for `AbstractSediment` - `particles`: slot for `BiogeochemicalParticles` -- `modifiers`: slot for components which modfiy the biogeochemistry when the tendencies have been calculated +- `modifiers`: slot for components which modfiy the biogeochemistry when the tendencies have been calculated or when the state is updated """ Biogeochemistry(underlying_biogeochemistry; light_attenuation = nothing, @@ -113,15 +113,17 @@ function update_tendencies!(bgc::Biogeochemistry, model) update_tendencies!(bgc, bgc.modifiers, model) end +update_tendencies!(bgc, modifier, model) = nothing update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modifier, model) for modifier in modifiers] @inline (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) function update_biogeochemical_state!(bgc::Biogeochemistry, model) + update_biogeochemical_state!(model, bgc.modifiers) update_biogeochemical_state!(model, bgc.light_attenuation) end -update_tendencies!(bgc, ::Nothing, model) = nothing +update_biogeochemical_state!(model, modifiers::Tuple) = [update_biogeochemical_state!(model, modifier) for modifier in modifiers] abstract type UnderlyingBiogeochemicalModel end diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 79492f75f..7ff506767 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -5,62 +5,24 @@ using Oceananigans.Utils: work_layout using Oceananigans.Architectures: device import Adapt: adapt_structure, adapt +import Oceananigans.Biogeochemistry: update_tendencies!, update_biogeochemical_state! """ - zero_negative_tracers!(sim; params = (exclude=(), )) - -Sets any tracers in `sim.model` which are negative to zero. Use like: -```julia -simulation.callbacks[:neg] = Callback(zero_negative_tracers!) -``` - -Tracers to exclude can be set in the parameters. + ZeroNegativeTracers(; exclude = ()) +Construct a modifier that zeros any negative tracers excluding those listed in `exclude`. !!! danger "Tracer conservation" This method is _not_ recommended as a way to preserve positivity of tracers since it does not conserve the total tracer. """ -function zero_negative_tracers!(model; params = (exclude=(), )) - @unroll for (tracer_name, tracer) in pairs(model.tracers) - if !(tracer_name in params.exclude) - parent(tracer) .= max.(0.0, parent(tracer)) - end - end -end -@inline zero_negative_tracers!(sim::Simulation; params = (exclude=(), )) = zero_negative_tracers!(sim.model; params) - -""" - error_on_neg(sim; params = (exclude=(), )) - -Throws an error if any tracers in `sim.model` are negative. Use like: -```julia -simulation.callbacks[:neg] = Callback(error_on_neg!) -``` - -Tracers to exclude can be set in the parameters. -""" -function error_on_neg!(sim; params = (exclude=(), )) - @unroll for (tracer_name, tracer) in pairs(sim.model.tracers) - if !(tracer_name in params.exclude) - if any(tracer .< 0.0) error("$tracer_name < 0") end - end - end +@kwdef struct ZeroNegativeTracers{E} + exclude :: E = () end -""" - warn_on_neg(sim; params = (exclude=(), )) - -Raises a warning if any tracers in `sim.model` are negative. Use like: -```julia -simulation.callbacks[:neg] = Callback(warn_on_neg!) -``` - -Tracers to exclude can be set in the parameters. -""" -function warn_on_neg!(sim; params = (exclude=(), )) - @unroll for (tracer_name, tracer) in pairs(sim.model.tracers) - if !(tracer_name in params.exclude) - if any(tracer .< 0.0) @warn "$tracer_name < 0" end +function update_biogeochemical_state!(model, zero::ZeroNegativeTracers) + @unroll for (tracer_name, tracer) in pairs(model.tracers) + if !(tracer_name in zero.exclude) + parent(tracer) .= max.(0.0, parent(tracer)) end end end @@ -69,32 +31,6 @@ end ##### Infastructure to rescale negative values ##### -@kernel function scale_for_negs!(fields, tracers, scalefactors) - i, j, k = @index(Global, NTuple) - t, p = 0.0, 0.0 - @unroll for (idx, tracer) in enumerate(tracers) - field = @inbounds fields[tracer][i, j, k] - scalefactor = @inbounds scalefactors[idx] - - t += field * scalefactor - if field > 0 - p += field * scalefactor - end - end - t < 0 && error("Cell total < 0, can not scale negative tracers.") - @unroll for tracer in tracers - field = @inbounds fields[tracer][i, j, k] - - if field > 0 - field *= t / p - else - field = 0 - end - - @inbounds fields[tracer][i, j, k] = field - end -end - struct ScaleNegativeTracers{FA, SA, W} tracers :: FA scalefactors :: SA @@ -122,43 +58,43 @@ github](https://github.com/OceanBioME/OceanBioME.jl/discussions/48). Future plans include implement a positivity-preserving timestepping scheme as the ideal alternative. """ -function ScaleNegativeTracers(; model, tracers, scalefactors = NamedTuple{tracers}(ones(length(tracers))), warn = false) +function ScaleNegativeTracers(tracers; scalefactors = NamedTuple{tracers}(ones(length(tracers))), warn = false) if length(scalefactors) != length(tracers) error("Incorrect number of scale factors provided") end - tracers = ntuple(n -> indexin([tracers[n]], [keys(model.tracers)...])[1], length(tracers)) - scalefactors = values(scalefactors) - return ScaleNegativeTracers(tracers, scalefactors, warn) end -@inline function (scale::ScaleNegativeTracers)(model) +function update_biogeochemical_state!(model, scale::ScaleNegativeTracers) workgroup, worksize = work_layout(model.grid, :xyz) dev = device(model.grid.architecture) scale_for_negs_kernel! = scale_for_negs!(dev, workgroup, worksize) - scale_for_negs_kernel!(values(model.tracers), scale.tracers, scale.scalefactors) + scale_for_negs_kernel!(model.tracers, scale.tracers, scale.scalefactors) end -@inline (scale::ScaleNegativeTracers)(sim::Simulation) = scale(sim.model) -@kernel function _remove_NaN_tendencies!(fields) +@kernel function scale_for_negs!(fields, tracers, scalefactors) i, j, k = @index(Global, NTuple) - for field in fields - if @inbounds isnan(field[i, j, k]) - @inbounds field[i, j, k] = 0.0 - end - end -end - -""" - remove_NaN_tendencies!(model) + t, p = 0.0, 0.0 + @unroll for (idx, tracer) in enumerate(tracers) + field = @inbounds fields[tracer][i, j, k] + scalefactor = @inbounds scalefactors[tracer] -Zeros any `NaN` value tendencies as a final protection against negative tracer run away. -""" -@inline function remove_NaN_tendencies!(model) - workgroup, worksize = work_layout(model.grid, :xyz) - remove_NaN_tendencies_kernel! = _remove_NaN_tendencies!(device(model.grid.architecture), workgroup, worksize) - remove_NaN_tendencies_kernel!(values(model.timestepper.Gⁿ)) + t += field * scalefactor + if field > 0 + p += field * scalefactor + end + end + t < 0 && error("Cell total < 0, can not scale negative tracers.") + @unroll for tracer in tracers + field = @inbounds fields[tracer][i, j, k] + + if field > 0 + field *= t / p + else + field = 0 + end - return nothing -end + @inbounds fields[tracer][i, j, k] = field + end +end \ No newline at end of file From 1f6a4c11166b4e3cdaef1b6c5ebeac349938bcd0 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 18:48:22 +0200 Subject: [PATCH 43/68] added method to automatically setup negative tracer scalking --- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 14 ++++++++++-- src/Models/AdvectedPopulations/NPZD.jl | 13 +++++++++-- src/OceanBioME.jl | 7 ++++++ src/Utils/negative_tracers.jl | 22 ++++++++++++++++--- 4 files changed, 49 insertions(+), 7 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 1e0ea9881..d23ef2280 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -47,11 +47,11 @@ using Oceananigans.Units using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation -using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, UnderlyingBiogeochemicalModel +using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, UnderlyingBiogeochemicalModel, ScaleNegativeTracers using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux -import OceanBioME: redfield +import OceanBioME: redfield, conserved_tracers import OceanBioME.BoxModels: update_boxmodel_state! import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, @@ -302,6 +302,8 @@ function LOBSTER(; grid, sinking_speeds = (sPOM = 3.47e-5, bPOM = 200/day), open_bottom::Bool = true, + scale_negatives = false, + particles::P = nothing, modifiers::M = nothing) where {FT, LA, S, P, M} @@ -347,6 +349,11 @@ function LOBSTER(; grid, sinking_velocities) + if scale_negatives + scaler = ScaleNegativeTracers(underlying_biogeochemistry) + modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler) + end + return Biogeochemistry(underlying_biogeochemistry; light_attenuation = light_attenuation_model, sediment = sediment_model, @@ -473,4 +480,7 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true sinking_flux(i, j, k, grid, advection, Val(:bPOC), bgc, tracers) @inline remineralisation_receiver(::LOBSTER) = :NH₄ + +@inline conserved_tracers(::LOBSTER) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) +@inline conserved_tracers(::lobster_variable_redfield) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) end # module diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index acf8cb0ee..ab6c4e651 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -16,7 +16,7 @@ module NPZDModel export NutrientPhytoplanktonZooplanktonDetritus, NPZD -using OceanBioME: Biogeochemistry, UnderlyingBiogeochemicalModel +using OceanBioME: Biogeochemistry, UnderlyingBiogeochemicalModel, ScaleNegativeTracers using Oceananigans.Units using Oceananigans.Fields: ZeroField @@ -26,7 +26,7 @@ using OceanBioME: setup_velocity_fields, show_sinking_velocities using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux -import OceanBioME: redfield +import OceanBioME: redfield, conserved_tracers import OceanBioME.BoxModels: update_boxmodel_state! import Base: show, summary @@ -175,6 +175,8 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, sinking_speeds = (P = 0.2551/day, D = 2.7489/day), open_bottom::Bool = true, + + scale_negatives = false, particles::P = nothing, modifiers::M = nothing) where {FT, LA, S, P, M} @@ -195,6 +197,11 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, remineralization_rate, sinking_velocities) + if scale_negatives + scaler = ScaleNegativeTracers(underlying_biogeochemistry) + modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler) + end + return Biogeochemistry(underlying_biogeochemistry; light_attenuation = light_attenuation_model, sediment = sediment_model, @@ -322,4 +329,6 @@ adapt_structure(to, npzd::NPZD) = @inline carbon_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * redfield(Val(:P), bgc) @inline remineralisation_receiver(::NPZD) = :N + +@inline conserved_tracers(::NPZD) = (:N, :P, :Z, :D) end # module diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index acf94c150..b24346309 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -150,6 +150,13 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros @inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) @inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) +""" + conserved_tracers(model::UnderlyingBiogeochemicalModel) + +Returns the names of tracers which together are conserved in `model` +""" +conserved_tracers() = nothing + summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") show(io::IO, model::Biogeochemistry) = print(io, show(model.underlying_biogeochemistry), " \n", diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 7ff506767..77dadb37c 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -47,16 +47,18 @@ adapt_structure(to, snt::ScaleNegativeTracers) = ScaleNegativeTracers(adapt(to, """ ScaleNegativeTracers(; tracers, scalefactors = ones(length(tracers)), warn = false) -Returns a callback that scales `tracers` so that none are negative. Use like: +Constructs a modifier to scale `tracers` so that none are negative. Use like: ```julia -negativity_protection! = ScaleNegativeTracers(; model, tracers = (:P, :Z, :N)) -simulation.callbacks[:neg] = Callback(negativity_protection!; callsite = UpdateStateCallsite()) +modifier = ScaleNegativeTracers((:P, :Z, :N)) +biogeochemistry = Biogeochemistry(...; modifier) ``` This method is better, though still imperfect, method to prevent numerical errors that lead to negative tracer values compared to [`zero_negative_tracers!`](@ref). Please see [discussion in github](https://github.com/OceanBioME/OceanBioME.jl/discussions/48). Future plans include implement a positivity-preserving timestepping scheme as the ideal alternative. + +If `warn` is true then scaling will raise a warning. """ function ScaleNegativeTracers(tracers; scalefactors = NamedTuple{tracers}(ones(length(tracers))), warn = false) if length(scalefactors) != length(tracers) @@ -66,6 +68,20 @@ function ScaleNegativeTracers(tracers; scalefactors = NamedTuple{tracers}(ones(l return ScaleNegativeTracers(tracers, scalefactors, warn) end +""" + ScaleNegativeTracers(model::UnderlyingBiogeochemicalModel; warn = false) + +Constructs a modifier to scale the conserved tracers in `model`. + +If `warn` is true then scaling will raise a warning. +""" +function ScaleNegativeTracers(model::UnderlyingBiogeochemicalModel; warn = false) + tracers = conserved_tracers(model) + scalefactors = NamedTuple{tracers}(ones(length(tracers))) + + return ScaleNegativeTracers(tracers, scalefactors, warn) +end + function update_biogeochemical_state!(model, scale::ScaleNegativeTracers) workgroup, worksize = work_layout(model.grid, :xyz) dev = device(model.grid.architecture) From ca4cf2e87aa943ded6f7f2547df42f94a72a3778 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 18:55:35 +0200 Subject: [PATCH 44/68] updated negative protection tests --- test/test_utils.jl | 35 +++++++++++------------------------ 1 file changed, 11 insertions(+), 24 deletions(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index 1245c8f69..86586a791 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -20,48 +20,35 @@ end function test_negative_scaling(arch) grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) - model = NonhydrostaticModel(; grid, tracers = (:A, :B)) + model = NonhydrostaticModel(; grid, biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid, scale_negatives = true)) - set!(model, A = 1, B = -0.5) + set!(model, N = 2, P = -1) - simulation = Simulation(model, Δt = 1) + simulation = Simulation(model, Δt = 1e-10, stop_iteration = 1) - negativity_protection! = ScaleNegativeTracers(; model, tracers = (:A, :B)) + run!(simulation) - negativity_protection!(simulation) - - return (model.tracers.A[1, 1, 1] == 0.5) && (model.tracers.B[1, 1, 1] == 0.0) -end - -function test_negative_scaling(arch) - grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) - - model = NonhydrostaticModel(; grid, tracers = :A) - - model.timestepper.Gⁿ.A .= NaN - - remove_NaN_tendencies!(model) - - return model.timestepper.Gⁿ.A[1, 1, 1] == 0.0 + return (model.tracers.N[1, 1, 1] ≈ 1) && (model.tracers.P[1, 1, 1] ≈ 0.0) end function test_negative_zeroing(arch) grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) - model = NonhydrostaticModel(; grid, tracers = (:A, :B)) + model = NonhydrostaticModel(; grid, biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid, modifiers = ZeroNegativeTracers(; exclude = (:Z, )))) - set!(model, A = -1, B = -0.5) + set!(model, N = 2, P = -1, Z = -1) - zero_negative_tracers!(model, params = (exclude = (:B, ), )) + simulation = Simulation(model, Δt = 1e-10, stop_iteration = 1) - return (model.tracers.A[1, 1, 1] == 0.0) && (model.tracers.B[1, 1, 1] == -0.5) + run!(simulation) + + return (model.tracers.N[1, 1, 1] ≈ 2) && (model.tracers.P[1, 1, 1] ≈ 0.0) && (model.tracers.Z[1, 1, 1] ≈ -1) end @testset "Test Utils" begin for arch in (CPU(), ) @test test_column_diffusion_timescale(arch) @test test_negative_scaling(arch) - @test test_negative_scaling(arch) @test test_negative_zeroing(arch) end end From c0ac3c16632614ecb4530bde86c3fb3d15d96028 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 18:55:52 +0200 Subject: [PATCH 45/68] updated examples --- examples/column.jl | 20 ++++++++++++-------- examples/data_forced.jl | 11 ++++++++--- examples/eady.jl | 6 ------ examples/kelp.jl | 3 +-- 4 files changed, 21 insertions(+), 19 deletions(-) diff --git a/examples/column.jl b/examples/column.jl index 201490814..412b459ca 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -40,14 +40,22 @@ nothing #hide # Define the grid. grid = RectilinearGrid(size=(1, 1, 50), extent=(20meters, 20meters, 200meters)) -# Specify the boundary conditions for DIC and O₂ based on the air-sea CO₂ and O₂ flux + +# ## Model +# First we will define the biogeochemical model including carbonate chemistry +# and scaling of negative tracers(see discussion in the [positivity preservation](@ref pos-preservation)) +# and then setup the Oceananigans model with the boundary condition for the DIC based on the air-sea CO₂ flux + +biogeochemistry = LOBSTER(; grid, + surface_phytosynthetically_active_radiation = PAR⁰, + carbonates = true, + scale_negatives = true) + CO₂_flux = GasExchange(; gas = :CO₂, temperature = temp, salinity = (args...) -> 35) model = NonhydrostaticModel(; grid, closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), - biogeochemistry = LOBSTER(; grid, - surface_phytosynthetically_active_radiation = PAR⁰, - carbonates = true), + biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), )) set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2409.0) @@ -56,7 +64,6 @@ set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2 # Next we setup a simulation and add some callbacks that: # - Show the progress of the simulation # - Store the model and particles output -# - Prevent the tracers from going negative from numerical error (see discussion in the [positivity preservation](@ref pos-preservation) implementation section) simulation = Simulation(model, Δt = 3minutes, stop_time = 100days) @@ -73,9 +80,6 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.trace filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing = true) - -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) nothing #hide # ## Run! diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 19e368283..2fa544b9b 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -65,12 +65,17 @@ grid = RectilinearGrid(size = (1, 1, Nz), x = (0, 20meters), y = (0, 20meters), # ## Biogeochemical and Oceananigans model # Here we instantiate the LOBSTER model with carbonate chemistry and a surface flux of DIC (CO₂) + +biogeochemistry = LOBSTER(; grid, + surface_phytosynthetically_active_radiation = surface_PAR, + carbonates = true, + scale_negatives = true) + CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = s_function) + model = NonhydrostaticModel(; grid, closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), - biogeochemistry = LOBSTER(; grid, - surface_phytosynthetically_active_radiation = surface_PAR, - carbonates = true), + biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),)) set!(model, P = 0.03, Z = 0.03, NO₃ = 11.0, NH₄ = 0.05, DIC = 2200.0, Alk = 2400.0) diff --git a/examples/eady.jl b/examples/eady.jl index 15a729dfe..260e907b8 100644 --- a/examples/eady.jl +++ b/examples/eady.jl @@ -119,12 +119,6 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.tracers filename = "eady_turbulence_bgc", overwrite_existing = true) -# Prevent the tracer values going negative - this is especially important in this model while no positivity -# preserving diffusion is implemented. -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) -simulation.callbacks[:nan_tendencies] = Callback(remove_NaN_tendencies!; callsite = TendencyCallsite()) -simulation.callbacks[:abort_zeros] = Callback(zero_negative_tracers!; callsite = UpdateStateCallsite()) nothing #hide # Run the simulation diff --git a/examples/kelp.jl b/examples/kelp.jl index 3a10a155c..76462e478 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -68,6 +68,7 @@ biogeochemistry = LOBSTER(; grid, surface_phytosynthetically_active_radiation = PAR⁰, carbonates = true, variable_redfield = true, + scale_negatives = true, particles) model = NonhydrostaticModel(; grid, @@ -103,8 +104,6 @@ simulation.output_writers[:particles] = JLD2OutputWriter(model, (; particles), schedule = TimeInterval(1day), overwrite_existing = true) -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) nothing #hide # ## Run! From 1a2b8232861644d14b87158cf4e13b8e9b766a40 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 30 Aug 2023 19:34:25 +0200 Subject: [PATCH 46/68] attempted to fix docs --- docs/src/model_components/utils.md | 17 +++++++---------- docs/src/model_implementation.md | 8 ++++---- examples/data_forced.jl | 4 +--- 3 files changed, 12 insertions(+), 17 deletions(-) diff --git a/docs/src/model_components/utils.md b/docs/src/model_components/utils.md index 8f542a928..01dba5a39 100644 --- a/docs/src/model_components/utils.md +++ b/docs/src/model_components/utils.md @@ -22,20 +22,17 @@ simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) ## Negative tracer detection As a temporary measure we have implemented a callback to either detect negative tracers and either scale a conserved group, force them back to zero, or throw an error. Please see the numerical implementations' page for details. This can be set up by: ```julia -negativity_protection! = ScaleNegativeTracers(tracers = (:P, :Z, :N)) -simulation.callbacks[:neg] = Callback(negativity_protection!; callsite = UpdateStateCallsite()) +negativity_protection = ScaleNegativeTracers((:P, :Z, :N)) +biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` You may also pass a scale factor for each component (e.g. in case they have different redfield ratios): ```julia -negativity_protection! = ScaleNegativeTracers(tracers = (:P, :Z, :D), scalefactors = (P = 1, Z = 1, D = 2)) -simulation.callbacks[:neg] = Callback(negativity_protection!; callsite = UpdateStateCallsite()) +negativity_protection = ScaleNegativeTracers((:P, :Z, :N); scalefactors = (1, 1, 2)) +biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` Here you should carefully consider which tracers form a conserved group (if at all). Alternatively, force to zero by: ```julia -simulation.callbacks[:neg] = Callback(OceanBioME.no_negative_tracers!, callsite = UpdateStateCallsite()) +negativity_protection = ZeroNegativeTracers() +biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` -or throw an error: -```julia -simulation.callbacks[:neg] = Callback(OceanBioME.error_on_neg!, callsite = UpdateStateCallsite()) -``` -The latter two both optionally take a named tuple of parameters which may include `exclude` which can be a tuple of tracer names (Symbols) which are allowed to be negative. +The latter optionally takes a named tuple of parameters which may include `exclude` which can be a tuple of tracer names (Symbols) which are allowed to be negative. diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 03b4968e8..4c10b85b2 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -253,9 +253,12 @@ w_sink(x, y, z) = 2 / day * tanh(z / 5) set!(sinking_velocity, w_sink) +negative_tracer_scaling = ScaleNegativeTracers((:N, :P)) + biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; sinking_velocity); light_attenuation, - sediment) + sediment, + modifiers = negative_tracer_scaling) # put the model together @@ -281,9 +284,6 @@ simulation.output_writers[:sediment] = JLD2OutputWriter(model, model.biogeochemi schedule = TimeInterval(1day), overwrite_existing = true) -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:N, :P)) -simulation.callbacks[:nan_tendencies] = Callback(remove_NaN_tendencies!; callsite = TendencyCallsite()) - run!(simulation) ``` diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 2fa544b9b..80363127d 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -106,13 +106,11 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, # TODO: make tendency callback to force no NaNs in tendencies -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) - wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 1.5, min_change = 0.75, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) + simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) nothing #hide From 6e37c82abbd07c23b291029a52091dca89752e10 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 30 Aug 2023 22:48:37 +0200 Subject: [PATCH 47/68] Update negative_tracers.jl --- src/Utils/negative_tracers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 77dadb37c..cbabb6141 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -53,7 +53,7 @@ modifier = ScaleNegativeTracers((:P, :Z, :N)) biogeochemistry = Biogeochemistry(...; modifier) ``` This method is better, though still imperfect, method to prevent numerical errors that lead to -negative tracer values compared to [`zero_negative_tracers!`](@ref). Please see [discussion in +negative tracer values compared to [`ZeroNegativeTracers`](@ref). Please see [discussion in github](https://github.com/OceanBioME/OceanBioME.jl/discussions/48). Future plans include implement a positivity-preserving timestepping scheme as the ideal alternative. From 3a71aaf5dd8b8106ec917c0937774c1d9a66e31b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 31 Aug 2023 09:57:13 +0200 Subject: [PATCH 48/68] hide output after loading exampel resoltes --- examples/column.jl | 2 +- examples/data_forced.jl | 4 ++-- examples/eady.jl | 2 +- examples/kelp.jl | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/column.jl b/examples/column.jl index 412b459ca..b005d491f 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -98,7 +98,7 @@ bPOM = FieldTimeSeries("$filename.jld2", "bPOM") Alk = FieldTimeSeries("$filename.jld2", "Alk") x, y, z = nodes(P) -times = P.times +times = P.times; # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 80363127d..bd7cf2145 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -110,7 +110,7 @@ wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 1.5, min_change = 0.75, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) - + simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) nothing #hide @@ -130,7 +130,7 @@ bPOM = FieldTimeSeries("$filename.jld2", "bPOM") Alk = FieldTimeSeries("$filename.jld2", "Alk") x, y, z = nodes(P) -times = P.times +times = P.times; # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth diff --git a/examples/eady.jl b/examples/eady.jl index 260e907b8..e528f0968 100644 --- a/examples/eady.jl +++ b/examples/eady.jl @@ -134,7 +134,7 @@ DIC = FieldTimeSeries("eady_turbulence_bgc.jld2", "DIC") times = ζ.times xζ, yζ, zζ = nodes(ζ) -xc, yc, zc = nodes(P) +xc, yc, zc = nodes(P); # and plot. diff --git a/examples/kelp.jl b/examples/kelp.jl index 76462e478..845a979c2 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -121,7 +121,7 @@ bPOC = FieldTimeSeries("$filename.jld2", "bPOC") Alk = FieldTimeSeries("$filename.jld2", "Alk") x, y, z = nodes(P) -times = P.times +times = P.times; # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth From 0deb0a149eba4fccfc0685d96db65588e67d6a93 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 31 Aug 2023 16:37:13 +0200 Subject: [PATCH 49/68] simplified LOBSTER tests --- src/OceanBioME.jl | 3 ++- test/test_LOBSTER.jl | 51 ++++++++++++++------------------------------ 2 files changed, 18 insertions(+), 36 deletions(-) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index b24346309..b0bfde7d5 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -149,13 +149,14 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros @inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) @inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) +@inline redfield(val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc) """ conserved_tracers(model::UnderlyingBiogeochemicalModel) Returns the names of tracers which together are conserved in `model` """ -conserved_tracers() = nothing +conserved_tracers(model::Biogeochemistry) = conserved_tracers(model.underlying_biogeochemistry) summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") show(io::IO, model::Biogeochemistry) = diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index ed667b850..6e3e4efdb 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -1,18 +1,11 @@ using Test -using OceanBioME: LOBSTER +using OceanBioME: LOBSTER, conserved_tracers, redfield using Oceananigans, CUDA -ΣN(model, variable_redfield) = variable_redfield ? ( - sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.sPON) + sum(model.tracers.bPON) + sum(model.tracers.DON)) : ( - sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.sPOM) + sum(model.tracers.bPOM) + sum(model.tracers.DOM)) +ΣN(model, biogeochemistry) = sum([sum(model.tracers[tracer_name]) for tracer_name in conserved_tracers(biogeochemistry)]) -function ΣC(model, carbonates, variable_redfield) - # *will only be conserved if carbonates on* - if variable_redfield - OC = sum(model.tracers.sPOC .+ model.tracers.bPOC .+ model.tracers.DOC) - else - OC = sum(model.tracers.sPOM .+ model.tracers.bPOM .+ model.tracers.DOM) * model.biogeochemistry.underlying_biogeochemistry.organic_redfield - end +function ΣC(model, carbonates, biogeochemistry) + OC = sum([sum(model.tracers[tracer_name] * redfield(Val(tracer_name), biogeochemistry, model.tracers)) for tracer_name in conserved_tracers(biogeochemistry)]) if carbonates IC = sum(model.tracers.DIC) @@ -20,23 +13,13 @@ function ΣC(model, carbonates, variable_redfield) IC = 0.0 end - LC = sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield - - return OC + IC + LC + return OC + IC end -ΣGⁿ(model, variable_redfield) = variable_redfield ? ( - sum(model.timestepper.Gⁿ.NO₃) + sum(model.timestepper.Gⁿ.NH₄) + sum(model.timestepper.Gⁿ.P) + sum(model.timestepper.Gⁿ.Z) + sum(model.timestepper.Gⁿ.sPON) + sum(model.timestepper.Gⁿ.bPON) + sum(model.timestepper.Gⁿ.DON)) : ( - sum(model.timestepper.Gⁿ.NO₃) + sum(model.timestepper.Gⁿ.NH₄) + sum(model.timestepper.Gⁿ.P) + sum(model.timestepper.Gⁿ.Z) + sum(model.timestepper.Gⁿ.sPOM) + sum(model.timestepper.Gⁿ.bPOM) + sum(model.timestepper.Gⁿ.DOM)) +ΣGⁿ(model, biogeochemistry) = sum([sum(model.timestepper.Gⁿ[tracer_name]) for tracer_name in conserved_tracers(biogeochemistry)]) - -function ΣGᶜ(model, carbonates, variable_redfield) - # *will only be conserved if carbonates on* - if variable_redfield - OC = sum(model.timestepper.Gⁿ.sPOC .+ model.timestepper.Gⁿ.bPOC .+ model.timestepper.Gⁿ.DOC) - else - OC = sum(model.timestepper.Gⁿ.sPOM .+ model.timestepper.Gⁿ.bPOM .+ model.timestepper.Gⁿ.DOM) * model.biogeochemistry.underlying_biogeochemistry.organic_redfield - end +function ΣGᶜ(model, carbonates, biogeochemistry) + OC = sum([sum(model.timestepper.Gⁿ[tracer_name] * redfield(Val(tracer_name), biogeochemistry, model.tracers)) for tracer_name in conserved_tracers(biogeochemistry)]) if carbonates IC = sum(model.timestepper.Gⁿ.DIC) @@ -44,9 +27,7 @@ function ΣGᶜ(model, carbonates, variable_redfield) IC = 0.0 end - LC = sum(model.timestepper.Gⁿ.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.timestepper.Gⁿ.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield - - return OC + IC + LC + return OC + IC end function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open_bottom, n_timesteps) @@ -59,7 +40,7 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open end # correct tracers and auxiliary fields have been setup, and order has not changed - required_tracers = variable_redfield ? (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) : (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) + required_tracers = conserved_tracers(biogeochemistry) if carbonates required_tracers = (required_tracers..., :DIC, :Alk) end @@ -105,25 +86,25 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open end end - ΣN₀ = CUDA.@allowscalar ΣN(model, variable_redfield) + ΣN₀ = CUDA.@allowscalar ΣN(model, model.biogeochemistry) - ΣC₀ = CUDA.@allowscalar ΣC(model, carbonates, variable_redfield) + ΣC₀ = CUDA.@allowscalar ΣC(model, carbonates, model.biogeochemistry) for _ in 1:n_timesteps time_step!(model, 1.0) end - ΣN₁ = CUDA.@allowscalar ΣN(model, variable_redfield) + ΣN₁ = CUDA.@allowscalar ΣN(model, model.biogeochemistry) CUDA.@allowscalar if !(sinking && open_bottom) #when we have open bottom sinking we won't conserve anything @test ΣN₀ ≈ ΣN₁ - @test ΣGⁿ(model, variable_redfield) ≈ 0.0 atol = 1e-15 # rtol=sqrt(eps) so is usually much larger than even this + @test ΣGⁿ(model, model.biogeochemistry) ≈ 0.0 atol = 1e-15 # rtol=sqrt(eps) so is usually much larger than even this if (carbonates && variable_redfield) - ΣC₁ = ΣC(model, carbonates, variable_redfield) + ΣC₁ = ΣC(model, carbonates, model.biogeochemistry) @test ΣC₀ ≈ ΣC₁# atol = 0.0001 # when we convert to and from - @test ΣGᶜ(model, carbonates, variable_redfield) ≈ 0.0 atol = 1e-15 + @test ΣGᶜ(model, carbonates, model.biogeochemistry) ≈ 0.0 atol = 1e-15 end end return model From 8cbfe2e60b61e20dd0c18f9d5521e5cd9b33014b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 31 Aug 2023 18:05:10 +0200 Subject: [PATCH 50/68] oops --- test/test_LOBSTER.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index 6e3e4efdb..c872a1fba 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -40,7 +40,7 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open end # correct tracers and auxiliary fields have been setup, and order has not changed - required_tracers = conserved_tracers(biogeochemistry) + required_tracers = conserved_tracers(model.biogeochemistry) if carbonates required_tracers = (required_tracers..., :DIC, :Alk) end From f58b7f9cd2faba50024442c961ef6aaaea740438 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 31 Aug 2023 20:18:44 +0200 Subject: [PATCH 51/68] fixed LOBSTER test --- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 4 +++ src/OceanBioME.jl | 5 +++- src/Utils/sinking_velocity_fields.jl | 9 ++++-- test/test_LOBSTER.jl | 28 ++++++++++--------- 4 files changed, 29 insertions(+), 17 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index d23ef2280..f34b13551 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -465,6 +465,10 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true @inline redfield(i, j, k, ::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.bPOC[i, j, k] / tracers.bPON[i, j, k] @inline redfield(i, j, k, ::Val{:DON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.DOC[i, j, k] / tracers.DON[i, j, k] +@inline redfield(::Val{:sPON}, bgc::lobster_variable_redfield, tracers) = tracers.sPOC / tracers.sPON +@inline redfield(::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = tracers.bPOC / tracers.bPON +@inline redfield(::Val{:DON}, bgc::lobster_variable_redfield, tracers) = tracers.DOC / tracers.DON + @inline nitrogen_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = sinking_flux(i, j, k, grid, advection, Val(:sPOM), bgc, tracers) + sinking_flux(i, j, k, grid, advection, Val(:bPOM), bgc, tracers) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index b0bfde7d5..fff4d6e37 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -147,9 +147,12 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros """ @inline redfield(val_tracer_name, bgc) = NaN +# fallbacks @inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) +@inline redfield(i, j, k, val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc, tracers) @inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) -@inline redfield(val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc) +@inline redfield(val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(val_tracer_name, bgc.underlying_biogeochemistry, tracers) +@inline redfield(val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc) """ conserved_tracers(model::UnderlyingBiogeochemicalModel) diff --git a/src/Utils/sinking_velocity_fields.jl b/src/Utils/sinking_velocity_fields.jl index e4d410eea..2a2805814 100644 --- a/src/Utils/sinking_velocity_fields.jl +++ b/src/Utils/sinking_velocity_fields.jl @@ -30,15 +30,18 @@ adapt_structure(to, velocities::NamedTuple{(:u, :v, :w), Tuple{AbstractField, Ab function show_sinking_velocities(sinking_velocities::NamedTuple{T, V}) where {T, V} str = "" - if length(T) == 1 - str = " └── $(T[1]): $(maximum_sinking(sinking_velocities[1])) to $(minimum_sinking(sinking_velocities[1])) m/s" + if length(T) == 0 + return str + elseif length(T) == 1 + return " └── $(T[1]): $(maximum_sinking(sinking_velocities[1])) to $(minimum_sinking(sinking_velocities[1])) m/s" else for idx in 1:length(T) - 1 str *= " ├── $(T[idx]): $(maximum_sinking(sinking_velocities[idx])) to $(minimum_sinking(sinking_velocities[idx])) m/s \n" end str *= " └── $(T[end]): $(maximum_sinking(sinking_velocities[end])) to $(minimum_sinking(sinking_velocities[end])) m/s" + + return str end - return str end maximum_sinking(velocity) = maximum(velocity) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index c872a1fba..db25b6d49 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -18,16 +18,18 @@ end ΣGⁿ(model, biogeochemistry) = sum([sum(model.timestepper.Gⁿ[tracer_name]) for tracer_name in conserved_tracers(biogeochemistry)]) -function ΣGᶜ(model, carbonates, biogeochemistry) - OC = sum([sum(model.timestepper.Gⁿ[tracer_name] * redfield(Val(tracer_name), biogeochemistry, model.tracers)) for tracer_name in conserved_tracers(biogeochemistry)]) - if carbonates - IC = sum(model.timestepper.Gⁿ.DIC) - else - IC = 0.0 - end +function ΣGᶜ(model, biogeochemistry) + # *will only be conserved if carbonates on* + # Leaving this as is for now since as the `redfield` is based on value not tendency then it returns the wrong results if + # the redfield is rapidly changing + OC = sum(model.timestepper.Gⁿ.sPOC .+ model.timestepper.Gⁿ.bPOC .+ model.timestepper.Gⁿ.DOC) - return OC + IC + IC = sum(model.timestepper.Gⁿ.DIC) + + LC = sum(model.timestepper.Gⁿ.P * (1 + biogeochemistry.organic_carbon_calcate_ratio) .+ model.timestepper.Gⁿ.Z) * biogeochemistry.phytoplankton_redfield + + return OC + IC + LC end function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open_bottom, n_timesteps) @@ -71,9 +73,9 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open model.tracers.bPON .= rand() model.tracers.DON .= rand() - model.tracers.sPOC .= model.tracers.sPON * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield - model.tracers.bPOC .= model.tracers.bPON * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield - model.tracers.DOC .= model.tracers.DON * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield + model.tracers.sPOC .= model.tracers.sPON * redfield(Val(:sPOM), model.biogeochemistry) + model.tracers.bPOC .= model.tracers.bPON * redfield(Val(:bPOM), model.biogeochemistry) + model.tracers.DOC .= model.tracers.DON * redfield(Val(:DOM), model.biogeochemistry) else model.tracers.sPOM .= rand() model.tracers.bPOM .= rand() @@ -102,9 +104,9 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open if (carbonates && variable_redfield) ΣC₁ = ΣC(model, carbonates, model.biogeochemistry) - @test ΣC₀ ≈ ΣC₁# atol = 0.0001 # when we convert to and from + @test ΣC₀ ≈ ΣC₁ - @test ΣGᶜ(model, carbonates, model.biogeochemistry) ≈ 0.0 atol = 1e-15 + @test ΣGᶜ(model, model.biogeochemistry.underlying_biogeochemistry) ≈ 0.0 atol = 1e-15 end end return model From ba0d49f3dda8b8c711f471d4c18bfc30621c0202 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 31 Aug 2023 21:45:24 +0200 Subject: [PATCH 52/68] Updated bgc docs page --- docs/src/model_components/biogeochemical/index.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/src/model_components/biogeochemical/index.md b/docs/src/model_components/biogeochemical/index.md index b3a74f495..8ad577778 100644 --- a/docs/src/model_components/biogeochemical/index.md +++ b/docs/src/model_components/biogeochemical/index.md @@ -5,7 +5,7 @@ Biogeochemical (BGC) models can be used within the [Oceananigans biogeochemistry For details of the BGC models currently implemented please see the following pages. ## Oceananigans setup -At the simplest level all that is required to setup a BGC model is to pass it to the Oceananigans model setup: +At the simplest level all that is required to setup an existint OceanBioME BGC model is to pass it to the Oceananigans model setup: ```julia model = NonhydrostaticModel(; grid, ..., @@ -19,3 +19,7 @@ MODEL_NAME(; grid, growth_rate = 10.0) This will set up the required tracers and auxiliary fields, and you may also set boundary conditions or additional forcing through the normal Oceananigans setup. Models usually have a default [light attenuation model](@ref light) specified, these may be substituted easily by passing different models as parameters as above. + +Our models are implemented in an abstract framework `Biogeochemistry` which contains `underlying_biogeochemistry`, `light_attenuation`, `sediment`, and `modifiers`. +This is automatically setup for existing BGC models, but may also be used to couple any BGC model with light attenuation and sediments. +See the [implementation page](@ref model_implementation) for some more information on how to couple other models. \ No newline at end of file From 87c7ae4a039c1c7bc7d81891693b80f461dc84b6 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 31 Aug 2023 22:34:23 +0200 Subject: [PATCH 53/68] Update index.md --- docs/src/model_components/biogeochemical/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_components/biogeochemical/index.md b/docs/src/model_components/biogeochemical/index.md index 8ad577778..191e58994 100644 --- a/docs/src/model_components/biogeochemical/index.md +++ b/docs/src/model_components/biogeochemical/index.md @@ -5,7 +5,7 @@ Biogeochemical (BGC) models can be used within the [Oceananigans biogeochemistry For details of the BGC models currently implemented please see the following pages. ## Oceananigans setup -At the simplest level all that is required to setup an existint OceanBioME BGC model is to pass it to the Oceananigans model setup: +At the simplest level all that is required to setup an existing OceanBioME BGC model is to pass it to the Oceananigans model setup: ```julia model = NonhydrostaticModel(; grid, ..., From 577f05c8691a1149c2b74bf6ce79ed1d5f052c8c Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 31 Aug 2023 23:50:26 +0200 Subject: [PATCH 54/68] Bump minor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9b501f041..cd1ca939c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright ", "John R Taylor ", "Si Chen "] -version = "0.5.0" +version = "0.6.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From a5fc7e2cd39e9371f7fc0861ff061893af30df4f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 1 Sep 2023 09:06:17 +0200 Subject: [PATCH 55/68] fixed julia version downgrade --- Manifest.toml | 6 +++--- src/OceanBioME.jl | 3 +-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index bcc3f050c..fc62c33f4 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.1" +julia_version = "1.9.3" manifest_format = "2.0" project_hash = "b8336a615981d764881d88f0d8193dd7f4667beb" @@ -140,7 +140,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.2+0" +version = "1.0.5+0" [[deps.ConstructionBase]] deps = ["LinearAlgebra"] @@ -572,7 +572,7 @@ version = "0.15.1" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.0" +version = "1.9.2" [[deps.PkgVersion]] deps = ["Pkg"] diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index fff4d6e37..d3dddac7a 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -138,7 +138,7 @@ struct BoxModelGrid end Returns the redfield ratio of `tracer_name` from `bgc` at `i`, `j`, `k`. """ -@inline redfield(i, j, k, val_tracer_name, bgc, tracers) = NaN +@inline redfield(i, j, k, val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc, tracers) """ redfield(val_tracer_name, bgc, tracers) @@ -149,7 +149,6 @@ Returns the redfield ratio of `tracer_name` from `bgc` when it is constant acros # fallbacks @inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) -@inline redfield(i, j, k, val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc, tracers) @inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) @inline redfield(val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(val_tracer_name, bgc.underlying_biogeochemistry, tracers) @inline redfield(val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc) From 0191dc103d23f8ed96cdb98fc1e0b8baf4995ae5 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 2 Sep 2023 13:59:49 +0200 Subject: [PATCH 56/68] dropped `UnderlyingBiogeochemistry` in favor of better using the Oceananigans Asbtract types This means any Ocenanigans BGC can go in the slot, and we don't need to define some extra things --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 6 ++++-- src/Models/AdvectedPopulations/NPZD.jl | 5 +++-- src/OceanBioME.jl | 10 ++++------ src/Utils/negative_tracers.jl | 3 ++- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index f34b13551..9c0a913a7 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -47,10 +47,12 @@ using Oceananigans.Units using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation -using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, UnderlyingBiogeochemicalModel, ScaleNegativeTracers +using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, ScaleNegativeTracers using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux +using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry + import OceanBioME: redfield, conserved_tracers import OceanBioME.BoxModels: update_boxmodel_state! @@ -65,7 +67,7 @@ import Base: show, summary import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver -struct LOBSTER{FT, B, W} <: UnderlyingBiogeochemicalModel +struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry phytoplankton_preference :: FT maximum_grazing_rate :: FT grazing_half_saturation :: FT diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index ab6c4e651..ed892a72e 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -16,7 +16,8 @@ module NPZDModel export NutrientPhytoplanktonZooplanktonDetritus, NPZD -using OceanBioME: Biogeochemistry, UnderlyingBiogeochemicalModel, ScaleNegativeTracers +using OceanBioME: Biogeochemistry, ScaleNegativeTracers +using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry using Oceananigans.Units using Oceananigans.Fields: ZeroField @@ -40,7 +41,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa import Adapt: adapt_structure, adapt -struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} <: UnderlyingBiogeochemicalModel +struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} <: AbstractContinuousFormBiogeochemistry # phytoplankton initial_photosynthetic_slope :: FT # α, 1/(W/m²)/s base_maximum_growth :: FT # μ₀, 1/s diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index d3dddac7a..e7093847f 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -39,7 +39,8 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, update_biogeochemical_state!, biogeochemical_drift_velocity, biogeochemical_auxiliary_fields, - update_tendencies! + update_tendencies!, + biogeochemical_transition import Adapt: adapt_structure import Base: show, summary @@ -116,7 +117,8 @@ end update_tendencies!(bgc, modifier, model) = nothing update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modifier, model) for modifier in modifiers] -@inline (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) +@inline biogeochemical_transition(i, j, k, grid, bgc::Biogeochemistry, val_tracer_name, clock, fields) = + biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) function update_biogeochemical_state!(bgc::Biogeochemistry, model) update_biogeochemical_state!(model, bgc.modifiers) @@ -125,10 +127,6 @@ end update_biogeochemical_state!(model, modifiers::Tuple) = [update_biogeochemical_state!(model, modifier) for modifier in modifiers] -abstract type UnderlyingBiogeochemicalModel end - -@inline (::UnderlyingBiogeochemicalModel)(val_tracer_name, x, y, z, t, fields...) = zero(x) - struct BoxModelGrid end @inline maximum_sinking_velocity(bgc) = 0.0 diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index cbabb6141..11839fd62 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -3,6 +3,7 @@ using KernelAbstractions using KernelAbstractions.Extras.LoopInfo: @unroll using Oceananigans.Utils: work_layout using Oceananigans.Architectures: device +using Oceananigans.Biogeochemistry: AbstractBiogeochemistry import Adapt: adapt_structure, adapt import Oceananigans.Biogeochemistry: update_tendencies!, update_biogeochemical_state! @@ -75,7 +76,7 @@ Constructs a modifier to scale the conserved tracers in `model`. If `warn` is true then scaling will raise a warning. """ -function ScaleNegativeTracers(model::UnderlyingBiogeochemicalModel; warn = false) +function ScaleNegativeTracers(model::AbstractBiogeochemistry; warn = false) tracers = conserved_tracers(model) scalefactors = NamedTuple{tracers}(ones(length(tracers))) From a9dac4b6e0f467328fc977a5e9aff307d640df00 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 2 Sep 2023 14:07:45 +0200 Subject: [PATCH 57/68] updated model implementation page --- docs/src/model_implementation.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 4c10b85b2..59fcd5185 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -18,7 +18,8 @@ For this example we are going to implement the simple Nutrient-Phytoplankton mod The first step is to import the abstract type from OceanBioME, some units from Oceananigans (for ease of parameter definition), and [`import`](https://stackoverflow.com/questions/27086159/what-is-the-difference-between-using-and-import-in-julia-when-building-a-mod) some functions from Oceananigans in order to add methods to: ```@example implementing -using OceanBioME: Biogeochemistry, UnderlyingBiogeochemicalModel +using OceanBioME: Biogeochemistry +using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry using Oceananigans.Units import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, @@ -29,7 +30,7 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, We then define our `struct` with the model parameters, as well as slots for the particles, light attenuation, and sediment models: ```@example implementing -@kwdef struct NutrientPhytoplankton{FT, W} <: UnderlyingBiogeochemicalModel +@kwdef struct NutrientPhytoplankton{FT, W} <: AbstractContinuousFormBiogeochemistry base_growth_rate :: FT = 1.27 / day # 1 / seconds nutrient_half_saturation :: FT = 0.025 * 1000 / 14 # mmol N / m³ light_half_saturation :: FT = 300.0 # micro einstein / m² / s From 3f6f2c37cea2e89906c379e3d06fa8e7c782729c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 3 Sep 2023 20:14:26 +1000 Subject: [PATCH 58/68] bit code cleanup --- docs/src/model_components/individuals/index.md | 4 ++-- examples/column.jl | 10 +++++----- examples/data_forced.jl | 12 ++++++------ 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index 95b39dc3b..be71f1e83 100644 --- a/docs/src/model_components/individuals/index.md +++ b/docs/src/model_components/individuals/index.md @@ -43,7 +43,7 @@ function update_lagrangian_particle_properties!(particles::GrowingParticles, mod return nothing end -nothing # hide +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: @@ -67,7 +67,7 @@ function update_tendencies!(bgc, particles::GrowingParticles, model) return nothing end -nothing # hide +nothing #hide ``` Now we can just plug this into any biogeochemical model setup to have particles (currently [NPZD](@ref NPZD) and [LOBSTER](@ref LOBSTER)): diff --git a/examples/column.jl b/examples/column.jl index b005d491f..0317fbc6e 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -54,7 +54,7 @@ biogeochemistry = LOBSTER(; grid, CO₂_flux = GasExchange(; gas = :CO₂, temperature = temp, salinity = (args...) -> 35) model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), + closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), )) @@ -68,10 +68,10 @@ set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2 simulation = Simulation(model, Δt = 3minutes, stop_time = 100days) progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", - iteration(sim), - prettytime(sim), - prettytime(sim.Δt), - prettytime(sim.run_wall_time)) + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10days)) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index bd7cf2145..5a57007ba 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -28,7 +28,7 @@ nothing #hide # Loading the forcing data from our online copy dd = DataDep( "example_data", - "example data from subpolar re analysis and observational products", + "example data from subpolar re analysis and observational products", "https://github.com/OceanBioME/OceanBioME_example_data/raw/main/subpolar.nc" ) register(dd) @@ -39,9 +39,9 @@ salinity = ncread(filename, "so") mld = ncread(filename, "mld") par = ncread(filename, "par") -temperature_itp = LinearInterpolation(times, temp) -salinity_itp = LinearInterpolation(times, salinity) -mld_itp = LinearInterpolation(times, mld) +temperature_itp = LinearInterpolation(times, temp) +salinity_itp = LinearInterpolation(times, salinity) +mld_itp = LinearInterpolation(times, mld) PAR_itp = LinearInterpolation(times, par) t_function(x, y, z, t) = temperature_itp(mod(t, 364days)) @@ -74,7 +74,7 @@ biogeochemistry = LOBSTER(; grid, CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = s_function) model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), + closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),)) @@ -133,7 +133,7 @@ x, y, z = nodes(P) times = P.times; # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and -# the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth +# the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth # that corresponds to `k = grid.Nz - 20`. air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) From 41be828536b5f670a003eaeabffa819c77620223 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 3 Sep 2023 20:18:19 +1000 Subject: [PATCH 59/68] avoid ; in literated examples --- examples/column.jl | 3 ++- examples/data_forced.jl | 3 ++- examples/kelp.jl | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/column.jl b/examples/column.jl index 0317fbc6e..30ac5b799 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -98,7 +98,8 @@ bPOM = FieldTimeSeries("$filename.jld2", "bPOM") Alk = FieldTimeSeries("$filename.jld2", "Alk") x, y, z = nodes(P) -times = P.times; +times = P.times +nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 5a57007ba..2f6fffb6f 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -130,7 +130,8 @@ bPOM = FieldTimeSeries("$filename.jld2", "bPOM") Alk = FieldTimeSeries("$filename.jld2", "Alk") x, y, z = nodes(P) -times = P.times; +times = P.times +nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth diff --git a/examples/kelp.jl b/examples/kelp.jl index 845a979c2..cb91ff6d6 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -121,7 +121,8 @@ bPOC = FieldTimeSeries("$filename.jld2", "bPOC") Alk = FieldTimeSeries("$filename.jld2", "Alk") x, y, z = nodes(P) -times = P.times; +times = P.times +nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth From 9a78a1b3fd7dab53871a1883a812ac9f2963abe7 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Sep 2023 12:05:30 +1000 Subject: [PATCH 60/68] update dependencies + bump compat for docs/Project --- Manifest.toml | 52 +++++++++++++++++++++-------------------------- docs/Project.toml | 2 +- 2 files changed, 24 insertions(+), 30 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index fc62c33f4..3a78e27c1 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -272,9 +272,9 @@ version = "1.3.1" [[deps.HDF5_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"] -git-tree-sha1 = "592e1c427983a465831fc73c5ae0ca5d0ac13a9e" +git-tree-sha1 = "10c72358aaaa5cd6bc7cc39b95e6eadf92f5a336" uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" -version = "1.14.1+0" +version = "1.14.2+0" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -345,15 +345,15 @@ version = "0.9.8" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "8695a49bfe05a2dc0feeefd06b4ca6361a018729" +git-tree-sha1 = "a9d2ce1d5007b1e8f6c5b89c5a31ff8bd146db5c" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "6.1.0" +version = "6.2.1" [[deps.LLVMExtra_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "c35203c1e1002747da220ffc3c0762ce7754b08c" +git-tree-sha1 = "7ca6850ae880cc99b59b88517545f91a52020afa" uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" -version = "0.0.23+0" +version = "0.0.25+0" [[deps.LLVMOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -388,10 +388,10 @@ version = "1.10.2+0" uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" [[deps.Libiconv_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.1+2" +version = "1.17.0+0" [[deps.LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] @@ -424,9 +424,9 @@ version = "2023.2.0+0" [[deps.MPI]] deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"] -git-tree-sha1 = "32cafbe56c7f0b7160a1a6c492773af66c0b722f" +git-tree-sha1 = "df53d0e1e0dbebf2315f4cd35e13e52ad43416c2" uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195" -version = "0.20.14" +version = "0.20.15" [deps.MPI.extensions] AMDGPUExt = "AMDGPU" @@ -500,9 +500,9 @@ version = "1.2.0" [[deps.Oceananigans]] deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] -git-tree-sha1 = "996eb159d0a1380d09beb7190ed0807d51ba0924" +git-tree-sha1 = "8746c619e15950de59160b2d99e61961cfa57fcc" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.87.1" +version = "0.87.2" [[deps.OffsetArrays]] deps = ["Adapt"] @@ -576,9 +576,9 @@ version = "1.9.2" [[deps.PkgVersion]] deps = ["Pkg"] -git-tree-sha1 = "f6cf8e7944e50901594838951729a1861e668cb8" +git-tree-sha1 = "f9501cc0430a26bc3d156ae1b5b0c1b47af4d6da" uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688" -version = "0.3.2" +version = "0.3.3" [[deps.PrecompileTools]] deps = ["Preferences"] @@ -697,12 +697,6 @@ git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" version = "1.1.1" -[[deps.SnoopPrecompile]] -deps = ["Preferences"] -git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" -uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" -version = "1.0.3" - [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -727,10 +721,10 @@ uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" version = "0.8.8" [[deps.StaticArrayInterface]] -deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "Static", "SuiteSparse"] -git-tree-sha1 = "33040351d2403b84afce74dae2e22d3f5b18edcb" +deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] +git-tree-sha1 = "03fec6800a986d191f64f5c0996b59ed526eda25" uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718" -version = "1.4.0" +version = "1.4.1" weakdeps = ["OffsetArrays", "StaticArrays"] [deps.StaticArrayInterface.extensions] @@ -764,9 +758,9 @@ version = "1.9.0" [[deps.StatsAPI]] deps = ["LinearAlgebra"] -git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7" +git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.6.0" +version = "1.7.0" [[deps.Strided]] deps = ["LinearAlgebra", "StridedViews", "TupleTools"] @@ -880,10 +874,10 @@ uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" version = "1.3.0" [[deps.XML2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] +git-tree-sha1 = "04a51d15436a572301b5abbb9d099713327e9fc4" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.10.3+0" +version = "2.10.4+0" [[deps.Zlib_jll]] deps = ["Libdl"] diff --git a/docs/Project.toml b/docs/Project.toml index a6258f013..a85e13173 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -17,4 +17,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" CairoMakie = "0.10" DocumenterCitations = "1" Literate = "≥2.9.0" -Oceananigans = "0.84.1, 0.85" +Oceananigans = "0.84.1, 0.85, 0.86, 0.87" From af43c8cd21417a5a0d9c3e9a79ab04faf22fba64 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 7 Sep 2023 10:24:15 +0100 Subject: [PATCH 61/68] Apply suggestions from code review Co-authored-by: Navid C. Constantinou --- Project.toml | 2 +- docs/src/model_components/biogeochemical/index.md | 2 +- docs/src/model_components/utils.md | 2 +- examples/column.jl | 4 ++-- examples/data_forced.jl | 4 ++-- examples/eady.jl | 3 ++- src/Boundaries/Sediments/Sediments.jl | 2 +- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 2 +- src/OceanBioME.jl | 1 + src/Utils/negative_tracers.jl | 5 +++-- 10 files changed, 15 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index cd1ca939c..4761ddc56 100644 --- a/Project.toml +++ b/Project.toml @@ -13,7 +13,7 @@ StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" [compat] Adapt = "3" -JLD2 = "^0.4" +JLD2 = "0.4" KernelAbstractions = "0.9" Oceananigans = "0.84.1, 0.85, 0.86, 0.87" Roots = "2" diff --git a/docs/src/model_components/biogeochemical/index.md b/docs/src/model_components/biogeochemical/index.md index 191e58994..2e28f7f0a 100644 --- a/docs/src/model_components/biogeochemical/index.md +++ b/docs/src/model_components/biogeochemical/index.md @@ -21,5 +21,5 @@ This will set up the required tracers and auxiliary fields, and you may also set Models usually have a default [light attenuation model](@ref light) specified, these may be substituted easily by passing different models as parameters as above. Our models are implemented in an abstract framework `Biogeochemistry` which contains `underlying_biogeochemistry`, `light_attenuation`, `sediment`, and `modifiers`. -This is automatically setup for existing BGC models, but may also be used to couple any BGC model with light attenuation and sediments. +This is automatically set up for existing BGC models, but may also be used to couple any BGC model with light attenuation and sediments. See the [implementation page](@ref model_implementation) for some more information on how to couple other models. \ No newline at end of file diff --git a/docs/src/model_components/utils.md b/docs/src/model_components/utils.md index 01dba5a39..1adfd85dc 100644 --- a/docs/src/model_components/utils.md +++ b/docs/src/model_components/utils.md @@ -35,4 +35,4 @@ Here you should carefully consider which tracers form a conserved group (if at a negativity_protection = ZeroNegativeTracers() biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` -The latter optionally takes a named tuple of parameters which may include `exclude` which can be a tuple of tracer names (Symbols) which are allowed to be negative. +The latter optionally takes a named tuple of parameters that may include `exclude`, which can be a tuple of tracer names (Symbols) which are allowed to be negative. diff --git a/examples/column.jl b/examples/column.jl index 30ac5b799..44acb5d6a 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -42,9 +42,9 @@ grid = RectilinearGrid(size=(1, 1, 50), extent=(20meters, 20meters, 200meters)) # ## Model -# First we will define the biogeochemical model including carbonate chemistry +# First we define the biogeochemical model including carbonate chemistry # and scaling of negative tracers(see discussion in the [positivity preservation](@ref pos-preservation)) -# and then setup the Oceananigans model with the boundary condition for the DIC based on the air-sea CO₂ flux +# and then setup the Oceananigans model with the boundary condition for the DIC based on the air-sea CO₂ flux. biogeochemistry = LOBSTER(; grid, surface_phytosynthetically_active_radiation = PAR⁰, diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 2f6fffb6f..5f7de5339 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -28,7 +28,7 @@ nothing #hide # Loading the forcing data from our online copy dd = DataDep( "example_data", - "example data from subpolar re analysis and observational products", + "example data from subpolar re-analysis and observational products", "https://github.com/OceanBioME/OceanBioME_example_data/raw/main/subpolar.nc" ) register(dd) @@ -134,7 +134,7 @@ times = P.times nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and -# the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth +# the carbon export by computing how much carbon sinks below some arbitrary depth; here we use depth # that corresponds to `k = grid.Nz - 20`. air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) diff --git a/examples/eady.jl b/examples/eady.jl index e528f0968..18f06783e 100644 --- a/examples/eady.jl +++ b/examples/eady.jl @@ -134,7 +134,8 @@ DIC = FieldTimeSeries("eady_turbulence_bgc.jld2", "DIC") times = ζ.times xζ, yζ, zζ = nodes(ζ) -xc, yc, zc = nodes(P); +xc, yc, zc = nodes(P) +nothing #hide # and plot. diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index aff6af773..d1fcc9767 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -11,7 +11,7 @@ using Oceananigans.Architectures: device using Oceananigans.Utils: launch! using Oceananigans.Advection: advective_tracer_flux_z using Oceananigans.Units: day -using Oceananigans.Fields: CenterField, Face, Center, ConstantField +using Oceananigans.Fields: ConstantField using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity using Oceananigans.Grids: zspacing using Oceananigans.Operators: volume diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 9c0a913a7..69fdbba6c 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -167,7 +167,7 @@ struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry zooplankton_calcite_dissolution, optionals, - + sinking_velocities) end end diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index e7093847f..a7b872642 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -139,6 +139,7 @@ Returns the redfield ratio of `tracer_name` from `bgc` at `i`, `j`, `k`. @inline redfield(i, j, k, val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc, tracers) """ + redfield(val_tracer_name, bgc) redfield(val_tracer_name, bgc, tracers) Returns the redfield ratio of `tracer_name` from `bgc` when it is constant across the domain. diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 11839fd62..a49229adc 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -10,7 +10,8 @@ import Oceananigans.Biogeochemistry: update_tendencies!, update_biogeochemical_s """ ZeroNegativeTracers(; exclude = ()) -Construct a modifier that zeros any negative tracers excluding those listed in `exclude`. + +Construct a modifier that zeroes any negative tracers excluding those listed in `exclude`. !!! danger "Tracer conservation" This method is _not_ recommended as a way to preserve positivity of tracers since @@ -72,7 +73,7 @@ end """ ScaleNegativeTracers(model::UnderlyingBiogeochemicalModel; warn = false) -Constructs a modifier to scale the conserved tracers in `model`. +Construct a modifier to scale the conserved tracers in `model`. If `warn` is true then scaling will raise a warning. """ From 89c55a5db20c72f611117bb7363bbfe2c339eea0 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 10:25:49 +0100 Subject: [PATCH 62/68] applied review suggesitons --- examples/data_forced.jl | 2 -- src/Light/Light.jl | 1 - 2 files changed, 3 deletions(-) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 5f7de5339..506befc75 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -104,8 +104,6 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, schedule = TimeInterval(1day), overwrite_existing = true) -# TODO: make tendency callback to force no NaNs in tendencies - wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 1.5, min_change = 0.75, cell_diffusion_timescale = column_diffusion_timescale, diff --git a/src/Light/Light.jl b/src/Light/Light.jl index 03addbb5e..2b36457da 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -26,7 +26,6 @@ import Base: show, summary import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields, update_biogeochemical_state!, required_biogeochemical_auxiliary_fields import Oceananigans.BoundaryConditions: _fill_top_halo! -# Fallback include("2band.jl") include("morel.jl") From ad12a2187c88c1af9a34035c6a26e2e12f9fa504 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 10:28:50 +0100 Subject: [PATCH 63/68] moved `Pkg.add`s out of `@example` blocks in quick start to see if that fixes the docs build issue --- docs/src/quick_start.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index 54eba8114..6f6f8bd69 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -3,7 +3,7 @@ OceanBioME provides biogeochemical models to plug into [Oceananigans](https://github.com/CliMA/Oceananigans.jl), for example this code will run one month of a single column, 7 variable (P, Z, sPOM, bPOM, DOM, NO₃, NH₄) biogeochemical situation with constant forcing. First we need to check we have the required dependencies: -```@example quickstart +```julia using Pkg Pkg.add(["OceanBioME", "Oceananigans"]) ``` From 78093d798570062a69fdc377dd9089a728d0fbda Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 13:12:01 +0100 Subject: [PATCH 64/68] fixed quick start --- docs/src/quick_start.md | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index 6f6f8bd69..854acb85d 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -30,13 +30,17 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, run!(simulation) ``` -This isn't quite as simple as it could be as it records the output so that we can visualize it: +This isn't quite as simple as it could be as it records the output so that we can visualize it, first check the required packages are installed: -```@example quickstart +```julia Pkg.add("CairoMakie") using CairoMakie +``` +and then load the data and plot: + +```@example quickstart phytoplankton = FieldTimeSeries("quickstart.jld2", "P") nitrates = FieldTimeSeries("quickstart.jld2", "NO₃") From cdbd356c0ead16dfce353f8dce8458f2b7bf18b6 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 13:51:08 +0100 Subject: [PATCH 65/68] rephrasing in quick start --- docs/src/quick_start.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index 854acb85d..aa82f893a 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -30,7 +30,7 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, run!(simulation) ``` -This isn't quite as simple as it could be as it records the output so that we can visualize it, first check the required packages are installed: +We can then visualize it, first check the required packages are installed: ```julia Pkg.add("CairoMakie") From 0cefee54d41156b1f4b27ba74666403d1df714d9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Sep 2023 16:43:10 +0300 Subject: [PATCH 66/68] Update quick_start.md --- docs/src/quick_start.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index aa82f893a..dbd6e7ee3 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -46,6 +46,8 @@ nitrates = FieldTimeSeries("quickstart.jld2", "NO₃") _, _, z = nodes(nitrates) +using CairoMakie + fig = Figure() axis_kwargs = (xlabel = "Day", ylabel = "Depth (m)") From a2570ed6df7ae72f9a9389e5126b3b3ba8342b05 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Sep 2023 16:44:50 +0300 Subject: [PATCH 67/68] Update quick_start.md --- docs/src/quick_start.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index dbd6e7ee3..8a28c4785 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -34,8 +34,6 @@ We can then visualize it, first check the required packages are installed: ```julia Pkg.add("CairoMakie") - -using CairoMakie ``` and then load the data and plot: From 8f2b56ec8f2837e136499e4abfd5cf4e24cb3b23 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 14:45:33 +0100 Subject: [PATCH 68/68] oops --- docs/src/quick_start.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index 8a28c4785..b27e43a7a 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -39,13 +39,13 @@ Pkg.add("CairoMakie") and then load the data and plot: ```@example quickstart +using CairoMakie + phytoplankton = FieldTimeSeries("quickstart.jld2", "P") nitrates = FieldTimeSeries("quickstart.jld2", "NO₃") _, _, z = nodes(nitrates) -using CairoMakie - fig = Figure() axis_kwargs = (xlabel = "Day", ylabel = "Depth (m)")