Skip to content

Commit

Permalink
Merge pull request #127 from OceanBioME/jsw/sediment_with_bathymetry
Browse files Browse the repository at this point in the history
  • Loading branch information
@jagoosw
jagoosw authored Aug 29, 2023
2 parents ef69e24 + 3a352ad commit ed7ec50
Show file tree
Hide file tree
Showing 10 changed files with 172 additions and 80 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "OceanBioME"
uuid = "a49af516-9db8-4be4-be45-1dad61c5a376"
authors = ["Jago Strong-Wright <[email protected]>", "John R Taylor <[email protected]>", "Si Chen <[email protected]>"]
version = "0.4.1"
version = "0.5.0"

[deps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
Expand Down
6 changes: 3 additions & 3 deletions docs/src/model_implementation.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -225,10 +225,10 @@ using OceanBioME.Boundaries.Sediments: sinking_flux
import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver
@inline nitrogen_flux(grid, advection, bgc::NutrientPhytoplankton, tracers, i, j) =
sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers)
@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(bgc, tracers, i, j) * 6.56
@inline carbon_flux(bgc::NutrientPhytoplankton, tracers, i, j) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56
@inline remineralisation_receiver(::NutrientPhytoplankton) = :N
```
Expand Down
47 changes: 39 additions & 8 deletions src/Boundaries/Sediments/Sediments.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ module Sediments
export SimpleMultiG, InstantRemineralisation

using KernelAbstractions
using OceanBioME: ContinuousFormBiogeochemistry
using OceanBioME: ContinuousFormBiogeochemistry, BoxModelGrid
using Oceananigans
using Oceananigans.Architectures: device
using Oceananigans.Utils: launch!
Expand All @@ -13,13 +13,16 @@ using Oceananigans.Fields: CenterField, Face
using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity
using Oceananigans.Grids: zspacing
using Oceananigans.Operators: volume
using Oceananigans.Fields: Center
using Oceananigans.Fields: Center, ConstantField
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, immersed_cell

import Adapt: adapt_structure, adapt

abstract type AbstractSediment end
abstract type FlatSediment <: AbstractSediment end

@inline bottom_index(i, j, sediment) = @inbounds bottom_index_array(sediment)[i, j]

sediment_fields(::AbstractSediment) = ()

@inline update_sediment_tendencies!(bgc, sediment, model) = nothing
Expand All @@ -43,14 +46,42 @@ end
@inbounds G⁻[i, j, 1] = G⁰[i, j, 1]
end

@inline nitrogen_flux() = 0
@inline carbon_flux() = 0
@inline remineralisation_receiver() = nothing

@inline advection_scheme(advection, val_tracer) = advection
@inline advection_scheme(advection::NamedTuple, val_tracer::Val{T}) where T = advection[T]

@inline function sinking_flux(i, j, k, grid, advection, val_tracer::Val{T}, bgc, tracers) where T
return - advective_tracer_flux_z(i, j, k, grid, advection_scheme(advection, val_tracer), biogeochemical_drift_velocity(bgc, val_tracer).w, tracers[T]) /
volume(i, j, k, grid, Center(), Center(), Center())
end

calculate_bottom_indices(::BoxModelGrid) = 1
calculate_bottom_indices(grid) = ones(Int, size(grid)[1:2]...)

@kernel function find_bottom_cell(grid, bottom_indices)
i, j = @index(Global, NTuple)

Nz = size(grid, 3)

@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
k_bottom = 1

@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]) /
volume(i, j, 1, grid, Center(), Center(), Center())
while immersed_cell(i, j, k_bottom, grid.underlying_grid, grid.immersed_boundary) && k_bottom < Nz
k_bottom += 1
end

@inbounds bottom_indices[i, j] = k_bottom
end

function calculate_bottom_indices(grid::ImmersedBoundaryGrid)
indices = zeros(Int, size(grid)[1:2]...)

launch!(grid.architecture, grid, :xy, find_bottom_cell, grid, indices)

return indices
end

include("coupled_timesteppers.jl")
include("simple_multi_G.jl")
Expand Down
36 changes: 23 additions & 13 deletions src/Boundaries/Sediments/instant_remineralization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,24 +7,27 @@ becoming N, and a fraction being perminantly burried.
Burial efficiency from [RemineralisationFraction](@citet).
"""
struct InstantRemineralisation{FT, F, TE} <: FlatSediment
struct InstantRemineralisation{FT, F, TE, B} <: FlatSediment
burial_efficiency_constant1 :: FT
burial_efficiency_constant2 :: FT
burial_efficiency_half_saturaiton :: FT

fields :: F
tendencies :: TE
bottom_indices :: B

InstantRemineralisation(burial_efficiency_constant1::FT,
burial_efficiency_constant2::FT,
burial_efficiency_half_saturaiton::FT,
fields::F,
tendencies::TE) where {FT, F, TE} =
new{FT, F, TE}(burial_efficiency_constant1,
burial_efficiency_constant2,
burial_efficiency_half_saturaiton,
fields,
tendencies)
tendencies::TE,
bottom_indices::B) where {FT, F, TE, B} =
new{FT, F, TE, B}(burial_efficiency_constant1,
burial_efficiency_constant2,
burial_efficiency_half_saturaiton,
fields,
tendencies,
bottom_indices)
end

"""
Expand Down Expand Up @@ -56,15 +59,18 @@ function InstantRemineralisation(; grid,
tracer_names = (:N_storage, )

# add field slicing back ( indices = (:, :, 1)) when output writer can cope
fields = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names))
tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)),
G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)))
fields = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names))
tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)),
G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)))

bottom_indices = calculate_bottom_indices(grid)

return InstantRemineralisation(burial_efficiency_constant1,
burial_efficiency_constant2,
burial_efficiency_half_saturaiton,
fields,
tendencies)
tendencies,
bottom_indices)
end

adapt_structure(to, sediment::InstantRemineralisation) =
Expand All @@ -77,12 +83,16 @@ adapt_structure(to, sediment::InstantRemineralisation) =
sediment_tracers(::InstantRemineralisation) = (:N_storage, )
sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, )

@inline bottom_index_array(sediment::InstantRemineralisation) = sediment.bottom_indices

@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, timestepper)
i, j = @index(Global, NTuple)

Δz = zspacing(i, j, 1, grid, Center(), Center(), Center())
k = bottom_index(i, j, sediment)

Δz = zspacing(i, j, k, grid, Center(), Center(), Center())

flux = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz
flux = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz

burial_efficiency = sediment.burial_efficiency_constant1 + sediment.burial_efficiency_constant2 * ((flux * 6.56) / (7 + flux * 6.56)) ^ 2

Expand Down
43 changes: 26 additions & 17 deletions src/Boundaries/Sediments/simple_multi_G.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ import Base: show, summary
Hold the parameters and fields for a simple "multi G" single-layer sediment model.
Based on the Level 3 model described by [Soetaert2000](@citet).
"""
struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE} <: FlatSediment
struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B} <: FlatSediment
fast_decay_rate :: FT
slow_decay_rate :: FT

Expand All @@ -24,6 +24,7 @@ struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE} <: FlatSediment

fields :: F
tendencies :: TE
bottom_indices :: B
end

"""
Expand Down Expand Up @@ -108,22 +109,26 @@ function SimpleMultiG(; grid,
tracer_names = (:C_slow, :C_fast, :N_slow, :N_fast, :C_ref, :N_ref)

# add field slicing back ( indices = (:, :, 1)) when output writer can cope
fields = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names))
tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)),
G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)))
fields = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names))
tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)),
G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)))

bottom_indices = calculate_bottom_indices(grid)

F = typeof(fields)
TE = typeof(tendencies)

return SimpleMultiG{FT, P1, P2, P3, P4, F, TE}(fast_decay_rate, slow_decay_rate,
fast_redfield, slow_redfield,
fast_fraction, slow_fraction, refactory_fraction,
nitrate_oxidation_params,
denitrification_params,
anoxic_params,
solid_dep_params,
fields,
tendencies)
B = typeof(bottom_indices)

return SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B}(fast_decay_rate, slow_decay_rate,
fast_redfield, slow_redfield,
fast_fraction, slow_fraction, refactory_fraction,
nitrate_oxidation_params,
denitrification_params,
anoxic_params,
solid_dep_params,
fields,
tendencies,
bottom_indices)
end

adapt_structure(to, sediment::SimpleMultiG) =
Expand All @@ -149,16 +154,20 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow,
C_ref = model.fields.C_ref,
N_ref = model.fields.N_ref)

@inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices

@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, timestepper)
i, j = @index(Global, NTuple)

Δz = zspacing(i, j, 1, grid, Center(), Center(), Center())
k = bottom_index(i, j, sediment)

Δz = zspacing(i, j, k, grid, Center(), Center(), Center())

@inbounds begin

carbon_deposition = carbon_flux(grid, advection, bgc, tracers, i, j) * Δz
carbon_deposition = carbon_flux(i, j, k, grid, advection, bgc, tracers) * Δz

nitrogen_deposition = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz
nitrogen_deposition = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz

# rates
C_min_slow = sediment.fields.C_slow[i, j, 1] * sediment.slow_decay_rate
Expand Down
4 changes: 2 additions & 2 deletions src/Light/2band.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,8 +105,8 @@ 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,
Expand Down
1 change: 1 addition & 0 deletions src/Light/Light.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ using Oceananigans.BoundaryConditions: fill_halo_regions!,
regularize_field_boundary_conditions,
ContinuousBoundaryFunction
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid

import Adapt: adapt_structure, adapt
import Base: show, summary
Expand Down
24 changes: 14 additions & 10 deletions src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,6 +321,10 @@ function LOBSTER(; grid,

particles::P = nothing) where {FT, LA, S, P}

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"
end

sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom)

optionals = Val((carbonates, oxygen, variable_redfield))
Expand Down Expand Up @@ -469,19 +473,19 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Any, <:Any, <:Val{(fals
LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, true)}, <:Any, <:Any}}


@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 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)

@inline carbon_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = nitrogen_flux(grid, advection, bgc, tracers, i, j) * bgc.organic_redfield
@inline carbon_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * bgc.organic_redfield

@inline nitrogen_flux(grid, advection, bgc::lobster_variable_redfield, tracers, i, j) =
sinking_flux(i, j, grid, advection, Val(:sPON), bgc, tracers) +
sinking_flux(i, j, grid, advection, Val(:bPON), bgc, tracers)
@inline nitrogen_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) =
sinking_flux(i, j, k, grid, advection, Val(:sPON), bgc, tracers) +
sinking_flux(i, j, k, grid, advection, Val(:bPON), bgc, tracers)

@inline carbon_flux(grid, advection, bgc::lobster_variable_redfield, tracers, i, j) =
sinking_flux(i, j, grid, advection, Val(:sPOC), bgc, tracers) +
sinking_flux(i, j, grid, advection, Val(:bPOC), bgc, tracers)
@inline carbon_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) =
sinking_flux(i, j, k, grid, advection, Val(:sPOC), bgc, tracers) +
sinking_flux(i, j, k, grid, advection, Val(:bPOC), bgc, tracers)

@inline remineralisation_receiver(::LOBSTER) = :NH₄
end # module
13 changes: 9 additions & 4 deletions src/Models/AdvectedPopulations/NPZD.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,10 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid,

sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom)

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"
end

return NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope,
base_maximum_growth,
nutrient_half_saturation,
Expand Down Expand Up @@ -339,9 +343,10 @@ adapt_structure(to, npzd::NPZD) =

adapt(to, npzd.particles))

@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 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 remineralisation_receiver(::NPZD) = :N
end # module
Loading

2 comments on commit ed7ec50

@jagoosw
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/90450

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.5.0 -m "<description of version>" ed7ec5006dc097efc81bc3ed75779f6bae3307da
git push origin v0.5.0

Please sign in to comment.