Merge pull request #220 from OceanBioME/jsw/giant-kelp-growth

Make particles easier to use

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "90c0672acfd9ac1461839f7013be62ced327f62a"
+project_hash = "dcef9b8c3a1cc421800a361070f72b558c4aafe7"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]

Project.toml

@@ -1,15 +1,17 @@
 name = "OceanBioME"
 uuid = "a49af516-9db8-4be4-be45-1dad61c5a376"
 authors = ["Jago Strong-Wright <[email protected]> and contributors"]
-version = "0.12.0"
+version = "0.13.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 GibbsSeaWater = "9a22fb26-0b63-4589-b28e-8f9d0b5c3d05"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
 

docs/make.jl

@@ -1,7 +1,7 @@
 using Documenter, DocumenterCitations, Literate
 
 using OceanBioME
-using OceanBioME: SLatissima, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus
+using OceanBioME: SugarKelp, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus
 using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation
 using OceanBioME: CarbonChemistry, GasExchange
 
@@ -53,7 +53,7 @@ if !isdir(OUTPUT_DIR) mkdir(OUTPUT_DIR) end
 
 model_parameters = (LOBSTER(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry,
                     NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid(), light_attenuation_model = nothing).underlying_biogeochemistry,
-                    SLatissima(),
+                    SugarKelp(),
                     TwoBandPhotosyntheticallyActiveRadiation(; grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))),
                     SimpleMultiG(; grid = BoxModelGrid()),
                     InstantRemineralisation(; grid = BoxModelGrid()),

docs/src/appendix/library.md

@@ -66,7 +66,7 @@ Modules = [OceanBioME.Models.PISCESModel.Nitrogen]
 ### Sugar kelp (Saccharina latissima)
 
 ```@autodocs
-Modules = [OceanBioME.Models.SLatissimaModel]
+Modules = [OceanBioME.Models.SugarKelpModel]
 ```
 
 ### Carbon Chemistry 
