Merge pull request #157 from OceanBioME/jsw/bump-oceananigans-0.90

Bump Oceananigans to 0.90

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.3"
 manifest_format = "2.0"
-project_hash = "abd05af8b9fe6476515dbe2a1c676a74ac1e6059"
+project_hash = "9c3a43fcc42b0e8b4f92f9c38f0c935f93e8165c"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -501,9 +501,15 @@ version = "1.2.0"
 
 [[deps.Oceananigans]]
 deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
-git-tree-sha1 = "e8caa94d03ec34434abc0fb036d43d34d1903411"
+git-tree-sha1 = "7ef7083f3fb79f225c50716bee81265c626ff98d"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.89.2"
+version = "0.90.1"
+
+    [deps.Oceananigans.extensions]
+    OceananigansEnzymeCoreExt = "EnzymeCore"
+
+    [deps.Oceananigans.weakdeps]
+    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]

Project.toml

@@ -16,7 +16,7 @@ Adapt = "3"
 Documenter = "1"
 JLD2 = "0.4"
 KernelAbstractions = "0.9"
-Oceananigans = "0.84.1, 0.85 - 0.89"
+Oceananigans = "0.90"
 Roots = "2"
 StructArrays = "0.4, 0.5, 0.6"
 julia = "1.9"

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/Project.toml

@@ -18,4 +18,4 @@ CairoMakie = "0.10"
 Documenter = "1"
 DocumenterCitations = "1"
 Literate = "≥2.9.0"
-Oceananigans = "^0.89.2"
+Oceananigans = "0.90"

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/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
 
@@ -251,7 +251,7 @@ sediment = InstantRemineralisation(; grid)
 
 sinking_velocity = ZFaceField(grid)
 
-w_sink(x, y, z) = 2 / day * tanh(z / 5)
+w_sink(z) = 2 / day * tanh(z / 5)
 
 set!(sinking_velocity, w_sink)
 
@@ -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

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
               disolved_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,
@@ -284,7 +284,7 @@ function LOBSTER(; grid,
                    disolved_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, 

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