Sediment model fixes + simpler model

Manifest.toml

@@ -92,10 +92,10 @@ uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 version = "1.16.0"
 
 [[deps.ChangesOfVariables]]
-deps = ["LinearAlgebra", "Test"]
-git-tree-sha1 = "f84967c4497e0e1955f9a582c232b02847c5f589"
+deps = ["InverseFunctions", "LinearAlgebra", "Test"]
+git-tree-sha1 = "2fba81a302a7be671aefe194f0525ef231104e7f"
 uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
-version = "0.1.7"
+version = "0.1.8"
 
 [[deps.CommonDataModel]]
 deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf"]
@@ -258,9 +258,9 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.InverseFunctions]]
 deps = ["Test"]
-git-tree-sha1 = "6667aadd1cdee2c6cd068128b3d226ebc4fb0c67"
+git-tree-sha1 = "eabe3125edba5c9c10b60a160b1779a000dc8b29"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.9"
+version = "0.1.11"
 
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"

Project.toml

@@ -24,6 +24,7 @@ julia = "^1.8.0"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
@@ -34,4 +35,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 examples = ["CairoMakie", "NetCDF", "Interpolations", "Statistics", "Printf", "JLD2", "DataDeps"]
-test = ["Test", "CUDA", "DataDeps", "Statistics", "JLD2"]
+test = ["Test", "CUDA", "DataDeps", "Documenter", "Statistics", "JLD2"]

docs/make.jl

@@ -3,7 +3,7 @@ using Documenter, DocumenterCitations, Literate
 using OceanBioME
 using OceanBioME.SLatissimaModel: SLatissima
 using OceanBioME.LOBSTERModel: LOBSTER
-using OceanBioME.Boundaries.Sediments: SimpleMultiG
+using OceanBioME.Boundaries.Sediments: SimpleMultiG, InstantRemineralisation
 using OceanBioME.Boundaries: OCMIP_default, GasExchange
 
 using CairoMakie
@@ -53,7 +53,8 @@ model_parameters = (LOBSTER(; grid = BoxModelGrid()),
                     NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid()),
                     SLatissima(),
                     TwoBandPhotosyntheticallyActiveRadiation(; grid = BoxModelGrid()),
-                    SimpleMultiG(BoxModelGrid(); depth = 1000),
+                    SimpleMultiG(; grid = BoxModelGrid(), depth = 1000),
+                    InstantRemineralisation(; grid = BoxModelGrid()),
                     OCMIP_default,
                     GasExchange(; gas = :CO₂).condition.parameters,
                     GasExchange(; gas = :O₂).condition.parameters)

docs/src/model_components/sediment.md

@@ -4,13 +4,18 @@ We currently have one single-layer sediment model implemented. The model, propos
 
 It is straightforward to set up:
 
-```@example
-using OceanBioME, OceanBioME.Sediments
-using Oceananigans
+```jldoctest simplemultig; filter = r".*@ OceanBioME.Boundaries.Sediments.*"
+using OceanBioME, Oceananigans, OceanBioME.Sediments
 
 grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))
 