@@ -98,3 +98,9 @@ Modules = [OceanBioME.Models.GasExchangeModel, OceanBioME.Models.GasExchangeMode
 ```@autodocs
 Modules = [OceanBioME.BoxModels]
 ```
+
+## Particles
+
+```@autodocs
+Modules = [OceanBioME.Particles]
+```

docs/src/index.md

@@ -4,7 +4,7 @@ OceanBioME.jl is a fast and flexible ocean biogeochemical modelling environment.
 
 OceanBioME.jl currently provides a core of several biogeochemical models Nutrient--Phytoplankton--Zooplankton--Detritus ([NPZD](@ref NPZD)), [LOBSTER](https://doi.org/10.1029/2004JC002588), a medium complexity model, and an early implementation of [PISCES](https://www.pisces-community.org/), a complex model. It also provides essential utilities like air-sea gas exchange models to provide appropriate top boundary conditions, a carbon chemistry model for computing the pCO₂, and sediment models to for the benthic boundary. 
 
-OceanBioME.jl includes a framework for integrating the growth of biological/active Lagrangian particles which move around and can interact with the (Eulerian) tracer fields - for example, consuming nutrients and carbon dioxide while releasing dissolved organic material. A growth model for sugar kelp is currently implemented using active particles, and this model can be used in a variety of dynamical scenarios including free-floating or bottom-attached particles.
+OceanBioME.jl includes a framework for integrating the growth of [biological/active particles](@ref individuals) which move around and can interact with the (Eulerian) tracer fields - for example, consuming nutrients and carbon dioxide while releasing dissolved organic material. A growth model for sugar kelp is currently implemented using active particles, and this model can be used in a variety of dynamical scenarios including free-floating or bottom-attached particles.
 
 ## Quick install
 

docs/src/model_components/biogeochemical/LOBSTER.md

@@ -95,15 +95,15 @@ When the carbonate chemistry is activated additional tracers ``DIC`` and ``Alk``
 
 ### Oxygen chemistry
 
-When the oxygen chemistry is activated additional tracer ``O_2`` evolve like:
+When the oxygen chemistry is activated, additional tracer ``O_2`` evolve like:
 
 ```math
 \frac{\partial O_2}{\partial t} = \mu_P L_{PAR}\left(L_{NO_3} + L_{NH_4}\right)R_{O_2}P - (R_{O_2} - R_{nit})\frac{\partial NH_4}{\partial t} - R_{O_2}\mu_nNH_4.
 ```
 
 ### Variable Redfield
 
-When the variable Redfield modification is activated the organic components are modified to evolve their nitrogen and carbon content separately. This means that the waste from non-Redfield models (e.g. loss from the [kelp](@ref SLatissima)) can be accounted for.
+When the variable Redfield modification is activated the organic components are modified to evolve their nitrogen and carbon content separately. This means that the waste from non-Redfield models (e.g. loss from the [kelp](@ref sugar-kelp)) can be accounted for.
 
 In this case the organic components are split into nitrogen and carbon compartments, so the tracers ``sPOM``, ``bPOM``, and ``DOM`` are replaced with ``sPON``, ``sPOC``, ``bPON``, ``bPOC``, ``DON``, and ``DOC``. The nitrogen compartments evolve as per the organic matter equations above (i.e. replacing each ``XOM`` with ``XON``), while the carbon compartments evolve like:
 

docs/src/model_components/individuals/index.md

@@ -1,92 +1,50 @@
 # [Individuals](@id individuals)
 
-The effects of individuals can be modelled in OceanBioME. We have implemented this through custom dynamics in the [Lagrangian Particle tracking feature of Oceananigans](https://clima.github.io/OceananigansDocumentation/stable/model_setup/lagrangian_particles/). We have extended these functionalities to make it easier to implement "active" particles which interact with the tracers. We have then implemented a model of [sugar kelp](@ref SLatissima) which can be followed as an example of using this functionality.
+The effects of individuals can be modelled in OceanBioME. We have implemented this through custom dynamics in the [Lagrangian Particle tracking feature of Oceananigans](https://clima.github.io/OceananigansDocumentation/stable/model_setup/lagrangian_particles/). We have extended these functionalities to make it easier to implement "active" particles which interact with the tracers. We have then implemented a model of [sugar kelp](@ref sugar-kelp) which can be followed as an example of using this functionality.
 
-To setup particles first create a particle type with the desired properties, e.g.:
+To setup particles first create a particle biogeochemistry, e.g.:
 
 ```@example particles
-using OceanBioME.Particles: BiogeochemicalParticles
 
-struct GrowingParticles{FT, VT} <: BiogeochemicalParticles 
+struct GrowingParticles{FT}
     nutrients_half_saturation :: FT
-
-    size :: VT
-    nitrate_uptake :: VT
-
-    x :: VT
-    y :: VT
-    z :: VT
 end
 ```
 
-You then need to overload particular functions to integrate the growth, so they need to first be `import`ed:
-
-```@example particles
-import Oceananigans.Biogeochemistry: update_tendencies!
-import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!
-```
-
-First, to integrate the particles properties we overload `update_lagrangian_particle_properties!`;
-in this fictitious case we will have a Mondo-quota nutrient uptake and growth:
+We then need to add some methods to tell `OceanBioME` what properties this particle has, and what tracers it interacts with:
 
 ```@example particles
-using Oceananigans.Fields: interpolate
-
-function update_lagrangian_particle_properties!(particles::GrowingParticles, model, bgc, Δt)
-    @inbounds for p in 1:length(particles)
-        nutrients = @inbounds interpolate(model.tracers.NO₃, particle.x[p], particle.y[p], particle.z[p])
 
-        uptake = nutrients / (particle.nutrients_half_saturation + nutrients)
+import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers
 
-        particles.size[p] += uptake * Δt
-        particles.nitrate_uptake[p] = uptake
-    end
-    return nothing
-end
+required_particle_fields(::GrowingParticles) = (:S, )
+required_tracers(::GrowingParticles) = (:N, )
+coupled_tracers(::GrowingParticles) = (:N, )
 
-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:
-
+So our model is going to track the `S`ize of the particles and take up `N`utrients. 
+Now we need to how this growth happens. 
+The forcing functions should be of the form `(particles::ParticleBiogeochemistry)(::Val{:PROPERTY}, t, required_particle_fields..., required_tracers...)`, so in this example:
 ```@example particles
-using OceanBioME.Particles: get_node
-
-function update_tendencies!(bgc, particles::GrowingParticles, model)
-    @inbounds for p in 1:length(particles)
-        i, j, k = fractional_indices((x, y, z), grid, Center(), Center(), Center())
-
-        # 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)
+(p::GrowingParticles)(::Val{:S}, t, S, N) = N / (N + p.nutrient_half_saturation)
+(p::GrowingParticles)(::Val{:N}, t, S, N) = - N / (N + p.nutrient_half_saturation)
+```
 
-        # 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))
+We can then create an instance of this particle model using `BiogeochemicalParticles`, and set their initial position and size:
+```@example particles
+using OceanBioME, Oceananigans
 
-        node_volume = volume(i, j, k, grid, Center(), Center(), Center())
+Lx, Ly, Lz = 100, 100, 100
+grid = RectilinearGrid(; size = (8, 8, 8), extent = (Lx, Ly, Lz))
 
-        model.timestepper.Gⁿ.NO₃[i, j, k] += particles.nitrate_uptake[p] / (d * node_volume)
-    end
-    return nothing
-end
+particles = BiogeochemicalParticles(10; grid, biogeochemistry = GrowingParticles(0.5))
 
-nothing #hide
+set!(particles, S = 0.1, x = rand(10) * Lx, y = rand(10) * Ly, z = rand(10) * Lz)
 ```
 
-Now we can just plug this into any biogeochemical model setup to have particles (currently [NPZD](@ref NPZD) and [LOBSTER](@ref LOBSTER)):
-
+We can then put these into a compatible biogeochemical model, for example:
 ```@example particles
-using OceanBioME, Oceananigans
-
-Lx, Ly, Lz = 1000, 1000, 100
-grid = RectilinearGrid(; size = (64, 64, 16), extent = (Lx, Ly, Lz))
-
-# Start the particles randomly distributed, floating on the surface
-particles = GrowingParticles(0.5, zeros(3), zeros(3), rand(3) * Lx, rand(3) * Ly, zeros(3))
 
-biogeochemistry = LOBSTER(; grid, particles)
+biogeochemistry = NPZD(; grid, particles)
 ```

docs/src/model_components/individuals/slatissima.md

@@ -1,4 +1,4 @@
-# [Sugar kelp (Saccharina latissima) individuals](@id SLatissima)
+# [Sugar kelp (Saccharina latissima) individuals](@id sugar-kelp)
 
 We have implemented a model of sugar kelp growth within this spatially infinitesimal Lagrangian particles framework originally based on the model of [Broch2012](@citet) and updated by [Broch2013](@citet), [Fossberg2018](@citet), and [Broch2019](@citet). This is the same model passively forced by [StrongWright2022](@citet).
 
@@ -7,6 +7,29 @@ The model tracks three variables, the frond area, A (dm²), carbon reserve, C (g
 Results could look something like this (from [StrongWright2022](@citet)):
 ![Example A, N, and C profiles from [StrongWright2022](@citet)](https://www.frontiersin.org/files/Articles/793977/fmars-08-793977-HTML/image_m/fmars-08-793977-g002.jpg)
 
+You can access the model biogeochemistry by setting up `SugarKelp`, i.e.:
+```jldoctest
+using OceanBioME
+
+kelp_bgc = SugarKelp()
+
+# output
+SugarKelp{FT} biogeochemistry (Broch & Slagstad, 2012) tracking the `N`itrogen and `C`arbon in a frond of `A`rea
+```
+which can be put into `BiogeochemicalParticles`, or you can directly manifest particles:
+```jldoctest
+using OceanBioME, Oceananigans
+
+grid = RectilinearGrid(size = (1, 1, 1), extent = (1, 1, 1));
+particles = SugarKelpParticles(10; grid)
+
+# output
+10 BiogeochemicalParticles with SugarKelp{FT} biogeochemistry:
+├── fields: (:A, :N, :C)
+└── coupled tracers: (:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON)
+
+```
+
 ## Model equations
 
 As per [Broch2012](@citet) this model variables evolve as: