Adds IronPhosphate sediment model

Project.toml

@@ -6,6 +6,9 @@ version = "0.10.3"
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+InteractiveErrors = "6cff3a6f-7015-4b3b-b269-4b312a73b9bf"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"

src/Boundaries/Sediments/Sediments.jl

@@ -1,6 +1,6 @@
 module Sediments
 
-export SimpleMultiG, InstantRemineralisation
+export SimpleMultiG, InstantRemineralisation, IronPhosphate
 
 using KernelAbstractions
 
@@ -62,6 +62,9 @@ end
 
 @inline nitrogen_flux() = 0
 @inline carbon_flux() = 0
+@inline iron_flux() = 0
+@inline oxygen_flux() = 0
+@inline poc_flux() = 0
 @inline remineralisation_receiver() = nothing
 @inline sinking_tracers() = nothing
 @inline required_tracers() = nothing
@@ -106,5 +109,6 @@ end
 include("coupled_timesteppers.jl")
 include("simple_multi_G.jl")
 include("instant_remineralization.jl")
+include("IronPhosphate.jl")
 
 end # module

src/LDNtesting.jl

@@ -0,0 +1,137 @@
+using OceanBioME, Test, CUDA, Oceananigans, JLD2, Documenter
+
+using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation
+using Oceananigans.Units
+
+using OceanBioME.Sediments: sediment_tracers, sediment_fields
+using Oceananigans: Field
+using Oceananigans.Fields: TracerFields
+
+using Oceananigans.Operators: volume, Azᶠᶜᶜ
+
+using OceanBioME.LOBSTERModel: VariableRedfieldLobster
+
+
+
+architecture = CUDA.has_cuda() ? GPU() : CPU()
+
+function intercept_tracer_tendencies!(model, intercepted_tendencies)
+    for (name, field) in enumerate(intercepted_tendencies)
+        field .= Array(interior(model.timestepper.Gⁿ[name + 3]))
+    end
+end
+
+function set_defaults!(sediment::SimpleMultiG)
+    set!(sediment.fields.N_fast, 0.0230)
+    set!(sediment.fields.N_slow, 0.0807)
+
+    set!(sediment.fields.C_fast, 0.5893)
+    set!(sediment.fields.C_slow, 0.1677)
+end 
+
+
+
+function set_defaults!(::VariableRedfieldLobster, model)
+
+    set!(model, Z = 0.5363, 
+                NO₃ = 2.3103, NH₄ = 0.0010, 
+                DIC = 2106.9, Alk = 2408.9, 
+                O₂ = 258.92, 
+                DOC = 5.3390, DON = 0.8115,
+                sPON = 0.2299, sPOC = 1.5080,
+                bPON = 0.0103, bPOC = 0.0781)
+
+
+
+    kick = 0.05
+    uᵢ(x, y, z) = kick * randn()
+    vᵢ(x, y, z) = kick * randn()
+    wᵢ(x, y, z) = kick * randn()
+    bᵢ(x, y, z) = kick * randn() + 1
+    Pᵢ(x, y, z) = (1000-z)/1500
+    #Pᵢ(x, y, z) = 0.4686 + exp(-((z - 500) / 50)^2)
+
+    set!(model, u = uᵢ, v = vᵢ, w = wᵢ, b = bᵢ, P = Pᵢ)
+    @info "flag"
+
+end
+
+
+total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) + 
+                                    sum(sed.fields.N_slow) + 
+                                    sum(sed.fields.N_ref)
+
+
+total_nitrogen(::VariableRedfieldLobster, model) = sum(model.tracers.NO₃) +
+                                                     sum(model.tracers.NH₄) +
+                                                     sum(model.tracers.P) +
+                                                     sum(model.tracers.Z) +
+                                                     sum(model.tracers.DON) +
+                                                     sum(model.tracers.sPON) +
+                                                     sum(model.tracers.bPON)
+
+
+
+function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAdamsBashforth2)
+    Re = 5000
+    model = model(; grid, 
+                    biogeochemistry,
+                    closure = nothing,
+                    buoyancy = Buoyancy(model=BuoyancyTracer()),
+                    tracers = :b)
+
+                    @info "flag1"
+    set_defaults!(model.biogeochemistry.sediment)
+
+    set_defaults!(biogeochemistry.underlying_biogeochemistry, model)
+
+    simulation = Simulation(model, Δt = 50, stop_time = 1day)
+                    @info "flag2"
+    intercepted_tendencies = Tuple(Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid)))
+
+    simulation.callbacks[:intercept_tendencies] = Callback(intercept_tracer_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies)
+
+    simulation.output_writers[:tracers] = JLD2OutputWriter(model, model.tracers,
+                                                       filename = "LDNtesting.jld2",
+                                                       schedule = TimeInterval(24minute),
+                                                       overwrite_existing = true)
+@info "flag3"
+
+    run!(simulation)
+    return nothing
+end
+
+bottom_height(x, y) = -1500 + 1000 * exp(- (x^2 + y^2) / 250) # a perfect hill
+
+display_name(::LOBSTER) = "LOBSTER"
+display_name(::NutrientPhytoplanktonZooplanktonDetritus) = "NPZD"
+display_name(::SimpleMultiG) = "Multi-G"
+display_name(::InstantRemineralisation) = "Instant remineralisation"
+display_name(::RectilinearGrid) = "Rectilinear grid"
+display_name(::LatitudeLongitudeGrid) = "Latitude longitude grid"
+display_name(::ImmersedBoundaryGrid) = "Immersed boundary grid"
+
+
+    grid = #RectilinearGrid(architecture; size=(3, 3, 50), extent=(10, 10, 500))
+    ImmersedBoundaryGrid(
+        LatitudeLongitudeGrid(architecture; size = (10, 10, 50), latitude = (-10, 10), longitude = (-10, 10), z = (-1000, 0)), 
+        GridFittedBottom(bottom_height))
+
+    timestepper = :QuasiAdamsBashforth2
+    sediment_model = SimpleMultiG(; grid)
+    model =  HydrostaticFreeSurfaceModel
+    biogeochemistry = LOBSTER(; grid,
+                                carbonates = ifelse(isa(sediment_model, SimpleMultiG), true, false), 
+                                oxygen = ifelse(isa(sediment_model, SimpleMultiG), true, false), 
+                                variable_redfield = ifelse(isa(sediment_model, SimpleMultiG), true, false), 
+                                sediment_model)
+
+    @info "Running sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))"
+    test_flat_sediment(grid, biogeochemistry, model; timestepper)                              
+
+
+
+
+
+
+

