Merge branch 'main' into main

README.md

@@ -46,9 +46,9 @@ model = NonhydrostaticModel(; grid, biogeochemistry,
 
 @inline front(x, z, μ, δ) = μ + δ * tanh((x - 7000 + 4 * z) / 500)
 
-Pᵢ(x, y, z) = ifelse(z > -50, 0.03, 0.01) 
-Nᵢ(x, y, z) = front(x, z, 2.5, -2)
-Tᵢ(x, y, z) = front(x, z, 9, 0.05)
+Pᵢ(x, z) = ifelse(z > -50, 0.03, 0.01)
+Nᵢ(x, z) = front(x, z, 2.5, -2)
+Tᵢ(x, z) = front(x, z, 9, 0.05)
 
 set!(model, N = Nᵢ, P = Pᵢ, Z = Pᵢ, T = Tᵢ)
 

docs/make.jl

@@ -59,10 +59,10 @@ model_parameters = (LOBSTER(; grid = BoxModelGrid()),
                     GasExchange(; gas = :CO₂).condition.func,
                     GasExchange(; gas = :O₂).condition.func)
 
-gas_exchange_gas(::Val{G}) where G = G
+exchanged_gas(::Val{G}) where G = G
 
 model_name(model) = if Base.typename(typeof(model)).wrapper == GasExchange
-                        "$(gas_exchange_gas(model.gas)) air-sea exchange"
+                        "$(exchanged_gas(model.gas)) air-sea exchange"
                     else
                         Base.typename(typeof(model)).wrapper
                     end

docs/src/index.md

@@ -44,9 +44,9 @@ model = NonhydrostaticModel(; grid, biogeochemistry,
 
 @inline front(x, z, μ, δ) = μ + δ * tanh((x - 7000 + 4 * z) / 500)
 
-Pᵢ(x, y, z) = ifelse(z > -50, 0.03, 0.01) 
-Nᵢ(x, y, z) = front(x, z, 2.5, -2)
-Tᵢ(x, y, z) = front(x, z, 9, 0.05)
+Pᵢ(x, z) = ifelse(z > -50, 0.03, 0.01)
+Nᵢ(x, z) = front(x, z, 2.5, -2)
+Tᵢ(x, z) = front(x, z, 9, 0.05)
 
 set!(model, N = Nᵢ, P = Pᵢ, Z = Pᵢ, T = Tᵢ)
 

docs/src/model_components/biogeochemical/LOBSTER.md

@@ -12,17 +12,11 @@ julia> using OceanBioME, Oceananigans
 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) 
- 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 
+LOBSTER{Float64} with carbonates ✅, oxygen ❌, variable Redfield ratio ❌and (:sPOM, :bPOM) sinking 
  Light attenuation: Two-band light attenuation model (Float64)
  Sediment: Nothing
  Particles: Nothing
+ Modifiers: Nothing
 ```
 
 ## Model equations

docs/src/model_components/individuals/index.md

@@ -54,16 +54,22 @@ using OceanBioME.Particles: get_node
 
 function update_tendencies!(bgc, particles::GrowingParticles, model)
     @inbounds for p in 1:length(particles)
-        # here we use an OceanBioME utility to find the nearest nodes to apply the tendency to
-        nodes, normfactor = @inbounds get_nearest_nodes(particles.x[p], particles.y[p], particles.z[p], model.grid, (Center(), Center(), Center()))
+        i, j, k = fractional_indices(x, y, z, (Center(), Center(), Center()), grid)
 
-        for (i, j, k, d) in nodes 
-            # Reflect back on Bounded boundaries or wrap around for Periodic boundaries
-            i, j, k = (get_node(TX(), i, grid.Nx), get_node(TY(), j, grid.Ny), get_node(TZ(), k, grid.Nz))
+        # Convert fractional indices to unit cell coordinates 0 ≤ (ξ, η, ζ) ≤ 1
+        # and integer indices (with 0-based indexing).
+        ξ, i = modf(i)
+        η, j = modf(j)
+        ζ, k = modf(k)
 
-            node_volume = volume(i, j, k, grid, LX(), LY(), LZ())
-            @inbounds model.timestepper.Gⁿ.NO₃[i, j, k] += particles.nitrate_uptake[p] / (d * node_volume)
-        end
+        # Round to nearest node and enforce boundary conditions
+        i, j, k = (get_node(TX(), Int(ifelse(ξ < 0.5, i + 1, i + 2)), grid.Nx),
+                   get_node(TY(), Int(ifelse(η < 0.5, j + 1, j + 2)), grid.Ny),
+                   get_node(TZ(), Int(ifelse(ζ < 0.5, k + 1, k + 2)), grid.Nz))
+
+        node_volume = volume(i, j, k, grid, Center(), Center(), Center())
+
+        model.timestepper.Gⁿ.NO₃[i, j, k] += particles.nitrate_uptake[p] / (d * node_volume)
     end
     return nothing
 end

docs/src/model_implementation.md

@@ -233,11 +233,11 @@ using OceanBioME.Sediments: InstantRemineralisation
 
 # define some simple forcing
 
-@inline surface_PAR(x, y, t) = 200 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
+@inline surface_PAR(t) = 200 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
 
-@inline temp(x, y, z, t) = cos(t * 2π / year + 50days) + 28
+@inline ∂ₜT(z, t) = - 2π / year * sin(t * 2π / year + 50days)
 
-@inline κₜ(x, y, z, t) = 1e-2 * (1 + tanh((z - 50) / 10)) / 2 + 1e-4
+@inline κₜ(z, t) = 1e-2 * (1 + tanh((z - 50) / 10)) / 2 + 1e-4
 
 # define the grid
 
@@ -267,7 +267,7 @@ biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; sinking_velocity);
 model = NonhydrostaticModel(; grid,
                               biogeochemistry,
                               closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), 
-                              forcing = (T = Relaxation(rate = 1/day, target = temp), ))
+                              forcing = (; T = ∂ₜT))
 
 set!(model, P = 0.01, N = 15, T = 28)
 

docs/src/quick_start.md

@@ -14,9 +14,7 @@ using Oceananigans.Units
 
 grid = RectilinearGrid(size = 10, extent = 200meters, topology = (Flat, Flat, Bounded))
 
-PAR = CenterField(grid)
-
-model = NonhydrostaticModel(; grid, biogeochemistry = LOBSTER(; grid), auxiliary_fields = (; PAR))
+model = NonhydrostaticModel(; grid, biogeochemistry = LOBSTER(; grid))
 
 set!(model, P = 0.001, Z = 0.001, NO₃ = 1, NH₄ = 0.01)
 

src/Light/2band.jl

@@ -1,9 +1,13 @@
 @kernel function update_TwoBandPhotosyntheticallyActiveRadiation!(PAR, grid, P, surface_PAR, t, PAR_model) 
     i, j = @index(Global, NTuple)
 
-    x, y, _ = node(i, j, 1, grid, Center(), Center(), Center())
+    LX, LY, _ = location(P)
+
+    k, k′ = domain_boundary_indices(RightBoundary(), grid.Nz)
+
+    X = z_boundary_node(i, j, k′, grid, LX(), LY())
 
-    PAR⁰ = surface_PAR(x, y, t)
+    PAR⁰ = surface_PAR(X..., t)
 
     kʳ = PAR_model.water_red_attenuation
     kᵇ = PAR_model.water_blue_attenuation
@@ -105,8 +109,7 @@ function TwoBandPhotosyntheticallyActiveRadiation(; grid,
 
     field = CenterField(grid; boundary_conditions = 
                             regularize_field_boundary_conditions(
-                            FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR)),
-                            grid, :PAR))
+                                FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR)), grid, :PAR))
 
     return TwoBandPhotosyntheticallyActiveRadiation(water_red_attenuation,
                                                     water_blue_attenuation,

src/Light/Light.jl

@@ -5,19 +5,21 @@ module Light
 
 export TwoBandPhotosyntheticallyActiveRadiation, update_PAR!
 
-using KernelAbstractions
+using KernelAbstractions, Oceananigans.Units
 using KernelAbstractions.Extras.LoopInfo: @unroll
 using Oceananigans.Architectures: device, architecture
 using Oceananigans.Utils: launch!
 using Oceananigans: Center, Face, fields
 using Oceananigans.Grids: node, znodes
-using Oceananigans.Fields: CenterField, TracerFields
+using Oceananigans.Fields: CenterField, TracerFields, location
 using Oceananigans.BoundaryConditions: fill_halo_regions!, 
                                        ValueBoundaryCondition, 
                                        FieldBoundaryConditions, 
                                        regularize_field_boundary_conditions, 
-                                       ContinuousBoundaryFunction
-using Oceananigans.Units
+                                       ContinuousBoundaryFunction,
+                                       z_boundary_node,
+                                       domain_boundary_indices,
+                                       RightBoundary
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid
 
 import Adapt: adapt_structure, adapt
@@ -29,4 +31,8 @@ import Oceananigans.BoundaryConditions: _fill_top_halo!
 include("2band.jl")
 include("morel.jl")
 
+default_surface_PAR(x, y, t) = default_surface_PAR(t)
+default_surface_PAR(x_or_y, t) = default_surface_PAR(t)
+default_surface_PAR(t) = 100 * max(0, cos(t * π / 12hours))
+
 end

src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl

@@ -46,7 +46,7 @@ export LOBSTER
 using Oceananigans.Units
 using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField
 
-using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation
+using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR
 using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, ScaleNegativeTracers
 using OceanBioME.BoxModels: BoxModel
 using OceanBioME.Boundaries.Sediments: sinking_flux
@@ -204,7 +204,7 @@ end
               dissolved_organic_breakdown_rate::FT = 3.86e-7, # 1/s
               zooplankton_calcite_dissolution::FT = 0.3,
 
-              surface_phytosynthetically_active_radiation::SPAR = (x, y, t) -> 100*max(0.0, cos(t*π/(12hours))),
+              surface_phytosynthetically_active_radiation::SPAR = default_surface_PAR,
 
               light_attenuation_model::LA =
                   TwoBandPhotosyntheticallyActiveRadiation(; grid,
@@ -246,17 +246,11 @@ julia> using Oceananigans
 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) 
- 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 
+LOBSTER{Float64} with carbonates ❌, oxygen ❌, variable Redfield ratio ❌and (:sPOM, :bPOM) sinking 
  Light attenuation: Two-band light attenuation model (Float64)
  Sediment: Nothing
  Particles: Nothing
+ Modifiers: Nothing
 ```
 """
 function LOBSTER(; grid,
@@ -290,7 +284,7 @@ function LOBSTER(; grid,
                    dissolved_organic_breakdown_rate::FT = 3.86e-7, # 1/s
                    zooplankton_calcite_dissolution::FT = 0.3,
 
-                   surface_phytosynthetically_active_radiation = (x, y, t) -> 100 * max(0.0, cos(t * π / (12hours))),
+                   surface_phytosynthetically_active_radiation = default_surface_PAR,
 
                    light_attenuation_model::LA =
                        TwoBandPhotosyntheticallyActiveRadiation(; grid, 
@@ -353,7 +347,13 @@ function LOBSTER(; grid,
 
     if scale_negatives
         scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid)
-        modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler)
+        if isnothing(modifiers)
+            modifiers = scaler
+        elseif modifiers isa Tuple
+            modifiers = (modifiers..., scaler)
+        else
+            modifiers = (modifiers, scaler)
+        end
     end
 
     return Biogeochemistry(underlying_biogeochemistry;
@@ -432,14 +432,14 @@ adapt_structure(to, lobster::LOBSTER) =
             adapt(to, lobster.optionals),
             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)")
+summary(::LOBSTER{FT, Val{B}, NamedTuple{K, V}}) where {FT, B, K, V} = string("LOBSTER{$FT} with carbonates $(B[1] ? :✅ : :❌), oxygen $(B[2] ? :✅ : :❌), variable Redfield ratio $(B[3] ? :✅ : :❌)and $K sinking")
 
-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))
+show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W}  = print(io, string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) 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)
 

src/Models/AdvectedPopulations/NPZD.jl

@@ -22,7 +22,7 @@ using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry
 using Oceananigans.Units
 using Oceananigans.Fields: ZeroField
 
-using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation
+using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR
 using OceanBioME: setup_velocity_fields, show_sinking_velocities
 using OceanBioME.BoxModels: BoxModel
 using OceanBioME.Boundaries.Sediments: sinking_flux
@@ -109,7 +109,7 @@ end
                                                zoo_base_mortality_rate::FT = 0.3395 / day, # 1/s/(mmol N / m³)²
                                                remineralization_rate::FT = 0.1213 / day, # 1/s
 
-                                               surface_phytosynthetically_active_radiation = (x, y, t) -> 100 * max(0.0, cos(t * π / 12hours)),
+                                               surface_phytosynthetically_active_radiation = default_surface_PAR,
                                                light_attenuation_model::LA =
                                                    TwoBandPhotosyntheticallyActiveRadiation(; grid,
                                                                                               surface_PAR = surface_phytosynthetically_active_radiation),
@@ -145,14 +145,11 @@ julia> using Oceananigans
 julia> grid = RectilinearGrid(size=(20, 30), extent=(200, 200), topology=(Bounded, Flat, Bounded));
 
 julia> model = NutrientPhytoplanktonZooplanktonDetritus(; grid)
-Nutrient Phytoplankton Zooplankton Detritus model (Float64) 
- Sinking Velocities:
-    ├── P: 0.0 to -2.9525462962962963e-6 m/s 
-    └── D: 0.0 to -3.181597222222222e-5 m/s 
+NutrientPhytoplanktonZooplanktonDetritus{Float64} model, with (:P, :D) sinking 
  Light attenuation: Two-band light attenuation model (Float64)
  Sediment: Nothing
  Particles: Nothing
-
+ Modifiers: Nothing
 ```
 """
 function NutrientPhytoplanktonZooplanktonDetritus(; grid,
@@ -168,7 +165,7 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid,
                                                     zoo_base_mortality_rate::FT = 0.3395 / day, # 1/s/(mmol N / m³)²
                                                     remineralization_rate::FT = 0.1213 / day, # 1/s
 
-                                                    surface_phytosynthetically_active_radiation = (x, y, t) -> 100 * max(0.0, cos(t * π / 12hours)),
+                                                    surface_phytosynthetically_active_radiation = default_surface_PAR,
                                                     light_attenuation_model::LA =
                                                         TwoBandPhotosyntheticallyActiveRadiation(; grid,
                                                                                                    surface_PAR = surface_phytosynthetically_active_radiation),
@@ -297,9 +294,9 @@ function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:NPZD}, <:Any
     getproperty(model.values, :T) .= model.forcing.T(model.clock.time)
 end
 
-summary(::NPZD{FT, W}) where {FT, W} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)")
-show(io::IO, model::NPZD) = string(summary(model), " \n",
-                                   " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities))
+summary(::NPZD{FT, NamedTuple{K, V}}) where {FT, K, V} = string("NutrientPhytoplanktonZooplanktonDetritus{$FT} model, with $K sinking")
+show(io::IO, model::NPZD{FT}) where {FT} = print(io, string("NutrientPhytoplanktonZooplanktonDetritus{$FT} model \n",
+                                                            "└── Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)))
 
 @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w)
 

src/OceanBioME.jl

@@ -163,10 +163,13 @@ conserved_tracers(model::Biogeochemistry) = conserved_tracers(model.underlying_b
 
 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, summary(model.underlying_biogeochemistry), " \n",
                 " Light attenuation: ", summary(model.light_attenuation), "\n",
                 " Sediment: ", summary(model.sediment), "\n",
-                " Particles: ", summary(model.particles))
+                " Particles: ", summary(model.particles), "\n",
+                " Modifiers: ", summary(model.modifiers))
+
+summary(modifiers::Tuple) = tuple([summary(modifier) for modifier in modifiers])
 
 include("Utils/Utils.jl")
 include("Boundaries/Boundaries.jl")

src/Utils/negative_tracers.jl

@@ -6,6 +6,7 @@ using Oceananigans.Architectures: device, architecture, arch_array
 using Oceananigans.Biogeochemistry: AbstractBiogeochemistry
 
 import Adapt: adapt_structure, adapt
+import Base: summary, show
 import Oceananigans.Biogeochemistry: update_tendencies!, update_biogeochemical_state!
 
 """
@@ -84,6 +85,10 @@ function ScaleNegativeTracers(bgc::AbstractBiogeochemistry, grid; warn = false)
     return ScaleNegativeTracers(tracers, scalefactors, warn)
 end
 
+summary(scaler::ScaleNegativeTracers) = string("Mass conserving negative scaling of $(scaler.tracers)")
+show(io::IO, scaler::ScaleNegativeTracers) = print(io, string(summary(scaler), "\n",
+                                                          "└── Scalefactors: $(scaler.scalefactors)"))
+
 function update_biogeochemical_state!(model, scale::ScaleNegativeTracers)
     workgroup, worksize = work_layout(model.grid, :xyz)