-sediment_model = SimpleMultiG(grid)
+sediment_model = SimpleMultiG(; grid)
+
+# output
+┌ Warning: Sediment models are an experimental feature and have not yet been validated.
+└ @ OceanBioME.Boundaries.Sediments ~/OceanBioME.jl/src/Boundaries/Sediments/simple_multi_G.jl:102
+[ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
+Single-layer multi-G sediment model (Float64)
 ```
 
 You may optionally specify the model parameters. This can then be passed in the setup of a BGC model:

src/Boundaries/Sediments/Sediments.jl

@@ -1,38 +1,40 @@
 module Sediments
 
-export SimpleMultiG
+export SimpleMultiG, InstantRemineralisation
 
 using KernelAbstractions
 using OceanBioME: ContinuousFormBiogeochemistry
 using Oceananigans
 using Oceananigans.Architectures: device
 using Oceananigans.Utils: launch!
-using Oceananigans.Advection: div_Uc
+using Oceananigans.Advection: advective_tracer_flux_z
 using Oceananigans.Units: day
 using Oceananigans.Fields: CenterField, Face
 using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity
+using Oceananigans.Grids: zspacing
+using Oceananigans.Operators: volume
+using Oceananigans.Fields: Center
 
-import Oceananigans.Biogeochemistry: update_tendencies!
 import Adapt: adapt_structure, adapt
 
 abstract type AbstractSediment end
 abstract type FlatSediment <: AbstractSediment end
 
 sediment_fields(::AbstractSediment) = ()
 
-@inline update_tendencies!(bgc::ContinuousFormBiogeochemistry{<:Any, <:FlatSediment, <:Any}, model) = update_tendencies!(bgc, bgc.sediment_model, model)
-@inline update_tendencies!(bgc, sediment, model) = nothing
+@inline update_sediment_tendencies!(bgc, sediment, model) = nothing
 
-function update_tendencies!(bgc, sediment::FlatSediment, model)
+function update_sediment_tendencies!(bgc, sediment::FlatSediment, model)
     arch = model.grid.architecture
 
-    for (i, tracer) in enumerate(sediment_tracers(sediment))    
-        launch!(arch, model.grid, :xy, store_flat_tendencies!, sediment.tendencies.Gⁿ[i], sediment.tendencies.G⁻[i])
+    for (i, tracer) in enumerate(sediment_tracers(sediment))
+        launch!(arch, model.grid, :xy, store_flat_tendencies!, sediment.tendencies.G⁻[i], sediment.tendencies.Gⁿ[i])
     end
 
     launch!(arch, model.grid, :xy,
             _calculate_tendencies!,
-            bgc.sediment_model, bgc, model.grid, model.tracers, model.timestepper)
+            bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper)
+
     return nothing
 end
 
@@ -41,7 +43,17 @@
     @inbounds G⁻[i, j, 1] = G⁰[i, j, 1]
 end
 
+
+@inline nitrogen_flux(grid, adveciton, bgc, tracers, i, j) = 0
+@inline carbon_flux(grid, adveciton, bgc, tracers, i, j) = 0
+@inline remineralisation_reciever(bgc, tendencies) = nothing
+
+@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())
+
 include("coupled_timesteppers.jl")
 include("simple_multi_G.jl")
+include("instant_remineralization.jl")
 
 end # module

src/Boundaries/Sediments/coupled_timesteppers.jl

@@ -1,9 +1,11 @@
-using Oceananigans: NonhydrostaticModel, prognostic_fields
+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!
 using Oceananigans.TurbulenceClosures: implicit_step!
+using Oceananigans.Models.HydrostaticFreeSurfaceModels: local_ab2_step!, ab2_step_free_surface!
+using Oceananigans.Architectures: AbstractArchitecture
 
 import Oceananigans.TimeSteppers: ab2_step!, rk3_substep!
 
@@ -44,6 +46,25 @@
     return nothing
 end
 
+@inline function ab2_step!(model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ)
+    # Step locally velocity and tracers
+    @apply_regionally local_ab2_step!(model, Δt, χ)
+
+    # blocking step for implicit free surface, non blocking for explicit
+    ab2_step_free_surface!(model.free_surface, model, Δt, χ)
+
+    sediment = model.biogeochemistry.sediment_model
+    arch = model.architecture
+
+    for (i, field) in enumerate(sediment_fields(sediment))
+        launch!(arch, model.grid, :xy, ab2_step_flat_field!, 
+                field, Δt, χ, 
+                sediment.tendencies.Gⁿ[i], 
+                sediment.tendencies.G⁻[i])
+
+    end
+end
+
 @kernel function ab2_step_flat_field!(u, Δt, χ, Gⁿ, G⁻)
     i, j = @index(Global, NTuple)
 
@@ -91,7 +112,6 @@
 
 @kernel function rk3_step_flat_field!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
     i, j = @index(Global, NTuple)
-
     @inbounds U[i, j, 1] += Δt * (γⁿ * Gⁿ[i, j, 1] + ζⁿ * G⁻[i, j, 1])
 end
 

src/Boundaries/Sediments/instant_remineralization.jl

@@ -0,0 +1,97 @@
+"""
+    struct InstantRemineralisation
+
+Hold the parameters and fields the simplest benthic boundary layer where
+organic carbon is assumed to remineralise instantly with some portion 
+becoming N, and a fraction being perminantly burried.
+
+Burial efficiency by 10.1029/2006GB002907, 2007
+"""
+struct InstantRemineralisation{FT, F, TE} <: FlatSediment
+          burial_efficiency_constant1 :: FT
+          burial_efficiency_constant2 :: FT
+    burial_efficiency_half_saturaiton :: FT
+
+                               fields :: F
+                           tendencies :: TE
+
+    function InstantRemineralisation(burial_efficiency_constant1::FT,
+                                     burial_efficiency_constant2::FT,
+                                     burial_efficiency_half_saturaiton::FT,
+                                     fields::F,
+                                     tendencies::TE) where {FT, F, TE}
+
+        return new{FT, F, TE}(burial_efficiency_constant1,
+                              burial_efficiency_constant2,
+                              burial_efficiency_half_saturaiton,
+                              fields,
+                              tendencies)
+    end
+end
+
+"""
+    InstantRemineralisation(; grid,
+        burial_efficiency_constant1::FT = 0.013,
+        burial_efficiency_constant2::FT = 0.53,
+        burial_efficiency_half_saturaiton::FT = 7) where FT
+
+Returns a single-layer instant remineralisaiton model for NPZD bgc models.
+
+Example
+=======
+
+```@example
+using OceanBioME, Oceananigans, OceanBioME.Sediments
+
+grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))
+
+sediment_model = InstantRemineralisation(; grid)
+```
+"""
+function InstantRemineralisation(; grid,
+            burial_efficiency_constant1 = 0.013,
+            burial_efficiency_constant2 = 0.53,
+            burial_efficiency_half_saturaiton = 7.0)
+
+    @warn "Sediment models are an experimental feature and have not yet been validated"
+
+    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)))
+
+    return InstantRemineralisation(burial_efficiency_constant1,
+                                   burial_efficiency_constant2,
+                                   burial_efficiency_half_saturaiton,
+                                   fields,
+                                   tendencies)
+end
+
+adapt_structure(to, sediment::InstantRemineralisation) = 
+    SimpleMultiG(sediment.burial_efficiency_constant1,
+                 sediment.burial_efficiency_constant2,
+                 sediment.burial_efficiency_half_saturaiton,
+                 adapt(to, sediment.fields),
+                 adapt(to, sediment.tendencies))
+
+sediment_tracers(::InstantRemineralisation) = (:N_storage, )
+sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, )
+
+@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, timestepper)
+    i, j = @index(Global, NTuple)
+
+    @inbounds begin                        
+        Δz = zspacing(i, j, 1, grid, Center(), Center(), Center())
+
+        flux = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz
+
+        burial_efficiency = sediment.burial_efficiency_constant1 + sediment.burial_efficiency_constant2 * ((flux / 6.56) / (7 + flux / 6.56)) ^ 2
+
+        # sediment evolution
+        sediment.tendencies.Gⁿ.N_storage[i, j, 1] = burial_efficiency * flux
+
+        remineralisation_reciever(bgc, timestepper.Gⁿ)[i, j, 1] += flux * (1 - burial_efficiency) / Δz
+    end
+end