src/LDNtesting2.jl

@@ -0,0 +1,43 @@
+T = FieldTimeSeries("buoyancy_front.jld2", "T")
+N = FieldTimeSeries("buoyancy_front.jld2", "N")
+P = FieldTimeSeries("buoyancy_front.jld2", "P")
+
+xc, yc, zc = nodes(T)
+
+times = T.times
+
+using CairoMakie
+
+n = Observable(1)
+
+T_lims = (8.94, 9.06)
+N_lims = (0, 4.5)
+P_lims = (0.007, 0.02)
+
+Tₙ = @lift interior(T[$n], :, 1, :)
+Nₙ = @lift interior(N[$n], :, 1, :)
+Pₙ = @lift interior(P[$n], :, 1, :)
+
+fig = Figure(size = (1000, 520), fontsize = 20)
+
+title = @lift "t = $(prettytime(times[$n]))"
+Label(fig[0, :], title)
+
+axis_kwargs = (xlabel = "x (m)", ylabel = "z (m)", width = 770, yticks = [-400, -200, 0])
+ax1 = Axis(fig[1, 1]; title = "Temperature (°C)", axis_kwargs...)
+ax2 = Axis(fig[2, 1]; title = "Nutrients concentration (mmol N / m³)",axis_kwargs...)
+ax3 = Axis(fig[3, 1]; title = "Phytoplankton concentration (mmol N / m³)", axis_kwargs...)
+
+hm1 = heatmap!(ax1, xc, zc, Tₙ, colorrange = T_lims, colormap = Reverse(:lajolla), interpolate = true)
+hm2 = heatmap!(ax2, xc, zc, Nₙ, colorrange = N_lims, colormap = Reverse(:bamako), interpolate = true)
+hm3 = heatmap!(ax3, xc, zc, Pₙ, colorrange = P_lims, colormap = Reverse(:bamako), interpolate = true)
+
+Colorbar(fig[1, 2], hm1, ticks = [8.95, 9.0, 9.05])
+Colorbar(fig[2, 2], hm2, ticks = [0, 2, 4])
+Colorbar(fig[3, 2], hm3, ticks = [0.01, 0.02, 0.03])
+
+rowgap!(fig.layout, 0)
+
+record(fig, "buoyancy_front.gif", 1:length(times)) do i
+    n[] = i
+end

