From 132523c2e40c149da5090415660b490c704b45b2 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 14:19:21 +0000 Subject: [PATCH 01/11] bumped version --- Manifest.toml | 12 +++++++++--- Project.toml | 2 +- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 27e9b019d..55fd8bd81 100644 --- a/Manifest.toml +++ b/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"] diff --git a/Project.toml b/Project.toml index b70ab4ffc..a28314a58 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,7 @@ Adapt = "3" Documenter = "0.27" 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" From c392154a2cfa40622143c7aa728b0f8b3155ee28 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 14:19:33 +0000 Subject: [PATCH 02/11] fixed light attenuation --- src/Light/2band.jl | 8 ++++++-- src/Light/Light.jl | 8 +++++--- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/Light/2band.jl b/src/Light/2band.jl index 37562902b..0decc1c9a 100644 --- a/src/Light/2band.jl +++ b/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 diff --git a/src/Light/Light.jl b/src/Light/Light.jl index 2b36457da..c0d216c65 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -5,7 +5,7 @@ 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! @@ -16,8 +16,10 @@ 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 From 9471702253b81ebdf25a0978987d3e8517a7664a Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 14:21:37 +0000 Subject: [PATCH 03/11] fixed old `individuals` docs page --- .../src/model_components/individuals/index.md | 22 ++++++++++++------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index aebc8e8e4..ce1b897ad 100644 --- a/docs/src/model_components/individuals/index.md +++ b/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 From 3cad4fe5ee241928959aff4bd52492c4bdb7d577 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 14:29:08 +0000 Subject: [PATCH 04/11] fixed a load of surface PAR functions --- docs/src/model_implementation.md | 6 +++--- docs/src/quick_start.md | 4 +--- src/Light/Light.jl | 4 ++++ src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 6 +++--- src/Models/AdvectedPopulations/NPZD.jl | 4 ++-- 5 files changed, 13 insertions(+), 11 deletions(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index db204e208..1098d21f7 100644 --- a/docs/src/model_implementation.md +++ b/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 temp(z, t) = cos(t * 2π / year + 50days) + 28 -@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 diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index b27e43a7a..a4a21dc5a 100644 --- a/docs/src/quick_start.md +++ b/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) diff --git a/src/Light/Light.jl b/src/Light/Light.jl index c0d216c65..c69e35f22 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -31,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.0, cos(t * π / 12hours)) + end diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 422bc2b7e..8e6770b92 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/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, @@ -290,7 +290,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, diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 1134967c4..9fed9ab3f 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -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), @@ -168,7 +168,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), From 34819bfefe123629a9e012487847c754e250c450 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 14:30:13 +0000 Subject: [PATCH 05/11] forgot import --- src/Light/Light.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Light/Light.jl b/src/Light/Light.jl index c69e35f22..677982789 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -11,7 +11,7 @@ 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, From 0b26ce38ee1d85048da36a06fa7bcb1e369872ca Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 14:30:51 +0000 Subject: [PATCH 06/11] updated docs compats --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index c7439da61..7fc56f152 100644 --- a/docs/Project.toml +++ b/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" From 517b777807a19ac25b24b5a720c957236597181e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 7 Nov 2023 10:45:46 -0800 Subject: [PATCH 07/11] drop y-dependance since it's flat --- README.md | 6 +++--- docs/src/index.md | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index fbf77378e..c1886edc8 100644 --- a/README.md +++ b/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ᵢ) diff --git a/docs/src/index.md b/docs/src/index.md index b6b9e25de..b8ff8ff54 100644 --- a/docs/src/index.md +++ b/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ᵢ) From 389e2396d2c6d21b0d40bb9d30f8d1027d98a99c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 7 Nov 2023 10:48:37 -0800 Subject: [PATCH 08/11] more type stable --- src/Light/Light.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Light/Light.jl b/src/Light/Light.jl index 677982789..b9cb7a6ba 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -33,6 +33,6 @@ 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.0, cos(t * π / 12hours)) +default_surface_PAR(t) = 100 * max(0, cos(t * π / 12hours)) end From a2bb048b6814831530793c9af5be9ab3e2837cd0 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 7 Nov 2023 10:48:45 -0800 Subject: [PATCH 09/11] add missing import --- src/Models/AdvectedPopulations/NPZD.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 9fed9ab3f..3dc9ef56e 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/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 From b62b8195c1f7def4439bcee0aa8cfd24ec218a90 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 7 Nov 2023 21:21:41 +0000 Subject: [PATCH 10/11] fix winking speed function --- docs/src/model_implementation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 1098d21f7..0dd37bc50 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -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) From 476ad87224a5ae8d26ab61ea76445fbe216c4936 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 8 Nov 2023 12:22:25 +0000 Subject: [PATCH 11/11] removed `Relaxation` since it temporarily doesn't work --- docs/src/model_implementation.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 0dd37bc50..11c5428c9 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -235,7 +235,7 @@ using OceanBioME.Sediments: InstantRemineralisation @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(z, t) = cos(t * 2π / year + 50days) + 28 +@inline ∂ₜT(z, t) = - 2π / year * sin(t * 2π / year + 50days) @inline κₜ(z, t) = 1e-2 * (1 + tanh((z - 50) / 10)) / 2 + 1e-4 @@ -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)