src/LDNtesting4.jl

@@ -0,0 +1,61 @@
+using OceanBioME, Test, CUDA, Oceananigans, JLD2, Documenter
+
+using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation
+using Oceananigans.Units
+
+using OceanBioME.Sediments: sediment_tracers, sediment_fields
+using Oceananigans: Field
+using Oceananigans.Fields: TracerFields
+
+using Oceananigans.Operators: volume, Azᶠᶜᶜ
+
+using OceanBioME.LOBSTERModel: VariableRedfieldLobster
+
+using CairoMakie
+using CairoMakie: record
+
+Z = FieldTimeSeries("LDNtesting.jld2", "b")
+N = FieldTimeSeries("LDNtesting.jld2", "Z")
+P = FieldTimeSeries("LDNtesting.jld2", "P")
+
+xc, yc, zc = nodes(Z)
+
+times = Z.times
+
+print(nodes(Z))
+9
+Z_lims = (findmin(Z.data)[1], findmax(Z.data)[1])
+N_lims = (findmin(N.data)[1], findmax(N.data)[1])
+P_lims = (findmin(P.data)[1], findmax(P.data)[1])
+
+n = Observable(1)
+
+Zₙ = @lift interior(Z[$n], :, 2, :)
+Nₙ = @lift interior(N[$n], :, 2, :)
+Pₙ = @lift interior(P[$n], :, 2, :)
+
+println(Pₙ)
+
+fig = Figure(size = (1000, 520), fontsize = 20)
+
+title = @lift "t = $(prettytime(times[$n]))"
+Label(fig[2, 3], title)
+
+axis_kwargs = (xlabel = "x (m)", ylabel = "z (m)", width = 150) #, yticks = [-400, -200, 0])
+ax1 = Axis(fig[1, 1]; title = "Zooplankton concentration\n (mmol N / m³)", axis_kwargs...)
+ax2 = Axis(fig[1, 3]; title = "NO₃ concentration\n (mmol N / m³)",axis_kwargs...)
+ax3 = Axis(fig[1, 5]; title = "Phytoplankton concentration\n (mmol N / m³)", axis_kwargs...)
+
+hm1 = heatmap!(ax1, xc, zc, Zₙ, colorrange = Z_lims, colormap = Reverse(:bamako), interpolate = true)
+hm2 = heatmap!(ax2, xc, zc, Nₙ, colorrange = N_lims, colormap = Reverse(:lajolla), interpolate = true)
+hm3 = heatmap!(ax3, xc, zc, Pₙ, colorrange = P_lims, colormap = Reverse(:bamako), interpolate = true)
+
+Colorbar(fig[1, 2], hm1)
+Colorbar(fig[1, 4], hm2)
+Colorbar(fig[1, 6], hm3)
+
+rowgap!(fig.layout, 0)
+
+record(fig, "LDNtesting.gif", 1:length(times)) do i
+    n[] = i
+end 

src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl

@@ -64,7 +64,7 @@ import OceanBioME: maximum_sinking_velocity
 import Adapt: adapt_structure, adapt
 import Base: show, summary
 
-import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers
+import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, iron_flux, oxygen_flux, poc_flux, remineralisation_receiver, sinking_tracers
 
 struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry
     phytoplankton_preference :: FT
@@ -486,6 +486,7 @@ const VariableRedfieldLobster = Union{LOBSTER{<:Any, <:Val{(false, false, true)}
     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₄
 
 @inline conserved_tracers(::LOBSTER) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)