From 658302791bac9da5600ff040dd6a75e8007e7dfb Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 17:12:56 +0100 Subject: [PATCH 001/112] fixed gpu incompatability in `calculate_bottom_indices` --- src/Boundaries/Sediments/Sediments.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index e295ab63c..8248f6228 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -5,7 +5,7 @@ export SimpleMultiG, InstantRemineralisation using KernelAbstractions using OceanBioME: ContinuousFormBiogeochemistry, BoxModelGrid using Oceananigans -using Oceananigans.Architectures: device +using Oceananigans.Architectures: device, architecture, arch_array using Oceananigans.Utils: launch! using Oceananigans.Advection: advective_tracer_flux_z using Oceananigans.Units: day @@ -76,7 +76,8 @@ calculate_bottom_indices(grid) = ones(Int, size(grid)[1:2]...) end function calculate_bottom_indices(grid::ImmersedBoundaryGrid) - indices = zeros(Int, size(grid)[1:2]...) + arch = architecture(grid) + indices = arch_array(arch, zeros(Int, size(grid)[1:2]...)) launch!(grid.architecture, grid, :xy, find_bottom_cell, grid, indices) From 20524a9840fddf6bc54013edff47d95eb898c3e9 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 17:18:49 +0100 Subject: [PATCH 002/112] fixed adapt method for instant remineralisation (oops) --- src/Boundaries/Sediments/instant_remineralization.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 1ea6e8d09..c5d969925 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -74,11 +74,11 @@ function InstantRemineralisation(; grid, 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)) + InstantRemineralisation(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, ) From cebc1fac32a9d231610214812db45081234aefea Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 17:22:21 +0100 Subject: [PATCH 003/112] and again --- src/Boundaries/Sediments/instant_remineralization.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index c5d969925..4b433e649 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -78,7 +78,8 @@ adapt_structure(to, sediment::InstantRemineralisation) = sediment.burial_efficiency_constant2, sediment.burial_efficiency_half_saturaiton, adapt(to, sediment.fields), - adapt(to, sediment.tendencies)) + adapt(to, sediment.tendencies), + adapt(to, sediment.bottom_indices)) sediment_tracers(::InstantRemineralisation) = (:N_storage, ) sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, ) From 13b040524c6736bc321671d0afe9290caa471da7 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 21:14:04 +0100 Subject: [PATCH 004/112] maybe fixed sediment --- src/Boundaries/Sediments/Sediments.jl | 2 +- src/Boundaries/Sediments/instant_remineralization.jl | 2 +- src/Boundaries/Sediments/simple_multi_G.jl | 10 +++++----- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 8248f6228..675ffa9f2 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -36,7 +36,7 @@ function update_sediment_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper) + bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ) return nothing end diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 4b433e649..8503a1f87 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -86,7 +86,7 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st @inline bottom_index_array(sediment::InstantRemineralisation) = sediment.bottom_indices -@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, timestepper) +@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, tendencies) i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 77e804910..49b2c0367 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -156,7 +156,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, @inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices -@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, timestepper) +@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, tendencies) i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) @@ -226,10 +226,10 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, Δz = grid.Δzᵃᵃᶜ[1] - timestepper.Gⁿ.NH₄[i, j, 1] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz - timestepper.Gⁿ.NO₃[i, j, 1] += Nᵐⁱⁿ * pₙᵢₜ / Δz - timestepper.Gⁿ.DIC[i, j, 1] += Cᵐⁱⁿ / Δz - timestepper.Gⁿ.O₂[i, j, 1] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly + tendencies.Gⁿ.NH₄[i, j, 1] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz + tendencies.Gⁿ.NO₃[i, j, 1] += Nᵐⁱⁿ * pₙᵢₜ / Δz + tendencies.Gⁿ.DIC[i, j, 1] += Cᵐⁱⁿ / Δz + tendencies.Gⁿ.O₂[i, j, 1] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly end end From 29514bf907f38df956d80bb8a518823297f7ca32 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 7 Sep 2023 21:19:49 +0100 Subject: [PATCH 005/112] maybe fixed sediment again --- src/Boundaries/Sediments/instant_remineralization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 8503a1f87..f25c5b614 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -100,5 +100,5 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st # sediment evolution @inbounds sediment.tendencies.Gⁿ.N_storage[i, j, 1] = burial_efficiency * flux - @inbounds timestepper.Gⁿ[remineralisation_receiver(bgc)][i, j, 1] += flux * (1 - burial_efficiency) / Δz + @inbounds tendencies[remineralisation_receiver(bgc)][i, j, 1] += flux * (1 - burial_efficiency) / Δz end From d9a210aa1a5237c4920b22bb6d97abdc1c9f443c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 8 Sep 2023 12:42:34 +0100 Subject: [PATCH 006/112] maybe fixed --- src/Boundaries/Sediments/Sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 7414fbb36..c322930e3 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -37,7 +37,7 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ) + bgc.sediment, bgc, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ) return nothing end From a67c79e0b0572d77b9149903fd0de701855cbe0f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 8 Sep 2023 14:53:33 +0100 Subject: [PATCH 007/112] fixed multi g --- src/Boundaries/Sediments/simple_multi_G.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 49b2c0367..f0919ab3d 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -226,10 +226,10 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, Δz = grid.Δzᵃᵃᶜ[1] - tendencies.Gⁿ.NH₄[i, j, 1] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz - tendencies.Gⁿ.NO₃[i, j, 1] += Nᵐⁱⁿ * pₙᵢₜ / Δz - tendencies.Gⁿ.DIC[i, j, 1] += Cᵐⁱⁿ / Δz - tendencies.Gⁿ.O₂[i, j, 1] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly + tendencies.NH₄[i, j, 1] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz + tendencies.NO₃[i, j, 1] += Nᵐⁱⁿ * pₙᵢₜ / Δz + tendencies.DIC[i, j, 1] += Cᵐⁱⁿ / Δz + tendencies.O₂[i, j, 1] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly end end From 125cf0a645d51017a48f697c2ca054dc51f71a2e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 11 Sep 2023 15:50:16 +0100 Subject: [PATCH 008/112] added GPU detection to tests --- test/runtests.jl | 4 +++- test/test_LOBSTER.jl | 15 +++++++------- test/test_NPZD.jl | 15 +++++++------- test/test_light.jl | 6 ++---- test/test_sediments.jl | 46 ++++++++++++++++++++--------------------- test/test_slatissima.jl | 3 +-- test/test_utils.jl | 8 +++---- 7 files changed, 45 insertions(+), 52 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 8dff45c8d..69a36e363 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,6 @@ -using OceanBioME, Documenter, Test +using OceanBioME, Documenter, Test, CUDA + +architecture = CUDA.has_cuda() ? GPU() : CPU() include("test_utils.jl") include("test_light.jl") diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index db25b6d49..c7390d616 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -114,14 +114,13 @@ end n_timesteps = 100 -for arch in (CPU(), ) - grid = RectilinearGrid(arch; size=(1, 1, 1), extent=(1, 1, 2)) - for open_bottom = (false, true), sinking = (false, true), variable_redfield = (false, true), oxygen = (false, true), carbonates = (false, true) - if !(sinking && open_bottom) # no sinking is the same with and without open bottom - @info "Testing on $(typeof(arch)) with carbonates $(carbonates ? :✅ : :❌), oxygen $(oxygen ? :✅ : :❌), variable redfield $(variable_redfield ? :✅ : :❌), sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" - @testset "$arch, $carbonates, $oxygen, $variable_redfield, $sinking, $open_bottom" begin - test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open_bottom, n_timesteps) - end +grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 2)) + +for open_bottom = (false, true), sinking = (false, true), variable_redfield = (false, true), oxygen = (false, true), carbonates = (false, true) + if !(sinking && open_bottom) # no sinking is the same with and without open bottom + @info "Testing on $(typeof(arch)) with carbonates $(carbonates ? :✅ : :❌), oxygen $(oxygen ? :✅ : :❌), variable redfield $(variable_redfield ? :✅ : :❌), sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" + @testset "$arch, $carbonates, $oxygen, $variable_redfield, $sinking, $open_bottom" begin + test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open_bottom, n_timesteps) end end end \ No newline at end of file diff --git a/test/test_NPZD.jl b/test/test_NPZD.jl index 465d41266..945070fb1 100644 --- a/test/test_NPZD.jl +++ b/test/test_NPZD.jl @@ -47,14 +47,13 @@ function test_NPZD(grid, sinking, open_bottom) return nothing end -for arch in (CPU(), ) - grid = RectilinearGrid(arch; size=(3, 3, 6), extent=(1, 1, 2)) - for sinking = (false, true), open_bottom = (false, true) - if !(sinking && open_bottom) # no sinking is the same with and without open bottom - @info "Testing on $(typeof(arch)) with sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" - @testset "$arch, $sinking, $open_bottom" begin - test_NPZD(grid, sinking, open_bottom) - end +grid = RectilinearGrid(architecture; size=(3, 3, 6), extent=(1, 1, 2)) + +for sinking = (false, true), open_bottom = (false, true) + if !(sinking && open_bottom) # no sinking is the same with and without open bottom + @info "Testing on $(typeof(arch)) with sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" + @testset "$arch, $sinking, $open_bottom" begin + test_NPZD(grid, sinking, open_bottom) end end end diff --git a/test/test_light.jl b/test/test_light.jl index fdcf99169..831c9dc8a 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -42,13 +42,11 @@ function test_two_band(grid, bgc, model_type) return all(results_PAR .≈ reverse(expected_PAR)) end -archs = (CPU(), ) @testset "Light attenuaiton model" begin for model in (NonhydrostaticModel, HydrostaticFreeSurfaceModel), - arch in archs, - grid in (RectilinearGrid(arch; size = (2, 2, 2), extent = (2, 2, 2)), - LatitudeLongitudeGrid(arch; size = (5, 5, 2), longitude = (-180, 180), latitude = (-85, 85), z = (-2, 0))), + grid in (RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)), + LatitudeLongitudeGrid(architecture; size = (5, 5, 2), longitude = (-180, 180), latitude = (-85, 85), z = (-2, 0))), bgc in (LOBSTER, NutrientPhytoplanktonZooplanktonDetritus) # this is now redundant since each model doesn't deal with the light separatly if !((model == NonhydrostaticModel) && ((grid isa LatitudeLongitudeGrid) | (grid isa OrthogonalSphericalShellGrid))) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 73c7f940a..14f3d2751 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -113,30 +113,28 @@ display_name(::ImmersedBoundaryGrid) = "Immersed boundary grid" bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill @testset "Sediment integration" begin - for architecture in (CPU(), ) - grids = [RectilinearGrid(architecture; size=(3, 3, 50), extent=(10, 10, 500)), - LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), - ImmersedBoundaryGrid( - LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), - GridFittedBottom(bottom_height))] - for grid in grids - for timestepper in (:QuasiAdamsBashforth2, :RungeKutta3), - sediment_model in (InstantRemineralisation(; grid), SimpleMultiG(; grid)), - model in (NonhydrostaticModel, HydrostaticFreeSurfaceModel) - for biogeochemistry in (NutrientPhytoplanktonZooplanktonDetritus(; grid, sediment_model), - 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)) - # get rid of incompatible combinations - run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || - (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || - (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) - if run - @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry, model; timestepper) - end + grids = [RectilinearGrid(architecture; size=(3, 3, 50), extent=(10, 10, 500)), + LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), + ImmersedBoundaryGrid( + LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), + GridFittedBottom(bottom_height))] + for grid in grids + for timestepper in (:QuasiAdamsBashforth2, :RungeKutta3), + sediment_model in (InstantRemineralisation(; grid), SimpleMultiG(; grid)), + model in (NonhydrostaticModel, HydrostaticFreeSurfaceModel) + for biogeochemistry in (NutrientPhytoplanktonZooplanktonDetritus(; grid, sediment_model), + 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)) + # get rid of incompatible combinations + run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || + (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || + (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) + if run + @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry, model; timestepper) end end end diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index ada407cd9..1483c89fd 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -10,8 +10,7 @@ function intercept_tendencies!(model, intercepted_tendencies) end @testset "SLatissima particle setup and conservations" begin - arch = CPU() - grid = RectilinearGrid(arch; size=(1, 1, 1), extent=(1, 1, 1)) + grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) # Initial properties diff --git a/test/test_utils.jl b/test/test_utils.jl index 86586a791..6baab9730 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -46,9 +46,7 @@ function test_negative_zeroing(arch) end @testset "Test Utils" begin - for arch in (CPU(), ) - @test test_column_diffusion_timescale(arch) - @test test_negative_scaling(arch) - @test test_negative_zeroing(arch) - end + @test test_column_diffusion_timescale(architecture) + @test test_negative_scaling(architecture) + @test test_negative_zeroing(architecture) end From ce514c0ae6fdbd29016579df0fc708721cb42ce8 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 11 Sep 2023 15:55:40 +0100 Subject: [PATCH 009/112] oops --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 69a36e363..e758c0874 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using OceanBioME, Documenter, Test, CUDA +using OceanBioME, Documenter, Test, CUDA, Oceananigans architecture = CUDA.has_cuda() ? GPU() : CPU() From 67f9365b92a3eb2b883a776f8cc787fe5af2804f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 11 Sep 2023 15:57:21 +0100 Subject: [PATCH 010/112] fixed test labelling --- test/test_LOBSTER.jl | 4 ++-- test/test_NPZD.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index c7390d616..2c4cfe305 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -118,8 +118,8 @@ grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 2)) for open_bottom = (false, true), sinking = (false, true), variable_redfield = (false, true), oxygen = (false, true), carbonates = (false, true) if !(sinking && open_bottom) # no sinking is the same with and without open bottom - @info "Testing on $(typeof(arch)) with carbonates $(carbonates ? :✅ : :❌), oxygen $(oxygen ? :✅ : :❌), variable redfield $(variable_redfield ? :✅ : :❌), sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" - @testset "$arch, $carbonates, $oxygen, $variable_redfield, $sinking, $open_bottom" begin + @info "Testing on $(typeof(architecture)) with carbonates $(carbonates ? :✅ : :❌), oxygen $(oxygen ? :✅ : :❌), variable redfield $(variable_redfield ? :✅ : :❌), sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" + @testset "$architecture, $carbonates, $oxygen, $variable_redfield, $sinking, $open_bottom" begin test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open_bottom, n_timesteps) end end diff --git a/test/test_NPZD.jl b/test/test_NPZD.jl index 945070fb1..d7e8d156a 100644 --- a/test/test_NPZD.jl +++ b/test/test_NPZD.jl @@ -51,8 +51,8 @@ grid = RectilinearGrid(architecture; size=(3, 3, 6), extent=(1, 1, 2)) for sinking = (false, true), open_bottom = (false, true) if !(sinking && open_bottom) # no sinking is the same with and without open bottom - @info "Testing on $(typeof(arch)) with sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" - @testset "$arch, $sinking, $open_bottom" begin + @info "Testing on $(typeof(architecture)) with sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" + @testset "$architecture, $sinking, $open_bottom" begin test_NPZD(grid, sinking, open_bottom) end end From ff0e34311f72fcb74d58fba8f36f6da6ea845727 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 13:18:33 +0100 Subject: [PATCH 011/112] maybe fixed scaling --- src/Utils/negative_tracers.jl | 39 +++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index a49229adc..01cee412f 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -39,7 +39,7 @@ struct ScaleNegativeTracers{FA, SA, W} warn :: W ScaleNegativeTracers(tracers::FA, scalefactors::SA, warn::W) where {FA, SA, W} = - new{FA, SA, W}(tracers, scalefactors, warn) + warn ? error("Warning not currently implemented") : new{FA, SA, W}(tracers, scalefactors, warn) end adapt_structure(to, snt::ScaleNegativeTracers) = ScaleNegativeTracers(adapt(to, snt.tracers), @@ -62,7 +62,7 @@ Future plans include implement a positivity-preserving timestepping scheme as th If `warn` is true then scaling will raise a warning. """ -function ScaleNegativeTracers(tracers; scalefactors = NamedTuple{tracers}(ones(length(tracers))), warn = false) +function ScaleNegativeTracers(tracers; scalefactors = ones(length(tracers)), warn = false) if length(scalefactors) != length(tracers) error("Incorrect number of scale factors provided") end @@ -79,40 +79,49 @@ If `warn` is true then scaling will raise a warning. """ function ScaleNegativeTracers(model::AbstractBiogeochemistry; warn = false) tracers = conserved_tracers(model) - scalefactors = NamedTuple{tracers}(ones(length(tracers))) + scalefactors = ones(length(tracers)) return ScaleNegativeTracers(tracers, scalefactors, warn) end function update_biogeochemical_state!(model, scale::ScaleNegativeTracers) workgroup, worksize = work_layout(model.grid, :xyz) + dev = device(model.grid.architecture) + scale_for_negs_kernel! = scale_for_negs!(dev, workgroup, worksize) - scale_for_negs_kernel!(model.tracers, scale.tracers, scale.scalefactors) + + tracers_to_scale = Tuple(model.tracers[tracer_name] for tracer_name in keys(scale.tracers)) + + scale_for_negs_kernel!(tracers_to_scale, scale.scalefactors) end -@kernel function scale_for_negs!(fields, tracers, scalefactors) +@kernel function scale_for_negs!(tracers, scalefactors) i, j, k = @index(Global, NTuple) + t, p = 0.0, 0.0 + @unroll for (idx, tracer) in enumerate(tracers) - field = @inbounds fields[tracer][i, j, k] - scalefactor = @inbounds scalefactors[tracer] + value = @inbounds tracer[i, j, k] + scalefactor = @inbounds scalefactors[idx] - t += field * scalefactor - if field > 0 - p += field * scalefactor + t += value * scalefactor + if value > 0 + p += value * scalefactor end end + t < 0 && error("Cell total < 0, can not scale negative tracers.") + @unroll for tracer in tracers - field = @inbounds fields[tracer][i, j, k] + value = @inbounds tracer[i, j, k] - if field > 0 - field *= t / p + if value > 0 + value *= t / p else - field = 0 + value = 0 end - @inbounds fields[tracer][i, j, k] = field + @inbounds tracer[i, j, k] = value end end \ No newline at end of file From 3c8edb6d6c06ab316139abbb0e0dc4a5fe68c28c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 15:27:58 +0100 Subject: [PATCH 012/112] tidying --- test/test_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index 6baab9730..91aace041 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -14,7 +14,7 @@ function test_column_diffusion_timescale(arch) carbonates = true), advection = nothing) - return column_diffusion_timescale(model) ≈ min_Δz^2 / κ + return column_diffusion_timescale(model) ≈ min_Δz ^ 2 / κ end function test_negative_scaling(arch) From e0d4e7d668d9bf52ce949aee9f15a5b3151e92ee Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 17:36:53 +0100 Subject: [PATCH 013/112] idea --- src/Boundaries/Sediments/Sediments.jl | 2 +- src/OceanBioME.jl | 9 +++------ src/Utils/timestep.jl | 1 + 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index c322930e3..8f80f958f 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -37,7 +37,7 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - bgc.sediment, bgc, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ) + sediment, bgc.underlying_biogeochemistry, model.grid, model.advection, model.tracers, model.timestepper) return nothing end diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index a7b872642..d11c517dd 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -102,11 +102,7 @@ biogeochemical_drift_velocity(bgc::Biogeochemistry, val_name) = biogeochemical_d biogeochemical_auxiliary_fields(bgc::Biogeochemistry) = merge(biogeochemical_auxiliary_fields(bgc.underlying_biogeochemistry), biogeochemical_auxiliary_fields(bgc.light_attenuation)) -adapt_structure(to, bgc::Biogeochemistry) = Biogeochemistry(adapt(to, bgc.underlying_biogeochemistry), - adapt(to, bgc.light_attenuation), - adapt(to, bgc.sediment), - adapt(to, bgc.particles), - adapt(to, bgc.modifiers)) +@inline adapt_structure(to, bgc::Biogeochemistry) = adapt(to, bgc.underlying_biogeochemistry) function update_tendencies!(bgc::Biogeochemistry, model) update_tendencies!(bgc, bgc.sediment, model) @@ -117,8 +113,9 @@ end update_tendencies!(bgc, modifier, model) = nothing update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modifier, model) for modifier in modifiers] +# do we still need this for CPU kernels??? @inline biogeochemical_transition(i, j, k, grid, bgc::Biogeochemistry, val_tracer_name, clock, fields) = - biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) + @show biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) function update_biogeochemical_state!(bgc::Biogeochemistry, model) update_biogeochemical_state!(model, bgc.modifiers) diff --git a/src/Utils/timestep.jl b/src/Utils/timestep.jl index d9d6096e7..910b88064 100644 --- a/src/Utils/timestep.jl +++ b/src/Utils/timestep.jl @@ -1,6 +1,7 @@ using Oceananigans.Advection: cell_advection_timescale using Oceananigans.Grids: Center, znodes, zspacing, minimum_zspacing +# NOT GPU SAFE @inline function column_diffusion_timescale(model) grid = model.grid From b32f4fd58a7245a8b75bd97ceaf25dc19ced6d2e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 17:44:32 +0100 Subject: [PATCH 014/112] particles fix --- src/Particles/Particles.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index a9f28472b..3c9e4d5eb 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -24,7 +24,7 @@ abstract type BiogeochemicalParticles end Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) @inline update_lagrangian_particle_properties!(model, bgc::Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, Δt) = - update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) + update_lagrangian_particle_properties!(bgc.particles, model, bgc.underlying_biogeochemistry, Δt) update_lagrangian_particle_properties!(::BiogeochemicalParticles, model, bgc, Δt) = nothing From df20a2f920b3833aedc74d3d2ae342c90402d6ea Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 17:57:17 +0100 Subject: [PATCH 015/112] maybe fixed --- src/Models/Individuals/SLatissima.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 1d8ead7c3..c05a8f53e 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -369,7 +369,7 @@ end t = clock.time @inbounds if p.A[idx] > 0 - NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers) + NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(bgc)) photo = photosynthesis(T, PAR, p) e = exudation(C, p) @@ -510,10 +510,8 @@ end @inline normed_day_length_change(n, ϕ) = (day_length(n, ϕ) - day_length(n - 1, ϕ)) / (day_length(76, ϕ) - day_length(75, ϕ)) -@inline function get_arguments(x, y, z, t, particles, bgc, grid, velocities, tracers) - +@inline function get_arguments(x, y, z, t, particles, bgc, grid, velocities, tracers, auxiliary_fields) bgc_tracers = required_biogeochemical_tracers(bgc) - auxiliary_fields = biogeochemical_auxiliary_fields(bgc) i, j, k = fractional_indices(x, y, z, (Center(), Center(), Center()), grid) From 5d4944f5f11595989fb3a23e1821a5148678343f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 18:01:19 +0100 Subject: [PATCH 016/112] a fix --- src/Utils/negative_tracers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 01cee412f..c6366d187 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -2,7 +2,7 @@ using Oceananigans: fields, Simulation using KernelAbstractions using KernelAbstractions.Extras.LoopInfo: @unroll using Oceananigans.Utils: work_layout -using Oceananigans.Architectures: device +using Oceananigans.Architectures: device, architecture, arch_array using Oceananigans.Biogeochemistry: AbstractBiogeochemistry import Adapt: adapt_structure, adapt @@ -79,7 +79,7 @@ If `warn` is true then scaling will raise a warning. """ function ScaleNegativeTracers(model::AbstractBiogeochemistry; warn = false) tracers = conserved_tracers(model) - scalefactors = ones(length(tracers)) + scalefactors = arch_array(architecture(model), ones(length(tracers))) return ScaleNegativeTracers(tracers, scalefactors, warn) end From 6b40a921fbab3a23ebf3a18e30b2d79add43c18c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 18:11:44 +0100 Subject: [PATCH 017/112] kelp please work --- src/Models/Individuals/SLatissima.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index c05a8f53e..91cb430f7 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -369,7 +369,7 @@ end t = clock.time @inbounds if p.A[idx] > 0 - NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(bgc)) + NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(bgc).PAR) photo = photosynthesis(T, PAR, p) e = exudation(C, p) @@ -510,7 +510,7 @@ end @inline normed_day_length_change(n, ϕ) = (day_length(n, ϕ) - day_length(n - 1, ϕ)) / (day_length(76, ϕ) - day_length(75, ϕ)) -@inline function get_arguments(x, y, z, t, particles, bgc, grid, velocities, tracers, auxiliary_fields) +@inline function get_arguments(x, y, z, t, particles, bgc, grid, velocities, tracers, PAR_field) bgc_tracers = required_biogeochemical_tracers(bgc) i, j, k = fractional_indices(x, y, z, (Center(), Center(), Center()), grid) @@ -522,7 +522,7 @@ end # TODO: ADD ALIASING/RENAMING OF TRACERS SO WE CAN USE E.G. N IN STEAD OF NO3 NO₃ = _interpolate(tracers.NO₃, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) - PAR = _interpolate(auxiliary_fields.PAR, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) * day / (3.99e-10 * 545e12) # W / m² / s to einstein / m² / day + PAR = _interpolate(PAR_field, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) * day / (3.99e-10 * 545e12) # W / m² / s to einstein / m² / day if :NH₄ in bgc_tracers NH₄ = _interpolate(tracers.NH₄, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) From 834d18ca94130193c307303efcd6658ad8da480a Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 18:15:44 +0100 Subject: [PATCH 018/112] oops wrong model --- src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl | 2 +- src/Models/AdvectedPopulations/NPZD.jl | 2 +- src/Utils/negative_tracers.jl | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 69fdbba6c..31fbee405 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -352,7 +352,7 @@ function LOBSTER(; grid, sinking_velocities) if scale_negatives - scaler = ScaleNegativeTracers(underlying_biogeochemistry) + scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid) modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler) end diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index ed892a72e..5d54b43f9 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -199,7 +199,7 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, sinking_velocities) if scale_negatives - scaler = ScaleNegativeTracers(underlying_biogeochemistry) + scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid) modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler) end diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index c6366d187..0e8b8a003 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -77,9 +77,9 @@ Construct a modifier to scale the conserved tracers in `model`. If `warn` is true then scaling will raise a warning. """ -function ScaleNegativeTracers(model::AbstractBiogeochemistry; warn = false) - tracers = conserved_tracers(model) - scalefactors = arch_array(architecture(model), ones(length(tracers))) +function ScaleNegativeTracers(bgc::AbstractBiogeochemistry, grid; warn = false) + tracers = conserved_tracers(bgc) + scalefactors = arch_array(architecture(grid), ones(length(tracers))) return ScaleNegativeTracers(tracers, scalefactors, warn) end From c5c956934cf7957366cc11fb671b664d5b6cb25b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 18:52:41 +0100 Subject: [PATCH 019/112] fixed utils test --- test/test_utils.jl | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index 91aace041..233c5dfdf 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -28,7 +28,10 @@ function test_negative_scaling(arch) run!(simulation) - return (model.tracers.N[1, 1, 1] ≈ 1) && (model.tracers.P[1, 1, 1] ≈ 0.0) + N = Array(interior(model.tracers.N))[1, 1, 1] + P = Array(interior(model.tracers.P))[1, 1, 1] + + return (N ≈ 1) && (P ≈ 0.0) end function test_negative_zeroing(arch) @@ -42,11 +45,15 @@ function test_negative_zeroing(arch) run!(simulation) - return (model.tracers.N[1, 1, 1] ≈ 2) && (model.tracers.P[1, 1, 1] ≈ 0.0) && (model.tracers.Z[1, 1, 1] ≈ -1) + N = Array(interior(model.tracers.N))[1, 1, 1] + P = Array(interior(model.tracers.P))[1, 1, 1] + Z = Array(interior(model.tracers.Z))[1, 1, 1] + + return (N ≈ 2) && (P ≈ 0.0) && (Z ≈ -1) end @testset "Test Utils" begin - @test test_column_diffusion_timescale(architecture) + @test test_column_diffusion_timescale(architecture) broken = architecture == GPU() @test test_negative_scaling(architecture) - @test test_negative_zeroing(architecture) + @test test_negative_zeroing(architecture) broken = architecture == GPU() end From 1065f22843b663b59c946bbc40834721cd23a521 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 18:59:20 +0100 Subject: [PATCH 020/112] not sure why this won't work --- src/Models/Individuals/SLatissima.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 91cb430f7..0004387e8 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -347,13 +347,13 @@ function update_lagrangian_particle_properties!(particles::SLatissima, model, bg update_particle_properties_kernel! = _update_lagrangian_particle_properties!(device(arch), workgroup, worksize) - update_particle_properties_kernel!(particles, bgc, model.grid, + update_particle_properties_kernel!(particles, light_attenuation, bgc.underlying, model.grid, model.velocities, model.tracers, model.clock, Δt) particles.custom_dynamics(particles, model, bgc, Δt) end -@kernel function _update_lagrangian_particle_properties!(p, bgc, grid, velocities, tracers, clock, Δt) +@kernel function _update_lagrangian_particle_properties!(p, light_attenuation, bgc, grid, velocities, tracers, clock, Δt) idx = @index(Global) @inbounds begin @@ -369,7 +369,7 @@ end t = clock.time @inbounds if p.A[idx] > 0 - NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(bgc).PAR) + NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(light_attenuation).PAR) photo = photosynthesis(T, PAR, p) e = exudation(C, p) From a149f8f99a194f699b42b8647c44f450c2ec6adc Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 19:00:57 +0100 Subject: [PATCH 021/112] that one works apparently --- test/test_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index 233c5dfdf..8f33d8bc9 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -55,5 +55,5 @@ end @testset "Test Utils" begin @test test_column_diffusion_timescale(architecture) broken = architecture == GPU() @test test_negative_scaling(architecture) - @test test_negative_zeroing(architecture) broken = architecture == GPU() + @test test_negative_zeroing(architecture) end From c090448b8da5d218e77b7612e08839936b036c6d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 19:03:33 +0100 Subject: [PATCH 022/112] oops --- src/Models/Individuals/SLatissima.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 0004387e8..b0f1c8e27 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -347,7 +347,7 @@ function update_lagrangian_particle_properties!(particles::SLatissima, model, bg update_particle_properties_kernel! = _update_lagrangian_particle_properties!(device(arch), workgroup, worksize) - update_particle_properties_kernel!(particles, light_attenuation, bgc.underlying, model.grid, + update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying, model.grid, model.velocities, model.tracers, model.clock, Δt) particles.custom_dynamics(particles, model, bgc, Δt) From ec50cf5d74687c3fdb61a8859d71cf845c7f4035 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 19:07:41 +0100 Subject: [PATCH 023/112] ah thats it --- src/Models/Individuals/SLatissima.jl | 2 +- src/Particles/Particles.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index b0f1c8e27..839a3f6ce 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -347,7 +347,7 @@ function update_lagrangian_particle_properties!(particles::SLatissima, model, bg update_particle_properties_kernel! = _update_lagrangian_particle_properties!(device(arch), workgroup, worksize) - update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying, model.grid, + update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying_biogeochemistry, model.grid, model.velocities, model.tracers, model.clock, Δt) particles.custom_dynamics(particles, model, bgc, Δt) diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 3c9e4d5eb..a9f28472b 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -24,7 +24,7 @@ abstract type BiogeochemicalParticles end Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) @inline update_lagrangian_particle_properties!(model, bgc::Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, Δt) = - update_lagrangian_particle_properties!(bgc.particles, model, bgc.underlying_biogeochemistry, Δt) + update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) update_lagrangian_particle_properties!(::BiogeochemicalParticles, model, bgc, Δt) = nothing From b658864620fe9aa9ea7388d933b512cca447f120 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 14 Sep 2023 19:21:44 +0100 Subject: [PATCH 024/112] fixed light test --- test/test_light.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_light.jl b/test/test_light.jl index 831c9dc8a..3b274be3e 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -37,7 +37,7 @@ function test_two_band(grid, bgc, model_type) expected_PAR = 100.0 .* [exp(- 0.5 * kʳ - ∫Chlʳ[1] * χʳ) + exp(- 0.5 * kᵇ - ∫Chlᵇ[1] * χᵇ), exp(- 1.5 * kʳ - ∫Chlʳ[2] * χʳ) + exp(- 1.5 * kᵇ - ∫Chlᵇ[2] * χᵇ)] ./ 2 - results_PAR = convert(Array, biogeochemical_auxiliary_fields(biogeochemistry).PAR)[1, 1, 1:2] + results_PAR = Array(interior(biogeochemical_auxiliary_fields(biogeochemistry).PAR))[1, 1, 1:2] return all(results_PAR .≈ reverse(expected_PAR)) end From a0d8da43b2396598901e348252f211b36e5a6a12 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 14:14:53 +0100 Subject: [PATCH 025/112] fixed npzd and losbter tests --- test/test_LOBSTER.jl | 2 +- test/test_NPZD.jl | 14 ++++++++++---- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index 2c4cfe305..df86a1e09 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -119,7 +119,7 @@ grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 2)) for open_bottom = (false, true), sinking = (false, true), variable_redfield = (false, true), oxygen = (false, true), carbonates = (false, true) if !(sinking && open_bottom) # no sinking is the same with and without open bottom @info "Testing on $(typeof(architecture)) with carbonates $(carbonates ? :✅ : :❌), oxygen $(oxygen ? :✅ : :❌), variable redfield $(variable_redfield ? :✅ : :❌), sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" - @testset "$architecture, $carbonates, $oxygen, $variable_redfield, $sinking, $open_bottom" begin + @testset "LOBSTER $architecture, $carbonates, $oxygen, $variable_redfield, $sinking, $open_bottom" begin test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open_bottom, n_timesteps) end end diff --git a/test/test_NPZD.jl b/test/test_NPZD.jl index d7e8d156a..e11503feb 100644 --- a/test/test_NPZD.jl +++ b/test/test_NPZD.jl @@ -26,7 +26,7 @@ function test_NPZD(grid, sinking, open_bottom) time_step!(model, 1.0) # and that they all return zero - @test all([all(values .== 0) for values in values(model.tracers)]) + @test all([all(Array(interior(values)) .== 0) for values in values(model.tracers)]) # mass conservation model.tracers.N .= rand() @@ -34,13 +34,19 @@ function test_NPZD(grid, sinking, open_bottom) model.tracers.Z .= rand() model.tracers.D .= rand() - ΣN₀ = sum(model.tracers.N) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.D) + ΣN₀ = sum(Array(interior(model.tracers.N))) + + sum(Array(interior(model.tracers.P))) + + sum(Array(interior(model.tracers.Z))) + + sum(Array(interior(model.tracers.D))) for n in 1:1000 time_step!(model, 1.0) end - ΣN₁ = sum(model.tracers.N) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.D) + ΣN₁ = sum(Array(interior(model.tracers.N))) + + sum(Array(interior(model.tracers.P))) + + sum(Array(interior(model.tracers.Z))) + + sum(Array(interior(model.tracers.D))) @test ΣN₀ ≈ ΣN₁ # guess this should actually fail with a high enough accuracy when sinking is on with an open bottom @@ -52,7 +58,7 @@ grid = RectilinearGrid(architecture; size=(3, 3, 6), extent=(1, 1, 2)) for sinking = (false, true), open_bottom = (false, true) if !(sinking && open_bottom) # no sinking is the same with and without open bottom @info "Testing on $(typeof(architecture)) with sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" - @testset "$architecture, $sinking, $open_bottom" begin + @testset "NPZD $architecture, $sinking, $open_bottom" begin test_NPZD(grid, sinking, open_bottom) end end From 058619b67e8505189dc8741d5ef7139b872de707 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 15:12:14 +0100 Subject: [PATCH 026/112] hunch about memory not being freed --- src/OceanBioME.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index d11c517dd..afd6b37f7 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -32,7 +32,9 @@ export ScaleNegativeTracers, ZeroNegativeTracers export ColumnField, isacolumn using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry +using Oceananigans.Architectures: architecture, device using Adapt +using KernelAbstractions: synchronize import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, @@ -108,6 +110,7 @@ function update_tendencies!(bgc::Biogeochemistry, model) update_tendencies!(bgc, bgc.sediment, model) update_tendencies!(bgc, bgc.particles, model) update_tendencies!(bgc, bgc.modifiers, model) + synchronize(device(architecture(model))) end update_tendencies!(bgc, modifier, model) = nothing From 1ff75e54b537d47ce1091961c8df92f06c297547 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 15:18:20 +0100 Subject: [PATCH 027/112] forgot to remove a show (but also revealed that a test wasn't using GPU --- src/OceanBioME.jl | 2 +- test/test_gasexchange.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index afd6b37f7..3fdf491ad 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -118,7 +118,7 @@ update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modi # do we still need this for CPU kernels??? @inline biogeochemical_transition(i, j, k, grid, bgc::Biogeochemistry, val_tracer_name, clock, fields) = - @show biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) + biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) function update_biogeochemical_state!(bgc::Biogeochemistry, model) update_biogeochemical_state!(model, bgc.modifiers) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 0c3dcad08..be7cd442b 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -52,7 +52,7 @@ end end @testset "Gas exchange coupling" begin - grid = RectilinearGrid(size=(1, 1, 2), extent=(1, 1, 1)) + grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) conc_field = CenterField(grid, indices=(:, :, grid.Nz)) conc_field .= 413.0 + 1.0 * rand() conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / (1year)) From 40f0d28f1d6e1005470bb1269bb3a9c1f4eef1d0 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 15:59:43 +0100 Subject: [PATCH 028/112] removed scalar indexing from kelp tests --- test/test_slatissima.jl | 40 ++++++++++++++++++++++++++-------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index 1483c89fd..69f1462ee 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -5,10 +5,24 @@ using Oceananigans.Fields: TracerFields function intercept_tendencies!(model, intercepted_tendencies) for tracer in keys(model.tracers) - copyto!(intercepted_tendencies[tracer], model.timestepper.Gⁿ[tracer]) + intercepted_tendencies[tracer] = Array(interior(model.timestepper.Gⁿ[tracer])) end end +sum_tracer_nitrogen(tracers) = sum(Array(interior(tracers.NO₃))) + + sum(Array(interior(tracers.NH₄))) + + sum(Array(interior(tracers.P))) + + sum(Array(interior(tracers.Z))) + + sum(Array(interior(tracers.sPON))) + + sum(Array(interior(tracers.bPON))) + + sum(Array(interior(tracers.DON))) + +sum_tracer_carbon(tracers, redfield) = sum(Array(interior(tracers.sPOC))) + + sum(Array(interior(tracers.bPOC))) + + sum(Array(interior(tracers.DOC))) + + sum(Array(interior(tracers.DIC))) + + sum(Array(interior(tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) + + model.tracers.Z))) * redfield @testset "SLatissima particle setup and conservations" begin grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) @@ -35,13 +49,11 @@ end set!(model, NO₃ = 10.0, NH₄ = 1.0, DIC = 2000, Alk = 2000) - initial_tracer_N = sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.sPON) + sum(model.tracers.bPON) + sum(model.tracers.DON) - initial_kelp_N = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.N .+ particles.structural_nitrogen)) ./ (14 * 0.001) + initial_tracer_N = sum_tracer_nitrogen(model.tracers) + initial_kelp_N = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.N) .+ particles.structural_nitrogen)) ./ (14 * 0.001) - initial_tracer_C = sum(model.tracers.sPOC) + sum(model.tracers.bPOC) + sum(model.tracers.DOC) + sum(model.tracers.DIC) + - sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield - - initial_kelp_C = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.C .+ particles.structural_carbon)) ./ (12 * 0.001) + initial_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield) + initial_kelp_C = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.C) .+ particles.structural_carbon)) ./ (12 * 0.001) model.clock.time = 60days # get to a high growth phase @@ -49,12 +61,11 @@ end time_step!(model, 1.0) end - final_tracer_N = sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.sPON) + sum(model.tracers.bPON) + sum(model.tracers.DON) - final_kelp_N = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.N .+ particles.structural_nitrogen)) ./ (14 * 0.001) + final_tracer_N = sum_tracer_nitrogen(model.tracers) + final_kelp_N = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.N) .+ particles.structural_nitrogen)) ./ (14 * 0.001) - final_tracer_C = sum(model.tracers.sPOC) + sum(model.tracers.bPOC) + sum(model.tracers.DOC) + sum(model.tracers.DIC) + - sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield - final_kelp_C = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.C .+ particles.structural_carbon)) ./ (12 * 0.001) + final_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield) + final_kelp_C = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.C) .+ particles.structural_carbon)) ./ (12 * 0.001) # kelp is being integrated @test initial_kelp_N != final_kelp_N @@ -66,12 +77,13 @@ end simulation = Simulation(model, Δt = 1.0, stop_iteration = 1) - intercepted_tendencies = TracerFields(keys(model.tracers), grid) + # slow but easiest to have this just done on CPU + intercepted_tendencies = (Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) run!(simulation) # the model is changing the tracer tendencies - not sure this test actually works as it didn't fail when it should have - @test any([any(intercepted_tendencies[tracer] .!= model.timestepper.Gⁿ[tracer]) for tracer in (:NO₃, :NH₄, :DIC, :DOC, :bPON, :bPOC)]) + @test any([any(intercepted_tendencies[tracer] .!= Array(interior(model.timestepper.Gⁿ[tracer]))) for tracer in (:NO₃, :NH₄, :DIC, :DOC, :bPON, :bPOC)]) end From 64fc80efd2eb6891c84fb0eb3144cbb20ec196d9 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 16:01:52 +0100 Subject: [PATCH 029/112] removed scalar indexing from sediment tests --- test/test_sediments.jl | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 14f3d2751..f988e6470 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -46,15 +46,30 @@ set_defaults!(::lobster_variable_redfield, model) = set_defaults!(::NutrientPhytoplanktonZooplanktonDetritus, model) = set!(model, N = 2.3, P = 0.4, Z = 0.5, D = 0.2) -total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) + - sum(sed.fields.N_slow) + - sum(sed.fields.N_ref) - -total_nitrogen(sed::InstantRemineralisation) = sum(sed.fields.N_storage) - -total_nitrogen(::LOBSTER, model) = sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.DOM) + sum(model.tracers.sPOM) + sum(model.tracers.bPOM) -total_nitrogen(::lobster_variable_redfield, 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) -total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(model.tracers.N) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.D) +total_nitrogen(sed::SimpleMultiG) = sum(Array(interior(sed.fields.N_fast))) + + sum(Array(interior(sed.fields.N_slow))) + + sum(Array(interior(sed.fields.N_ref))) + +total_nitrogen(sed::InstantRemineralisation) = sum(Array(interior(sed.fields.N_storage))) + +total_nitrogen(::LOBSTER, model) = sum(Array(interior(model.tracers.NO₃))) + + sum(Array(interior(model.tracers.NH₄))) + + sum(Array(interior(model.tracers.P))) + + sum(Array(interior(model.tracers.Z))) + + sum(Array(interior(model.tracers.DOM))) + + sum(Array(interior(model.tracers.sPOM))) + + sum(Array(interior(model.tracers.bPOM) +total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tracers.NO₃))) + + sum(Array(interior(model.tracers.NH₄))) + + sum(Array(interior(model.tracers.P))) + + sum(Array(interior(model.tracers.Z))) + + sum(Array(interior(model.tracers.DON))) + + sum(Array(interior(model.tracers.sPON))) + + sum(Array(interior(model.tracers.bPON) +total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(Array(interior(model.tracers.N))) + + sum(Array(interior(model.tracers.P))) + + sum(Array(interior(model.tracers.Z))) + + sum(Array(interior(model.tracers.D) function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAdamsBashforth2) model = isa(model, NonhydrostaticModel) ? model(; grid, From e0ceb49331e763567b14ad574b767f35222a13b0 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 16:16:44 +0100 Subject: [PATCH 030/112] fixed (?) gas exchange test --- test/test_gasexchange.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index be7cd442b..097465541 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -24,8 +24,10 @@ function test_gas_exchange_model(grid, air_concentration) set!(model, T = 15.0, S = 35.0, DIC = 2220, Alk = 2500) - # is everything communicating properly? - @test Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1, 1, grid, model.clock, fields(model)) ≈ -0.0003 atol = 0.0001 + # is everything communicating properly? (can't think of a way to not use allow scalar here) + value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1, 1, grid, model.clock, fields(model)) + + @test value ≈ -0.0003 atol = 0.0001 # just incase we broke something @test isnothing(time_step!(model, 1.0)) From 49282b317277d120a0fcc24abe587eaba65ab3b3 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 16:20:08 +0100 Subject: [PATCH 031/112] don't run doctests on GPU --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index e758c0874..6a721a0e4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,6 @@ include("test_gasexchange.jl") include("test_slatissima.jl") include("test_sediments.jl") -@testset "Doctests" begin +architecture == CPU() && @testset "Doctests" begin doctest(OceanBioME) end From 8c880c34ffec460baeca9160bb58850ca727d73f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 16:21:44 +0100 Subject: [PATCH 032/112] gas exchange test still not running --- test/test_gasexchange.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 097465541..7cf556d12 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -53,11 +53,12 @@ end @test (mean(pCO₂_err) < 10 && std(pCO₂_err) < 15) end +@inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / (1year)) + @testset "Gas exchange coupling" begin grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) conc_field = CenterField(grid, indices=(:, :, grid.Nz)) conc_field .= 413.0 + 1.0 * rand() - conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / (1year)) for air_concentration in [413.1, conc_function, conc_field] @info "Testing with $(typeof(air_concentration))" test_gas_exchange_model(grid, air_concentration) From 956734fde5d6d870084f2f0c9ba83720a082f46e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Fri, 15 Sep 2023 16:37:47 +0100 Subject: [PATCH 033/112] fix gas exchange test --- test/test_gasexchange.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 7cf556d12..d3ddf6443 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -3,7 +3,7 @@ using OceanBioME: Boundaries, GasExchange, LOBSTER using Oceananigans, DataDeps, JLD2, Statistics using Oceananigans.Units -year = years = 365days # just for the idealised case below +const year = years = 365days # just for the idealised case below dd = DataDep( "test_data", @@ -53,7 +53,7 @@ end @test (mean(pCO₂_err) < 10 && std(pCO₂_err) < 15) end -@inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / (1year)) +@inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / year) @testset "Gas exchange coupling" begin grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) From e414c19908936dee46176fd9fba5ab59dd916d60 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 16 Sep 2023 18:10:33 +0100 Subject: [PATCH 034/112] fixed gas exchange on fields --- src/Boundaries/gasexchange.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index badb7610a..3cde0e6e1 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -3,6 +3,8 @@ ##### As per OCMIP Phase 2 http://ocmip5.ipsl.jussieu.fr/OCMIP/phase2/simulations/Abiotic/Cchem/co2calc.f simplified as in PISCESv2 ##### +using Oceananigans.Fields: interpolate + struct pCO₂{P0, P1, P2, PB, PW, FT} solubility :: P0 bicarbonate_dissociation :: P1 @@ -257,9 +259,4 @@ end @inline get_value(x, y, t, air_concentration::Number) = air_concentration @inline get_value(x, y, t, air_concentration::Function) = air_concentration(x, y, t) - -# interpolation does not work on 2d grids as trilinear interpolation fails but we're only ever going to be asking for values on grid points -@inline function get_value(x, y, t, air_concentration::Field{Center, Center, Center}) - (ξ, i), (η, j), (ζ, k) = modf.(fractional_indices(x, y, znodes(air_concentration)[end], (Center(), Center(), Center()), air_concentration.grid)) - return air_concentration[floor(Int, i) + 1, floor(Int, j) + 1, floor(Int, k) + 1] -end \ No newline at end of file +@inline get_value(x, y, t, conc::Field{Center, Center, Center}) = interpolate(conc, Center(), Center(), Center(), conc.grid, x, y, 0) \ No newline at end of file From de5f4d3898d0d39738f6deb5812a2da9198658f0 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 16 Sep 2023 18:27:03 +0100 Subject: [PATCH 035/112] fix gas exchange test --- test/test_gasexchange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index d3ddf6443..573ce8316 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -25,7 +25,7 @@ function test_gas_exchange_model(grid, air_concentration) set!(model, T = 15.0, S = 35.0, DIC = 2220, Alk = 2500) # is everything communicating properly? (can't think of a way to not use allow scalar here) - value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1, 1, grid, model.clock, fields(model)) + value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1.0, 1.0, grid, model.clock, fields(model)) @test value ≈ -0.0003 atol = 0.0001 From 3d409b8a109da5d0535fa3aff5379e1412e85844 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 16 Sep 2023 18:47:20 +0100 Subject: [PATCH 036/112] fixed again --- src/Boundaries/gasexchange.jl | 2 +- test/test_gasexchange.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index 3cde0e6e1..d8e4f3781 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -259,4 +259,4 @@ end @inline get_value(x, y, t, air_concentration::Number) = air_concentration @inline get_value(x, y, t, air_concentration::Function) = air_concentration(x, y, t) -@inline get_value(x, y, t, conc::Field{Center, Center, Center}) = interpolate(conc, Center(), Center(), Center(), conc.grid, x, y, 0) \ No newline at end of file +@inline get_value(x, y, t, conc::Field{Center, Center, Center}) = interpolate(conc, Center(), Center(), Center(), conc.grid, x, y, 0.0) \ No newline at end of file diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 573ce8316..d3ddf6443 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -25,7 +25,7 @@ function test_gas_exchange_model(grid, air_concentration) set!(model, T = 15.0, S = 35.0, DIC = 2220, Alk = 2500) # is everything communicating properly? (can't think of a way to not use allow scalar here) - value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1.0, 1.0, grid, model.clock, fields(model)) + value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1, 1, grid, model.clock, fields(model)) @test value ≈ -0.0003 atol = 0.0001 From 670076a01eedb258d3215e909a71de9eb3aaaffe Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 12:41:57 +0100 Subject: [PATCH 037/112] gas exhcnage boundary value fix --- src/Boundaries/gasexchange.jl | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index d8e4f3781..7af97c4f0 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -3,7 +3,7 @@ ##### As per OCMIP Phase 2 http://ocmip5.ipsl.jussieu.fr/OCMIP/phase2/simulations/Abiotic/Cchem/co2calc.f simplified as in PISCESv2 ##### -using Oceananigans.Fields: interpolate +using Oceananigans.Fields: fractional_indices, _interpolate struct pCO₂{P0, P1, P2, PB, PW, FT} solubility :: P0 @@ -259,4 +259,15 @@ end @inline get_value(x, y, t, air_concentration::Number) = air_concentration @inline get_value(x, y, t, air_concentration::Function) = air_concentration(x, y, t) -@inline get_value(x, y, t, conc::Field{Center, Center, Center}) = interpolate(conc, Center(), Center(), Center(), conc.grid, x, y, 0.0) \ No newline at end of file + +# interpolate doesn't really work on 2D fields +@inline function get_value(x, y, t, conc::Field{LX, LY, LZ}) where {LX, LY, LZ} + grid = conc.grid + + i, j, _ = fractional_indices(x, y, 0.0, (LX(), LY(), Center()), conc.grid) + + ξ, i = mod(i, 1), Base.unsafe_trunc(Int, i) + η, j = mod(j, 1), Base.unsafe_trunc(Int, j) + + return _interpolate(conc, ξ, η, 0, i+1, j+1, grid.Nz) +end \ No newline at end of file From 19f6e3bf7b4450afe20f429d2d14ad2e3f70fbde Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 12:42:16 +0100 Subject: [PATCH 038/112] corrected gas exchange test --- test/test_gasexchange.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index d3ddf6443..c5ce7af2e 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -54,11 +54,11 @@ end end @inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / year) +conc_field = CenterField(grid, indices=(:, :, grid.Nz)) +set!(conc_field, 413.0) @testset "Gas exchange coupling" begin grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) - conc_field = CenterField(grid, indices=(:, :, grid.Nz)) - conc_field .= 413.0 + 1.0 * rand() for air_concentration in [413.1, conc_function, conc_field] @info "Testing with $(typeof(air_concentration))" test_gas_exchange_model(grid, air_concentration) From c78d917b0065cbd34675ac720302581e08d5876c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 14:19:12 +0100 Subject: [PATCH 039/112] oops --- test/test_gasexchange.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index c5ce7af2e..bc0287fcf 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -53,12 +53,14 @@ end @test (mean(pCO₂_err) < 10 && std(pCO₂_err) < 15) end +grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) + @inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / year) + conc_field = CenterField(grid, indices=(:, :, grid.Nz)) set!(conc_field, 413.0) @testset "Gas exchange coupling" begin - grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) for air_concentration in [413.1, conc_function, conc_field] @info "Testing with $(typeof(air_concentration))" test_gas_exchange_model(grid, air_concentration) From d80cf05046a535347cb9cd7974490fd28355acda Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 14:45:04 +0100 Subject: [PATCH 040/112] fixed kelp test --- test/test_slatissima.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index 69f1462ee..661523e9a 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -17,12 +17,13 @@ sum_tracer_nitrogen(tracers) = sum(Array(interior(tracers.NO₃))) + sum(Array(interior(tracers.bPON))) + sum(Array(interior(tracers.DON))) -sum_tracer_carbon(tracers, redfield) = sum(Array(interior(tracers.sPOC))) + - sum(Array(interior(tracers.bPOC))) + - sum(Array(interior(tracers.DOC))) + - sum(Array(interior(tracers.DIC))) + - sum(Array(interior(tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) - + model.tracers.Z))) * redfield +sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = + sum(Array(interior(tracers.sPOC))) + + sum(Array(interior(tracers.bPOC))) + + sum(Array(interior(tracers.DOC))) + + sum(Array(interior(tracers.DIC))) + + sum(Array(interior(tracers.P * (1 + organic_carbon_calcate_ratio) + tracers.Z))) * redfield + @testset "SLatissima particle setup and conservations" begin grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) @@ -52,7 +53,7 @@ sum_tracer_carbon(tracers, redfield) = sum(Array(interior(tracers.sPOC))) + initial_tracer_N = sum_tracer_nitrogen(model.tracers) initial_kelp_N = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.N) .+ particles.structural_nitrogen)) ./ (14 * 0.001) - initial_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield) + initial_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield, model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) initial_kelp_C = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.C) .+ particles.structural_carbon)) ./ (12 * 0.001) model.clock.time = 60days # get to a high growth phase @@ -64,7 +65,7 @@ sum_tracer_carbon(tracers, redfield) = sum(Array(interior(tracers.sPOC))) + final_tracer_N = sum_tracer_nitrogen(model.tracers) final_kelp_N = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.N) .+ particles.structural_nitrogen)) ./ (14 * 0.001) - final_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield) + final_tracer_C = sum_tracer_carbon(model.tracers, model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield, model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) final_kelp_C = sum(Array(particles.A) .* particles.structural_dry_weight_per_area .* (Array(particles.C) .+ particles.structural_carbon)) ./ (12 * 0.001) # kelp is being integrated From 0f05b2198105c0764f8468a94c4246007a1dbbb3 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 14:47:21 +0100 Subject: [PATCH 041/112] fixed kelp test --- test/test_slatissima.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index 661523e9a..b68bb63b3 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -86,5 +86,5 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = run!(simulation) # the model is changing the tracer tendencies - not sure this test actually works as it didn't fail when it should have - @test any([any(intercepted_tendencies[tracer] .!= Array(interior(model.timestepper.Gⁿ[tracer]))) for tracer in (:NO₃, :NH₄, :DIC, :DOC, :bPON, :bPOC)]) + @test any([any(intercepted_tendencies[idx] .!= Array(interior(model.timestepper.Gⁿ[tracer]))) for (idx, tracer) in enumerate((:NO₃, :NH₄, :DIC, :DOC, :bPON, :bPOC))]) end From b33a2beb9043f6fb62e823afcb622ea1937fe65e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 14:50:21 +0100 Subject: [PATCH 042/112] fixed kelp test --- test/test_slatissima.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index b68bb63b3..01a66ab16 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -79,7 +79,7 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = simulation = Simulation(model, Δt = 1.0, stop_iteration = 1) # slow but easiest to have this just done on CPU - intercepted_tendencies = (Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) + intercepted_tendencies = Tuple(Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) From 72250523247836e182515cc50980f24fa2310b5c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sun, 17 Sep 2023 18:43:50 +0100 Subject: [PATCH 043/112] maybe finally fixed gas exchange --- src/Boundaries/gasexchange.jl | 11 ++++++++--- test/test_gasexchange.jl | 1 + 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index 7af97c4f0..632e0609b 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -3,7 +3,7 @@ ##### As per OCMIP Phase 2 http://ocmip5.ipsl.jussieu.fr/OCMIP/phase2/simulations/Abiotic/Cchem/co2calc.f simplified as in PISCESv2 ##### -using Oceananigans.Fields: fractional_indices, _interpolate +using Oceananigans.Fields: fractional_indices struct pCO₂{P0, P1, P2, PB, PW, FT} solubility :: P0 @@ -264,10 +264,15 @@ end @inline function get_value(x, y, t, conc::Field{LX, LY, LZ}) where {LX, LY, LZ} grid = conc.grid - i, j, _ = fractional_indices(x, y, 0.0, (LX(), LY(), Center()), conc.grid) + i, j, _ = fractional_indices(x, y, 0.0, (LX(), LY(), Center()), grid) ξ, i = mod(i, 1), Base.unsafe_trunc(Int, i) η, j = mod(j, 1), Base.unsafe_trunc(Int, j) - return _interpolate(conc, ξ, η, 0, i+1, j+1, grid.Nz) + _, _, ks = eachindex(conc).indices + + return @inbounds ((1 - ξ) * (1 - η) * conc[i, j, ks[1]] + + (1 - ξ) * η * conc[i, j+1, ks[1]] + + ξ * (1 - η) * conc[i+1, j, ks[1]] + + ξ * η * conc[i+1, j+1, ks[1]]) end \ No newline at end of file diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index bc0287fcf..18442b5e5 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -59,6 +59,7 @@ grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) conc_field = CenterField(grid, indices=(:, :, grid.Nz)) set!(conc_field, 413.0) +Oceananigans.fill_halo_regions!(conc_field) @testset "Gas exchange coupling" begin for air_concentration in [413.1, conc_function, conc_field] From 2d6616376e8658d99515264e9964e4220209b80b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 12:45:36 +0100 Subject: [PATCH 044/112] typos --- test/test_sediments.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index f988e6470..cf00d10f7 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -58,18 +58,18 @@ total_nitrogen(::LOBSTER, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.Z))) + sum(Array(interior(model.tracers.DOM))) + sum(Array(interior(model.tracers.sPOM))) + - sum(Array(interior(model.tracers.bPOM) + sum(Array(interior(model.tracers.bPOM))) total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.NH₄))) + sum(Array(interior(model.tracers.P))) + sum(Array(interior(model.tracers.Z))) + sum(Array(interior(model.tracers.DON))) + sum(Array(interior(model.tracers.sPON))) + - sum(Array(interior(model.tracers.bPON) + sum(Array(interior(model.tracers.bPON))) total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(Array(interior(model.tracers.N))) + sum(Array(interior(model.tracers.P))) + sum(Array(interior(model.tracers.Z))) + - sum(Array(interior(model.tracers.D) + sum(Array(interior(model.tracers.D))) function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAdamsBashforth2) model = isa(model, NonhydrostaticModel) ? model(; grid, From 689db1b001dff761e3e21e03eefea1f9dcf3606d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:24:02 +0100 Subject: [PATCH 045/112] more GPU stuff --- src/Boundaries/gasexchange.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index 632e0609b..ab43665f1 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -3,7 +3,7 @@ ##### As per OCMIP Phase 2 http://ocmip5.ipsl.jussieu.fr/OCMIP/phase2/simulations/Abiotic/Cchem/co2calc.f simplified as in PISCESv2 ##### -using Oceananigans.Fields: fractional_indices +using Oceananigans.Fields: fractional_indices, indices struct pCO₂{P0, P1, P2, PB, PW, FT} solubility :: P0 @@ -269,7 +269,7 @@ end ξ, i = mod(i, 1), Base.unsafe_trunc(Int, i) η, j = mod(j, 1), Base.unsafe_trunc(Int, j) - _, _, ks = eachindex(conc).indices + _, _, ks = indices(conc) return @inbounds ((1 - ξ) * (1 - η) * conc[i, j, ks[1]] + (1 - ξ) * η * conc[i, j+1, ks[1]] From 71978f2b96f02434c76f9d26044d2ff53ecc6d8e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:44:15 +0100 Subject: [PATCH 046/112] dropped support for using fields directly in gas exchange boudnaries --- src/Boundaries/gasexchange.jl | 7 +++++-- test/test_gasexchange.jl | 6 +++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index ab43665f1..dbb0ad06c 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -259,7 +259,9 @@ end @inline get_value(x, y, t, air_concentration::Number) = air_concentration @inline get_value(x, y, t, air_concentration::Function) = air_concentration(x, y, t) - +#=It is not possible for this to work on GPU since we need the grid and locs before (since fields are reduced to vectors) + We probably need a more generic and better way todo this but here is not the place. + For now if you wanted to use data you could do similar to https://github.com/OceanBioME/GlobalOceanBioME.jl/blob/main/src/one_degree_surface_par.jl # interpolate doesn't really work on 2D fields @inline function get_value(x, y, t, conc::Field{LX, LY, LZ}) where {LX, LY, LZ} grid = conc.grid @@ -275,4 +277,5 @@ end + (1 - ξ) * η * conc[i, j+1, ks[1]] + ξ * (1 - η) * conc[i+1, j, ks[1]] + ξ * η * conc[i+1, j+1, ks[1]]) -end \ No newline at end of file +end +=# \ No newline at end of file diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 18442b5e5..8018c5a2d 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -57,9 +57,9 @@ grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) @inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / year) -conc_field = CenterField(grid, indices=(:, :, grid.Nz)) -set!(conc_field, 413.0) -Oceananigans.fill_halo_regions!(conc_field) +#conc_field = CenterField(grid, indices=(:, :, grid.Nz)) +#set!(conc_field, 413.0) +#Oceananigans.fill_halo_regions!(conc_field) @testset "Gas exchange coupling" begin for air_concentration in [413.1, conc_function, conc_field] From b09592f779c8ab574082d8ad6a890d49bceab14e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:48:45 +0100 Subject: [PATCH 047/112] oops --- test/test_gasexchange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 8018c5a2d..08717ecce 100644 --- a/test/test_gasexchange.jl +++ b/test/test_gasexchange.jl @@ -62,7 +62,7 @@ grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) #Oceananigans.fill_halo_regions!(conc_field) @testset "Gas exchange coupling" begin - for air_concentration in [413.1, conc_function, conc_field] + for air_concentration in [413.1, conc_function]#, conc_field] @info "Testing with $(typeof(air_concentration))" test_gas_exchange_model(grid, air_concentration) end From 6b550b25ec65225bdf9b94c409d821a3a9d8739c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:52:24 +0100 Subject: [PATCH 048/112] possibly fixed kelp test --- test/test_slatissima.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index 01a66ab16..03019c411 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -22,7 +22,7 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = sum(Array(interior(tracers.bPOC))) + sum(Array(interior(tracers.DOC))) + sum(Array(interior(tracers.DIC))) + - sum(Array(interior(tracers.P * (1 + organic_carbon_calcate_ratio) + tracers.Z))) * redfield + sum(Array(interior(tracers.P)) .* (1 + organic_carbon_calcate_ratio) .+ Array(interior(tracers.Z))) .* redfield @testset "SLatissima particle setup and conservations" begin grid = RectilinearGrid(architecture; size=(1, 1, 1), extent=(1, 1, 1)) From a72e9c0ae0f0b5488ec7629d61d4bc78c80b5886 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:55:07 +0100 Subject: [PATCH 049/112] added adapt structures for kelp particles --- test/test_slatissima.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index 03019c411..c527e6d63 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -29,10 +29,11 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = # Initial properties - particles = SLatissima(; x = ones(Float64, 2), - A = ones(Float64, 2) .* 5, - N = ones(Float64, 2), - C = ones(Float64, 2), + particles = SLatissima(; architecture + x = arch_array(architecture, ones(Float64, 2)), + A = arch_array(architecture, ones(Float64, 2) .* 5), + N = arch_array(architecture, ones(Float64, 2)), + C = arch_array(architecture, ones(Float64, 2)), latitude = 1.0, pescribed_temperature = (args...) -> 10.0, pescribed_salinity = (args...) -> 35.0) From 3ff0ea932706b5b81f3092d1aa8f681e12596f26 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:55:45 +0100 Subject: [PATCH 050/112] typo --- test/test_slatissima.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index c527e6d63..bfdda15ad 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -29,7 +29,7 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = # Initial properties - particles = SLatissima(; architecture + particles = SLatissima(; architecture, x = arch_array(architecture, ones(Float64, 2)), A = arch_array(architecture, ones(Float64, 2) .* 5), N = arch_array(architecture, ones(Float64, 2)), From f5ccfd1837b00b77320d0a862bb60d0188feb573 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 13:56:20 +0100 Subject: [PATCH 051/112] forgot to add arch array --- test/test_slatissima.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index bfdda15ad..c9547ae43 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -2,6 +2,7 @@ using Test, OceanBioME, Oceananigans using OceanBioME.SLatissimaModel: SLatissima using Oceananigans.Units using Oceananigans.Fields: TracerFields +using Oceananigans.Architectures: arch_array function intercept_tendencies!(model, intercepted_tendencies) for tracer in keys(model.tracers) From 9f93ebe8724375e2348ac38d21b5c887a48b8fe4 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 14:16:35 +0100 Subject: [PATCH 052/112] corrected kelp test tollerance for GPU --- test/test_slatissima.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index c9547ae43..93bd2cca7 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -75,8 +75,12 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = @test initial_kelp_C != final_kelp_C # conservaitons - @test initial_tracer_N + initial_kelp_N ≈ final_tracer_N + final_kelp_N - @test initial_tracer_C + initial_kelp_C ≈ final_tracer_C + final_kelp_C + # (GPU eps is much larger (~10⁻⁷) than on CPU) + rtol = ifelse(isa(architecture, CPU), max(eps(initial_tracer_N + initial_kelp_N), eps(final_tracer_N + final_kelp_N)), 2e-7) + @test isapprox(initial_tracer_N + initial_kelp_N, final_tracer_N + final_kelp_N; rtol) + + rtol = ifelse(isa(architecture, CPU), max(eps(initial_tracer_C + initial_kelp_C), eps(final_tracer_C + final_kelp_C)), 7e-7) + @test isapprox(initial_tracer_C + initial_kelp_C, final_tracer_C + final_kelp_C; rtol) simulation = Simulation(model, Δt = 1.0, stop_iteration = 1) From d3485cc257104d51293a78ffd2e6910699b7feed Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 14:20:21 +0100 Subject: [PATCH 053/112] arch array for sediments --- src/Boundaries/Sediments/instant_remineralization.jl | 2 +- src/Boundaries/Sediments/simple_multi_G.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index f053fac39..7ef1d0d05 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -63,7 +63,7 @@ function InstantRemineralisation(; grid, 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))) - bottom_indices = calculate_bottom_indices(grid) + bottom_indices = arch_array(architecture(grid), calculate_bottom_indices(grid)) return InstantRemineralisation(burial_efficiency_constant1, burial_efficiency_constant2, diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index f0919ab3d..57ac8e9ff 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -113,7 +113,7 @@ function SimpleMultiG(; grid, 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))) - bottom_indices = calculate_bottom_indices(grid) + bottom_indices = arch_array(architecture(grid), calculate_bottom_indices(grid)) F = typeof(fields) TE = typeof(tendencies) @@ -144,7 +144,8 @@ adapt_structure(to, sediment::SimpleMultiG) = sediment.anoxic_params, sediment.solid_dep_params, adapt(to, sediment.fields), - adapt(to, sediment.tendencies)) + adapt(to, sediment.tendencies), + adapt(to, sediment.bottom_indices)) sediment_tracers(::SimpleMultiG) = (:C_slow, :C_fast, :C_ref, :N_slow, :N_fast, :N_ref) sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, From 043a2bebddffa4d7c5fb22138b4d9c310a785179 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:05:52 +0100 Subject: [PATCH 054/112] attempt to fix sediment tests --- test/test_sediments.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index cf00d10f7..02ead51bf 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -13,7 +13,7 @@ using OceanBioME.LOBSTERModel: lobster_variable_redfield function intercept_tendencies!(model, intercepted_tendencies) for tracer in keys(model.tracers) - copyto!(intercepted_tendencies[tracer], model.timestepper.Gⁿ[tracer]) + intercepted_tendencies[tracer] = Array(interior(model.timestepper.Gⁿ[tracer])) end end @@ -89,7 +89,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd simulation = Simulation(model, Δt = 50, stop_time = 1day) - intercepted_tendencies = TracerFields(keys(model.tracers), grid) + intercepted_tendencies = Tuple(Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) @@ -98,11 +98,11 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd run!(simulation) # the model is changing the tracer tendencies - @test any([any(intercepted_tendencies[tracer] .!= model.timestepper.Gⁿ[tracer]) for tracer in keys(model.tracers)]) + @test any([any(intercepted_tendencies[idx] .!= Array(interior(model.timestepper.Gⁿ[tracer]))) for (idx, tracer) in enumerate(keys(model.tracers))]) # the sediment tendencies are being updated - @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.Gⁿ]) - @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.G⁻]) + @test all([any(Array(interior(tend)) .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.Gⁿ]) + @test all([any(Array(interior(tend)) .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.G⁻]) # the sediment values are being integrated initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) @@ -111,7 +111,8 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd N₁ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) # conservations - @test N₁ ≈ N₀ + rtol = ifelse(isa(architecture, CPU), max(eps(N₀), eps(N₁)), 5e-7) + @test isapprox(N₀, N₁; rtol) return nothing end From 90d2ddeb8b81decff2ba027c5fde2ce172c70211 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:08:23 +0100 Subject: [PATCH 055/112] this shouldn't have worked anywhere before --- src/Boundaries/Sediments/Sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 8f80f958f..120ab9f3d 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -37,7 +37,7 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - sediment, bgc.underlying_biogeochemistry, model.grid, model.advection, model.tracers, model.timestepper) + sediment, bgc.underlying_biogeochemistry, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ) return nothing end From d2e8cc29a009649ecc718da27a47146aa346827e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:09:42 +0100 Subject: [PATCH 056/112] sediment test fix --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 02ead51bf..2d97e057c 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -106,7 +106,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd # the sediment values are being integrated initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) - @test all([any(field .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) + @test all([any(Array(interior(field)) .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) N₁ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) From 3a95540ed336c3c831a6aee144fda16f1f6e5231 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:11:19 +0100 Subject: [PATCH 057/112] sediment test feedback fix --- test/test_sediments.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 2d97e057c..249d0529c 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -149,8 +149,8 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) if run - @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry, model; timestepper) + @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))" + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" test_flat_sediment(grid, biogeochemistry, model; timestepper) end end end From 185e7b1d947986241a7cd372f6dea40188827534 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:15:33 +0100 Subject: [PATCH 058/112] sediment adapt fix --- src/Boundaries/Sediments/simple_multi_G.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 57ac8e9ff..178c1f348 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -25,6 +25,25 @@ struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B} <: FlatSediment fields :: F tendencies :: TE bottom_indices :: B + + SimpleMultiG(fast_decay_rate::FT, slow_decay_rate::FT, + fast_redfield::FT, slow_redfield::FT, + fast_fraction::FT, slow_fraction::FT, refactory_fraction::FT, + nitrate_oxidation_params::P1, + denitrification_params::P2, + anoxic_params::P3, + solid_dep_params::P4, + fields::F, tendencies::TE, + bottom_indices::B) where {FT, P1, P2, P3, P4, F, TE, B} = + new{FT, P1, P2, P3, P4, F, TE, B}(fast_decay_rate, slow_decay_rate, + fast_redfield, slow_redfield, + fast_fraction, slow_fraction, refactory_fraction, + nitrate_oxidation_params, + denitrification_params, + anoxic_params, + solid_dep_params, + fields, tendencies, + bottom_indices) end """ From 02eb7b56f9500fd4684da7a637bd8900dc7db4c6 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:17:52 +0100 Subject: [PATCH 059/112] another multi g fix --- src/Boundaries/Sediments/simple_multi_G.jl | 70 ++++++++++------------ 1 file changed, 33 insertions(+), 37 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 178c1f348..b9a75a29d 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -48,18 +48,18 @@ end """ SimpleMultiG(; grid - fast_decay_rate::FT = 2/day, - slow_decay_rate::FT = 0.2/day, - fast_redfield::FT = 0.1509, - slow_redfield::FT = 0.13, - fast_fraction::FT = 0.74, - slow_fraction::FT = 0.26, - refactory_fraction::FT = 0.1, - nitrate_oxidation_params::P1 = (A = - 1.9785, B = 0.2261, C = -0.0615, D = -0.0289, E = - 0.36109, F = - 0.0232), - denitrification_params::P2 = (A = - 3.0790, B = 1.7509, C = 0.0593, D = - 0.1923, E = 0.0604, F = 0.0662), - anoxic_params::P3 = (A = - 3.9476, B = 2.6269, C = - 0.2426, D = -1.3349, E = 0.1826, F = - 0.0143), + fast_decay_rate = 2/day, + slow_decay_rate = 0.2/day, + fast_redfield = 0.1509, + slow_redfield = 0.13, + fast_fraction = 0.74, + slow_fraction = 0.26, + refactory_fraction = 0.1, + nitrate_oxidation_params = (A = - 1.9785, B = 0.2261, C = -0.0615, D = -0.0289, E = - 0.36109, F = - 0.0232), + denitrification_params = (A = - 3.0790, B = 1.7509, C = 0.0593, D = - 0.1923, E = 0.0604, F = 0.0662), + anoxic_params = (A = - 3.9476, B = 2.6269, C = - 0.2426, D = -1.3349, E = 0.1826, F = - 0.0143), depth = abs(znodes(grid, Face())[1]), - solid_dep_params::P4 = (A = 0.233, B = 0.336, C = 982, D = - 1.548, depth = depth)) + solid_dep_params = (A = 0.233, B = 0.336, C = 982, D = - 1.548, depth = depth)) Return a single-layer "multi G" sediment model (`SimpleMultiG`) on `grid`, where parameters can be optionally specified. @@ -90,37 +90,37 @@ Single-layer multi-G sediment model (Float64) ``` """ function SimpleMultiG(; grid, - fast_decay_rate::FT = 2/day, - slow_decay_rate::FT = 0.2/day, - fast_redfield::FT = 0.1509, - slow_redfield::FT = 0.13, - fast_fraction::FT = 0.74, - slow_fraction::FT = 0.26, - refactory_fraction::FT = 0.1, - nitrate_oxidation_params::P1 = (A = - 1.9785, + fast_decay_rate = 2/day, + slow_decay_rate = 0.2/day, + fast_redfield = 0.1509, + slow_redfield = 0.13, + fast_fraction = 0.74, + slow_fraction = 0.26, + refactory_fraction = 0.1, + nitrate_oxidation_params = (A = - 1.9785, B = 0.2261, C = -0.0615, D = -0.0289, E = - 0.36109, F = - 0.0232), - denitrification_params::P2 = (A = - 3.0790, + denitrification_params = (A = - 3.0790, B = 1.7509, C = 0.0593, D = - 0.1923, E = 0.0604, F = 0.0662), - anoxic_params::P3 = (A = - 3.9476, + anoxic_params = (A = - 3.9476, B = 2.6269, C = - 0.2426, D = -1.3349, E = 0.1826, F = - 0.0143), depth = abs(znodes(grid, Face())[1]), - solid_dep_params::P4 = (A = 0.233, + solid_dep_params = (A = 0.233, B = 0.336, C = 982, D = - 1.548, - depth = depth)) where {FT, P1, P2, P3, P4} + depth = depth)) @warn "Sediment models are an experimental feature and have not yet been validated." @info "This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC." @@ -134,20 +134,16 @@ function SimpleMultiG(; grid, bottom_indices = arch_array(architecture(grid), calculate_bottom_indices(grid)) - F = typeof(fields) - TE = typeof(tendencies) - B = typeof(bottom_indices) - - return SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B}(fast_decay_rate, slow_decay_rate, - fast_redfield, slow_redfield, - fast_fraction, slow_fraction, refactory_fraction, - nitrate_oxidation_params, - denitrification_params, - anoxic_params, - solid_dep_params, - fields, - tendencies, - bottom_indices) + return SimpleMultiG(fast_decay_rate, slow_decay_rate, + fast_redfield, slow_redfield, + fast_fraction, slow_fraction, refactory_fraction, + nitrate_oxidation_params, + denitrification_params, + anoxic_params, + solid_dep_params, + fields, + tendencies, + bottom_indices) end adapt_structure(to, sediment::SimpleMultiG) = From f107df0a9fd4aac7a66597db5e74d0ad19f90b46 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 15:23:24 +0100 Subject: [PATCH 060/112] immersed boundaries fix for multi g --- src/Boundaries/Sediments/simple_multi_G.jl | 40 ++++++++++------------ 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index b9a75a29d..26242c875 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -186,30 +186,30 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, nitrogen_deposition = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz # rates - C_min_slow = sediment.fields.C_slow[i, j, 1] * sediment.slow_decay_rate - C_min_fast = sediment.fields.C_fast[i, j, 1] * sediment.fast_decay_rate + C_min_slow = sediment.fields.C_slow[i, j, k] * sediment.slow_decay_rate + C_min_fast = sediment.fields.C_fast[i, j, k] * sediment.fast_decay_rate - N_min_slow = sediment.fields.N_slow[i, j, 1] * sediment.slow_decay_rate - N_min_fast = sediment.fields.N_fast[i, j, 1] * sediment.fast_decay_rate + N_min_slow = sediment.fields.N_slow[i, j, k] * sediment.slow_decay_rate + N_min_fast = sediment.fields.N_fast[i, j, k] * sediment.fast_decay_rate Cᵐⁱⁿ = C_min_slow + C_min_fast Nᵐⁱⁿ = N_min_slow + N_min_fast - k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1]) + k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, k] + sediment.fields.C_fast[i, j, k]) # sediment evolution - sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow - sediment.tendencies.Gⁿ.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast - sediment.tendencies.Gⁿ.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition + sediment.tendencies.Gⁿ.C_slow[i, j, k] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow + sediment.tendencies.Gⁿ.C_fast[i, j, k] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast + sediment.tendencies.Gⁿ.C_ref[i, j, k] = sediment.refactory_fraction * carbon_deposition - sediment.tendencies.Gⁿ.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow - sediment.tendencies.Gⁿ.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast - sediment.tendencies.Gⁿ.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition + sediment.tendencies.Gⁿ.N_slow[i, j, k] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow + sediment.tendencies.Gⁿ.N_fast[i, j, k] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast + sediment.tendencies.Gⁿ.N_ref[i, j, k] = sediment.refactory_fraction * nitrogen_deposition # efflux/influx - O₂ = tracers.O₂[i, j, 1] - NO₃ = tracers.NO₃[i, j, 1] - NH₄ = tracers.NH₄[i, j, 1] + O₂ = tracers.O₂[i, j, k] + NO₃ = tracers.NO₃[i, j, k] + NH₄ = tracers.NH₄[i, j, k] pₙᵢₜ = exp(sediment.nitrate_oxidation_params.A + sediment.nitrate_oxidation_params.B * log(Cᵐⁱⁿ * day) * log(O₂) + @@ -234,18 +234,16 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, sediment.anoxic_params.F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) if isnan(pₐₙₒₓ) - println("$(Cᵐⁱⁿ), $(k), $(O₂), $(NO₃)") + #println("$(Cᵐⁱⁿ), $(k), $(O₂), $(NO₃)") error("Sediment anoxia has caused model failure") end pₛₒₗᵢ = sediment.solid_dep_params.A * (sediment.solid_dep_params.C * sediment.solid_dep_params.depth ^ sediment.solid_dep_params.D) ^ sediment.solid_dep_params.B - Δz = grid.Δzᵃᵃᶜ[1] - - tendencies.NH₄[i, j, 1] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz - tendencies.NO₃[i, j, 1] += Nᵐⁱⁿ * pₙᵢₜ / Δz - tendencies.DIC[i, j, 1] += Cᵐⁱⁿ / Δz - tendencies.O₂[i, j, 1] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly + tendencies.NH₄[i, j, k] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz + tendencies.NO₃[i, j, k] += Nᵐⁱⁿ * pₙᵢₜ / Δz + tendencies.DIC[i, j, k] += Cᵐⁱⁿ / Δz + tendencies.O₂[i, j, k] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly end end From 97253fd2d3f0f49a3848a0b665910659e890d9a7 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 16:17:24 +0100 Subject: [PATCH 061/112] some sediment bugs --- .../Sediments/instant_remineralization.jl | 2 +- src/Boundaries/Sediments/simple_multi_G.jl | 25 +++++++++---------- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 7ef1d0d05..e4c76337b 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -100,7 +100,7 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st # sediment evolution @inbounds sediment.tendencies.Gⁿ.N_storage[i, j, 1] = burial_efficiency * flux - @inbounds tendencies[remineralisation_receiver(bgc)][i, j, 1] += flux * (1 - burial_efficiency) / Δz + @inbounds tendencies[remineralisation_receiver(bgc)][i, j, k] += flux * (1 - burial_efficiency) / Δz end summary(::InstantRemineralisation{FT}) where {FT} = string("Single-layer instant remineralisaiton ($FT)") diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 26242c875..1918cea33 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -180,31 +180,30 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, Δz = zspacing(i, j, k, grid, Center(), Center(), Center()) @inbounds begin - carbon_deposition = carbon_flux(i, j, k, grid, advection, bgc, tracers) * Δz nitrogen_deposition = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz # rates - C_min_slow = sediment.fields.C_slow[i, j, k] * sediment.slow_decay_rate - C_min_fast = sediment.fields.C_fast[i, j, k] * sediment.fast_decay_rate + C_min_slow = sediment.fields.C_slow[i, j, 1] * sediment.slow_decay_rate + C_min_fast = sediment.fields.C_fast[i, j, 1] * sediment.fast_decay_rate - N_min_slow = sediment.fields.N_slow[i, j, k] * sediment.slow_decay_rate - N_min_fast = sediment.fields.N_fast[i, j, k] * sediment.fast_decay_rate + N_min_slow = sediment.fields.N_slow[i, j, 1] * sediment.slow_decay_rate + N_min_fast = sediment.fields.N_fast[i, j, 1] * sediment.fast_decay_rate Cᵐⁱⁿ = C_min_slow + C_min_fast Nᵐⁱⁿ = N_min_slow + N_min_fast - k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, k] + sediment.fields.C_fast[i, j, k]) + k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1]) # sediment evolution - sediment.tendencies.Gⁿ.C_slow[i, j, k] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow - sediment.tendencies.Gⁿ.C_fast[i, j, k] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast - sediment.tendencies.Gⁿ.C_ref[i, j, k] = sediment.refactory_fraction * carbon_deposition + sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow + sediment.tendencies.Gⁿ.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast + sediment.tendencies.Gⁿ.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition - sediment.tendencies.Gⁿ.N_slow[i, j, k] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow - sediment.tendencies.Gⁿ.N_fast[i, j, k] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast - sediment.tendencies.Gⁿ.N_ref[i, j, k] = sediment.refactory_fraction * nitrogen_deposition + sediment.tendencies.Gⁿ.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow + sediment.tendencies.Gⁿ.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast + sediment.tendencies.Gⁿ.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition # efflux/influx O₂ = tracers.O₂[i, j, k] @@ -243,7 +242,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, tendencies.NH₄[i, j, k] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz tendencies.NO₃[i, j, k] += Nᵐⁱⁿ * pₙᵢₜ / Δz tendencies.DIC[i, j, k] += Cᵐⁱⁿ / Δz - tendencies.O₂[i, j, k] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly + tendencies.O₂[i, j, k] -= max(0, ((1 - pₐₙₒₓ * pₛₒₗᵢ) * Cᵐⁱⁿ + 2 * Nᵐⁱⁿ * pₙᵢₜ)/ Δz) # this seems dodge but this model doesn't cope with anoxia properly (I think) end end From a7ce62fba8be196b28f20d9f83ebf1d07c5205fb Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 16:36:53 +0100 Subject: [PATCH 062/112] spotted the bug finally --- src/Boundaries/Sediments/simple_multi_G.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 1918cea33..4fe0320a2 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -194,7 +194,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, Cᵐⁱⁿ = C_min_slow + C_min_fast Nᵐⁱⁿ = N_min_slow + N_min_fast - k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1]) + reactivity = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1]) # sediment evolution sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow @@ -213,7 +213,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, pₙᵢₜ = exp(sediment.nitrate_oxidation_params.A + sediment.nitrate_oxidation_params.B * log(Cᵐⁱⁿ * day) * log(O₂) + sediment.nitrate_oxidation_params.C * log(Cᵐⁱⁿ * day) ^ 2 + - sediment.nitrate_oxidation_params.D * log(k) * log(NH₄) + + sediment.nitrate_oxidation_params.D * log(reactivity) * log(NH₄) + sediment.nitrate_oxidation_params.E * log(Cᵐⁱⁿ * day) + sediment.nitrate_oxidation_params.F * log(Cᵐⁱⁿ * day) * log(NH₄)) / (Nᵐⁱⁿ * day) @@ -222,18 +222,18 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, sediment.denitrification_params.B * log(Cᵐⁱⁿ * day) + sediment.denitrification_params.C * log(NO₃) ^ 2 + sediment.denitrification_params.D * log(Cᵐⁱⁿ * day) ^ 2 + - sediment.denitrification_params.E * log(k) ^ 2 + - sediment.denitrification_params.F * log(O₂) * log(k)) / (Cᵐⁱⁿ * day) + sediment.denitrification_params.E * log(reactivity) ^ 2 + + sediment.denitrification_params.F * log(O₂) * log(reactivity)) / (Cᵐⁱⁿ * day) =# pₐₙₒₓ = exp(sediment.anoxic_params.A + sediment.anoxic_params.B * log(Cᵐⁱⁿ * day) + sediment.anoxic_params.C * log(Cᵐⁱⁿ * day) ^ 2 + - sediment.anoxic_params.D * log(k) + - sediment.anoxic_params.E * log(O₂) * log(k) + + sediment.anoxic_params.D * log(reactivity) + + sediment.anoxic_params.E * log(O₂) * log(reactivity) + sediment.anoxic_params.F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) if isnan(pₐₙₒₓ) - #println("$(Cᵐⁱⁿ), $(k), $(O₂), $(NO₃)") + println("$(Cᵐⁱⁿ), $(reactivity), $(O₂), $(NO₃)") error("Sediment anoxia has caused model failure") end From cfb2175a6ee216a9d35ec44c3f76b4a116c7c4bb Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 16:56:09 +0100 Subject: [PATCH 063/112] multi g sediment will now compile but parameter space too big --- src/Boundaries/Sediments/simple_multi_G.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 4fe0320a2..3dc8a96a8 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -233,7 +233,6 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, sediment.anoxic_params.F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) if isnan(pₐₙₒₓ) - println("$(Cᵐⁱⁿ), $(reactivity), $(O₂), $(NO₃)") error("Sediment anoxia has caused model failure") end From cf328849c722a25c8c4dc63ac04864569264dfda Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 16:59:35 +0100 Subject: [PATCH 064/112] marked multi G on GPU as broken --- test/test_sediments.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 249d0529c..1fb9d0532 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -148,6 +148,14 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) + + if isa(sediment_model, SimpleMultiG) && isa(architecture, GPU) + run = false + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin + @test false broken = true + end + end + if run @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))" @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" test_flat_sediment(grid, biogeochemistry, model; timestepper) From df04f8e110635e44bf2433e40bcb0ec7a72e3190 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 17:26:22 +0100 Subject: [PATCH 065/112] bump minor version This PR breaks the `ScaleNegativeTracers` API --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4761ddc56..9ab81d17c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright ", "John R Taylor ", "Si Chen "] -version = "0.6.0" +version = "0.7.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From ec474eab5f7e0d0d2a94d16c545bd8cff106b72c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 18:45:12 +0100 Subject: [PATCH 066/112] rtols were incorrect --- test/test_sediments.jl | 2 +- test/test_slatissima.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 1fb9d0532..7cf1e902e 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -111,7 +111,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd N₁ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) # conservations - rtol = ifelse(isa(architecture, CPU), max(eps(N₀), eps(N₁)), 5e-7) + rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) @test isapprox(N₀, N₁; rtol) return nothing diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index 93bd2cca7..9cd4bdc48 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -76,10 +76,10 @@ sum_tracer_carbon(tracers, redfield, organic_carbon_calcate_ratio) = # conservaitons # (GPU eps is much larger (~10⁻⁷) than on CPU) - rtol = ifelse(isa(architecture, CPU), max(eps(initial_tracer_N + initial_kelp_N), eps(final_tracer_N + final_kelp_N)), 2e-7) + rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_N + initial_kelp_N), √eps(final_tracer_N + final_kelp_N)), 2e-7) @test isapprox(initial_tracer_N + initial_kelp_N, final_tracer_N + final_kelp_N; rtol) - rtol = ifelse(isa(architecture, CPU), max(eps(initial_tracer_C + initial_kelp_C), eps(final_tracer_C + final_kelp_C)), 7e-7) + rtol = ifelse(isa(architecture, CPU), max(√eps(initial_tracer_C + initial_kelp_C), √eps(final_tracer_C + final_kelp_C)), 7e-7) @test isapprox(initial_tracer_C + initial_kelp_C, final_tracer_C + final_kelp_C; rtol) simulation = Simulation(model, Δt = 1.0, stop_iteration = 1) From 54404a993447730c526bae6242f727d14437fbde Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 18:49:57 +0100 Subject: [PATCH 067/112] Documenter 1.0 just dropped and it breaks a lot of stuff that I will fix in a different PR --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 9ab81d17c..cbb241cdf 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ Oceananigans = "0.84.1, 0.85, 0.86, 0.87" Roots = "2" StructArrays = "0.4, 0.5, 0.6" julia = "1.9" +Documenter = "0.27" [extras] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" From 8f133dfe5301e9241eddb3a5e5a74b0c0559142e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 19:25:52 +0100 Subject: [PATCH 068/112] add method for box model grid --- src/BoxModel/boxmodel.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/BoxModel/boxmodel.jl b/src/BoxModel/boxmodel.jl index fdac2c50e..935e944db 100644 --- a/src/BoxModel/boxmodel.jl +++ b/src/BoxModel/boxmodel.jl @@ -10,7 +10,7 @@ using Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields -using Oceananigans: Clock, prettytime +using Oceananigans: Clock, prettytime, CPU using Oceananigans.TimeSteppers: tick! using OceanBioME: BoxModelGrid @@ -19,6 +19,7 @@ using StructArrays, JLD2 import Oceananigans.Simulations: run! import Oceananigans: set! import Oceananigans.Fields: CenterField, regularize_field_boundary_conditions +import Oceananigans.Architectures: architecture import Base: show, summary @inline no_func(args...) = 0.0 @@ -188,5 +189,6 @@ show(io::IO, model::BoxModel{B, V, FT, F, TS, C}) where {B, V, FT, F, TS, C} = CenterField(::BoxModelGrid, args...; kwargs...) = nothing regularize_field_boundary_conditions(::Any, ::BoxModelGrid, ::Any) = nothing +architecture(::BoxModelGrid) = CPU() end # module From 30cdfba55bf60dce235a7794996b4d00a6609a06 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 21:37:34 +0100 Subject: [PATCH 069/112] tried to simplify multi G adapted to GPU [skip ci] --- src/Boundaries/Sediments/simple_multi_G.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 3dc8a96a8..404c73a11 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -159,7 +159,7 @@ adapt_structure(to, sediment::SimpleMultiG) = sediment.anoxic_params, sediment.solid_dep_params, adapt(to, sediment.fields), - adapt(to, sediment.tendencies), + adapt(to, sediment.tendencies.Gⁿ), adapt(to, sediment.bottom_indices)) sediment_tracers(::SimpleMultiG) = (:C_slow, :C_fast, :C_ref, :N_slow, :N_fast, :N_ref) @@ -197,13 +197,13 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, reactivity = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1]) # sediment evolution - sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow - sediment.tendencies.Gⁿ.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast - sediment.tendencies.Gⁿ.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition + sediment.tendencies.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow + sediment.tendencies.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast + sediment.tendencies.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition - sediment.tendencies.Gⁿ.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow - sediment.tendencies.Gⁿ.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast - sediment.tendencies.Gⁿ.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition + sediment.tendencies.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow + sediment.tendencies.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast + sediment.tendencies.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition # efflux/influx O₂ = tracers.O₂[i, j, k] From 215e0a2712c63e2f75f66521780b0aa29fba6864 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 21:43:47 +0100 Subject: [PATCH 070/112] removed errors from kernels --- src/Boundaries/Sediments/simple_multi_G.jl | 4 ---- src/Utils/negative_tracers.jl | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 404c73a11..15c278827 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -232,10 +232,6 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, sediment.anoxic_params.E * log(O₂) * log(reactivity) + sediment.anoxic_params.F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) - if isnan(pₐₙₒₓ) - error("Sediment anoxia has caused model failure") - end - pₛₒₗᵢ = sediment.solid_dep_params.A * (sediment.solid_dep_params.C * sediment.solid_dep_params.depth ^ sediment.solid_dep_params.D) ^ sediment.solid_dep_params.B tendencies.NH₄[i, j, k] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 0e8b8a003..ad9940cc1 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -111,7 +111,7 @@ end end end - t < 0 && error("Cell total < 0, can not scale negative tracers.") + t < 0 && t = NaN @unroll for tracer in tracers value = @inbounds tracer[i, j, k] From 19a40f189b95bad505da500ebf4fb2697a61fdb1 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 22:06:47 +0100 Subject: [PATCH 071/112] removed column diffusion stuff --- src/Utils/negative_tracers.jl | 2 +- src/Utils/timestep.jl | 16 ---------------- test/test_utils.jl | 18 ------------------ 3 files changed, 1 insertion(+), 35 deletions(-) diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index ad9940cc1..8cc839ac9 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -111,7 +111,7 @@ end end end - t < 0 && t = NaN + t < 0 && (t = NaN) @unroll for tracer in tracers value = @inbounds tracer[i, j, k] diff --git a/src/Utils/timestep.jl b/src/Utils/timestep.jl index 910b88064..f6aa7db6e 100644 --- a/src/Utils/timestep.jl +++ b/src/Utils/timestep.jl @@ -1,22 +1,6 @@ using Oceananigans.Advection: cell_advection_timescale using Oceananigans.Grids: Center, znodes, zspacing, minimum_zspacing -# NOT GPU SAFE -@inline function column_diffusion_timescale(model) - grid = model.grid - - z = znodes(grid, Center(), Center(), Center()) - - Δz2_ν = zeros(grid.Nz) - - for k in 1:grid.Nz - Δzₖ = zspacing(1, 1, k, grid, Center(), Center(), Center()) - @inbounds Δz2_ν[k] = Δzₖ^2 / model.closure.κ[1](0.0, 0.0, z[k], model.clock.time) # assumes all tracer closures are the same and x/y invariant - end - - return minimum(Δz2_ν) -end - @inline column_advection_timescale(model) = minimum_zspacing(model.grid) / maximum_sinking_velocity(model.biogeochemistry) @inline sinking_advection_timescale(model) = min(minimum_zspacing(model.grid) / maximum_sinking_velocity(model.biogeochemistry), cell_advection_timescale(model)) diff --git a/test/test_utils.jl b/test/test_utils.jl index 8f33d8bc9..780e99030 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -1,22 +1,5 @@ using Test, OceanBioME, Oceananigans -function test_column_diffusion_timescale(arch) - κ = 1e-3 - @inline κₜ(x, y, z, t) = κ - - grid = RectilinearGrid(arch, size=(1, 1, 5), x=(0, 10), y=(0, 3), z=[-1, -0.6, -0.5, -0.2, -0.19, 0]) - min_Δz = minimum_zspacing(grid) - - model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), - biogeochemistry = LOBSTER(; grid, - surface_phytosynthetically_active_radiation = (x, y, t) -> 100, - carbonates = true), - advection = nothing) - - return column_diffusion_timescale(model) ≈ min_Δz ^ 2 / κ -end - function test_negative_scaling(arch) grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) @@ -53,7 +36,6 @@ function test_negative_zeroing(arch) end @testset "Test Utils" begin - @test test_column_diffusion_timescale(architecture) broken = architecture == GPU() @test test_negative_scaling(architecture) @test test_negative_zeroing(architecture) end From fb35422a548578f66c956886d010a8e7d6515e1d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 22:10:04 +0100 Subject: [PATCH 072/112] missed an allow scalar --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 7cf1e902e..5300d4966 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -93,7 +93,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) - N₀ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + N₀ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) run!(simulation) From f10ff1084ab92f45fd86e92490103853016426e9 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 22:12:39 +0100 Subject: [PATCH 073/112] fixed docs compats --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index a85e13173..be993e459 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,3 +18,4 @@ CairoMakie = "0.10" DocumenterCitations = "1" Literate = "≥2.9.0" Oceananigans = "0.84.1, 0.85, 0.86, 0.87" +Documenter = "0.27" From 41188736162c0d069d5bff82ba6792e91040d96b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 22:16:47 +0100 Subject: [PATCH 074/112] removed broken flag for multi g test --- test/test_sediments.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 5300d4966..9d014fc6a 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -148,13 +148,6 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) - - if isa(sediment_model, SimpleMultiG) && isa(architecture, GPU) - run = false - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin - @test false broken = true - end - end if run @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))" From 61e10f40d8a1c87b9d52047285aa8d647696200b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 22:27:07 +0100 Subject: [PATCH 075/112] missed an allow scalar --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 9d014fc6a..6bfb76aa5 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -108,7 +108,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) @test all([any(Array(interior(field)) .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) - N₁ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + N₁ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) # conservations rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) From 31cb1b7de162a1e5e750c92ba17ca542c91c5ab9 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 23:00:29 +0100 Subject: [PATCH 076/112] use simpler LOBSTER model with Multi G to fit parameter space --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 6bfb76aa5..4f00aee91 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -142,7 +142,7 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill 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), + variable_redfield = ifelse(isa(sediment_model, SimpleMultiG), false, true), sediment_model)) # get rid of incompatible combinations run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || From f5d6f4a6e539692f77832df31d8ba9d2aa57b7b6 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 23:12:45 +0100 Subject: [PATCH 077/112] fix LOBSTER sediment test [skip ci] --- test/test_sediments.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 4f00aee91..9a34a1e7d 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -34,7 +34,7 @@ set_defaults!(::LOBSTER, model) = sPOM = 0.2299, bPOM = 0.0103) -set_defaults!(::lobster_variable_redfield, model) = +set_defaults!(::LOBSTER{<:Any, Val{(true, true, true)}}, model) = set!(model, P = 0.4686, Z = 0.5363, NO₃ = 2.3103, NH₄ = 0.0010, DIC = 2106.9, Alk = 2408.9, @@ -43,6 +43,13 @@ set_defaults!(::lobster_variable_redfield, model) = sPON = 0.2299, sPOC = 1.5080, bPON = 0.0103, bPOC = 0.0781) +set_defaults!(::LOBSTER{<:Any, Val{(false, false, true)}}, model) = + set!(model, P = 0.4686, Z = 0.5363, + NO₃ = 2.3103, NH₄ = 0.0010, + DOC = 5.3390, DON = 0.8115, + sPON = 0.2299, sPOC = 1.5080, + bPON = 0.0103, bPOC = 0.0781) + set_defaults!(::NutrientPhytoplanktonZooplanktonDetritus, model) = set!(model, N = 2.3, P = 0.4, Z = 0.5, D = 0.2) @@ -59,6 +66,7 @@ total_nitrogen(::LOBSTER, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.DOM))) + sum(Array(interior(model.tracers.sPOM))) + sum(Array(interior(model.tracers.bPOM))) + total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.NH₄))) + sum(Array(interior(model.tracers.P))) + @@ -66,6 +74,7 @@ total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tr sum(Array(interior(model.tracers.DON))) + sum(Array(interior(model.tracers.sPON))) + sum(Array(interior(model.tracers.bPON))) + total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(Array(interior(model.tracers.N))) + sum(Array(interior(model.tracers.P))) + sum(Array(interior(model.tracers.Z))) + From ebea1e06abe1504930a41d18370340951ad19aa3 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 18 Sep 2023 23:35:16 +0100 Subject: [PATCH 078/112] lat/lon test still not working properly --- test/test_sediments.jl | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 9a34a1e7d..aeb6fef31 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -96,14 +96,24 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd set_defaults!(biogeochemistry.underlying_biogeochemistry, model) - simulation = Simulation(model, Δt = 50, stop_time = 1day) + N₀ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + + [time_step!(model, 1) for i in 1:3] + + N₁ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + + # conservations + + rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) + + @test isapprox(N₀, N₁; rtol) + + simulation = Simulation(model, Δt = 1, stop_iteration = 1) intercepted_tendencies = Tuple(Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) - N₀ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) - run!(simulation) # the model is changing the tracer tendencies @@ -116,13 +126,6 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd # the sediment values are being integrated initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) @test all([any(Array(interior(field)) .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) - - N₁ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) - - # conservations - rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) - @test isapprox(N₀, N₁; rtol) - return nothing end From 3339990522ad04a89bd6c0e297bd16122d08465e Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 11:23:27 +0100 Subject: [PATCH 079/112] half changing the sediment adapts made it not work for CPU --- src/Boundaries/Sediments/Sediments.jl | 2 +- .../Sediments/instant_remineralization.jl | 6 +++--- src/Boundaries/Sediments/simple_multi_G.jl | 16 ++++++++-------- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 120ab9f3d..43fbd1314 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -37,7 +37,7 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - sediment, bgc.underlying_biogeochemistry, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ) + sediment, bgc.underlying_biogeochemistry, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ, sediment.tendencies.Gⁿ) return nothing end diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index e4c76337b..7270160b5 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -78,7 +78,7 @@ adapt_structure(to, sediment::InstantRemineralisation) = sediment.burial_efficiency_constant2, sediment.burial_efficiency_half_saturaiton, adapt(to, sediment.fields), - adapt(to, sediment.tendencies), + nothing, adapt(to, sediment.bottom_indices)) sediment_tracers(::InstantRemineralisation) = (:N_storage, ) @@ -86,7 +86,7 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st @inline bottom_index_array(sediment::InstantRemineralisation) = sediment.bottom_indices -@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, tendencies) +@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, tendencies, sediment_tendencies) i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) @@ -98,7 +98,7 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st burial_efficiency = sediment.burial_efficiency_constant1 + sediment.burial_efficiency_constant2 * ((flux * 6.56) / (7 + flux * 6.56)) ^ 2 # sediment evolution - @inbounds sediment.tendencies.Gⁿ.N_storage[i, j, 1] = burial_efficiency * flux + @inbounds sediment_tendencies.N_storage[i, j, 1] = burial_efficiency * flux @inbounds tendencies[remineralisation_receiver(bgc)][i, j, k] += flux * (1 - burial_efficiency) / Δz end diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 15c278827..3d9f228f7 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -159,7 +159,7 @@ adapt_structure(to, sediment::SimpleMultiG) = sediment.anoxic_params, sediment.solid_dep_params, adapt(to, sediment.fields), - adapt(to, sediment.tendencies.Gⁿ), + nothing, adapt(to, sediment.bottom_indices)) sediment_tracers(::SimpleMultiG) = (:C_slow, :C_fast, :C_ref, :N_slow, :N_fast, :N_ref) @@ -172,7 +172,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, @inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices -@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, tendencies) +@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, tendencies, sediment_tendencies) i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) @@ -197,13 +197,13 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, reactivity = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1]) # sediment evolution - sediment.tendencies.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow - sediment.tendencies.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast - sediment.tendencies.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition + sediment_tendencies.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow + sediment_tendencies.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast + sediment_tendencies.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition - sediment.tendencies.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow - sediment.tendencies.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast - sediment.tendencies.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition + sediment_tendencies.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow + sediment_tendencies.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast + sediment_tendencies.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition # efflux/influx O₂ = tracers.O₂[i, j, k] From 98e47e6e36eb9a4beae10b273819d79529d821c5 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 11:30:30 +0100 Subject: [PATCH 080/112] remove column diffusion timescale references --- docs/src/model_components/utils.md | 5 ---- examples/data_forced.jl | 1 - src/OceanBioME.jl | 2 +- test/test_sediments.jl | 42 +++++++++++++----------------- 4 files changed, 19 insertions(+), 31 deletions(-) diff --git a/docs/src/model_components/utils.md b/docs/src/model_components/utils.md index 1adfd85dc..70870a951 100644 --- a/docs/src/model_components/utils.md +++ b/docs/src/model_components/utils.md @@ -8,11 +8,6 @@ We have added a few additional utilities which extend the capabilities of Oceana wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_advection_timescale = column_advection_timescale) simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) ``` -Additionally, in a column model you may have a functional definition for the viscosity, so we define an additional diffusion timescale function: -```julia -wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) -simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) -``` Finally, sinking may be more limiting than the normal advective CFL conditions so, we have an additional cell advection timescale defined for 3D models: ```julia wizard = TimeStepWizard(cfl = 0.6, diffusive_cfl = 0.5, max_change = 1.5, min_change = 0., cell_advection_timescale = sinking_advection_timescale) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 506befc75..eb00268a5 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -106,7 +106,6 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 1.5, min_change = 0.75, - cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 3fdf491ad..92ade5aed 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -23,7 +23,7 @@ export TwoBandPhotosyntheticallyActiveRadiation export Boundaries, Sediments, GasExchange, FlatSediment # Utilities -export column_advection_timescale, column_diffusion_timescale, sinking_advection_timescale, Budget +export column_advection_timescale, sinking_advection_timescale, Budget # Positivity preservaiton utilities export ScaleNegativeTracers, ZeroNegativeTracers diff --git a/test/test_sediments.jl b/test/test_sediments.jl index aeb6fef31..aeaf02839 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -34,7 +34,7 @@ set_defaults!(::LOBSTER, model) = sPOM = 0.2299, bPOM = 0.0103) -set_defaults!(::LOBSTER{<:Any, Val{(true, true, true)}}, model) = +set_defaults!(::lobster_variable_redfield, model) = set!(model, P = 0.4686, Z = 0.5363, NO₃ = 2.3103, NH₄ = 0.0010, DIC = 2106.9, Alk = 2408.9, @@ -43,13 +43,6 @@ set_defaults!(::LOBSTER{<:Any, Val{(true, true, true)}}, model) = sPON = 0.2299, sPOC = 1.5080, bPON = 0.0103, bPOC = 0.0781) -set_defaults!(::LOBSTER{<:Any, Val{(false, false, true)}}, model) = - set!(model, P = 0.4686, Z = 0.5363, - NO₃ = 2.3103, NH₄ = 0.0010, - DOC = 5.3390, DON = 0.8115, - sPON = 0.2299, sPOC = 1.5080, - bPON = 0.0103, bPOC = 0.0781) - set_defaults!(::NutrientPhytoplanktonZooplanktonDetritus, model) = set!(model, N = 2.3, P = 0.4, Z = 0.5, D = 0.2) @@ -66,7 +59,6 @@ total_nitrogen(::LOBSTER, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.DOM))) + sum(Array(interior(model.tracers.sPOM))) + sum(Array(interior(model.tracers.bPOM))) - total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.NH₄))) + sum(Array(interior(model.tracers.P))) + @@ -74,7 +66,6 @@ total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tr sum(Array(interior(model.tracers.DON))) + sum(Array(interior(model.tracers.sPON))) + sum(Array(interior(model.tracers.bPON))) - total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(Array(interior(model.tracers.N))) + sum(Array(interior(model.tracers.P))) + sum(Array(interior(model.tracers.Z))) + @@ -96,24 +87,14 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd set_defaults!(biogeochemistry.underlying_biogeochemistry, model) - N₀ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) - - [time_step!(model, 1) for i in 1:3] - - N₁ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) - - # conservations - - rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) - - @test isapprox(N₀, N₁; rtol) - - simulation = Simulation(model, Δt = 1, stop_iteration = 1) + simulation = Simulation(model, Δt = 50, stop_time = 1day) intercepted_tendencies = Tuple(Array(interior(field)) for field in values(TracerFields(keys(model.tracers), grid))) simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) + N₀ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + run!(simulation) # the model is changing the tracer tendencies @@ -126,6 +107,13 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd # the sediment values are being integrated initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) @test all([any(Array(interior(field)) .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) + + N₁ = CUDA.@allowscalar total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) + + # conservations + rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) + @test isapprox(N₀, N₁; rtol) + return nothing end @@ -154,13 +142,19 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill 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), false, true), + variable_redfield = ifelse(isa(sediment_model, SimpleMultiG), true, false), sediment_model)) # get rid of incompatible combinations run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) + if isa(sediment_model, SimpleMultiG) && isa(grid, LatitudeLongitudeGrid) && isa(architecture, CPU) + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin + @test false broken = true + end + end + if run @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))" @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" test_flat_sediment(grid, biogeochemistry, model; timestepper) From 181aa5efe58b380b059b30deaafd8f6089e801e6 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 12:10:48 +0100 Subject: [PATCH 081/112] fixed sediment test and data forced example --- examples/data_forced.jl | 13 +++++++++++-- test/test_sediments.jl | 2 +- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index eb00268a5..a6abb4fd0 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -20,6 +20,11 @@ using OceanBioME using Oceananigans, Random, Printf, NetCDF, Interpolations, DataDeps using Oceananigans.Units +using Oceananigans.Fields: FunctionField + +import Oceananigans.TurbulenceClosures: maximum_numeric_diffusivity + +maximum_numeric_diffusivity(κ::NamedTuple) = maximum(maximum.(values(κ))) const year = years = 365days # just for these idealised cases nothing #hide @@ -73,8 +78,12 @@ biogeochemistry = LOBSTER(; grid, CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = s_function) -model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), +clock = Clock(; time = 0.0) + +κ = FunctionField{Center, Center, Center}(κₜ, grid; clock) + +model = NonhydrostaticModel(; grid, clock, + closure = ScalarDiffusivity(ν = κ, κ = κ), biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),)) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index aeaf02839..86f47f303 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -149,7 +149,7 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) - if isa(sediment_model, SimpleMultiG) && isa(grid, LatitudeLongitudeGrid) && isa(architecture, CPU) + if isa(sediment_model, SimpleMultiG) && isa(grid, LatitudeLongitudeGrid) && isa(architecture, GPU) @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin @test false broken = true end From 2187c9cdc486adca3c797488a27f3cae0c1a7e2f Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 12:14:08 +0100 Subject: [PATCH 082/112] missed one --- examples/data_forced.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index a6abb4fd0..3ed220a1a 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -25,6 +25,7 @@ using Oceananigans.Fields: FunctionField import Oceananigans.TurbulenceClosures: maximum_numeric_diffusivity maximum_numeric_diffusivity(κ::NamedTuple) = maximum(maximum.(values(κ))) +maximum_numeric_diffusivity(κ::FunctionField) = maximum(κ) const year = years = 365days # just for these idealised cases nothing #hide @@ -53,6 +54,11 @@ t_function(x, y, z, t) = temperature_itp(mod(t, 364days)) s_function(x, y, z, t) = salinity_itp(mod(t, 364days)) surface_PAR(x, y, t) = PAR_itp(mod(t, 364days)) κₜ(x, y, z, t) = 2e-2 * max(1 - (z + mld_itp(mod(t, 364days)) / 2)^2 / (mld_itp(mod(t, 364days)) / 2)^2, 0) + 1e-4 + +clock = Clock(; time = 0.0) + +κ = FunctionField{Center, Center, Center}(κₜ, grid; clock) + nothing #hide # ## Grid and PAR field @@ -78,10 +84,6 @@ biogeochemistry = LOBSTER(; grid, CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = s_function) -clock = Clock(; time = 0.0) - -κ = FunctionField{Center, Center, Center}(κₜ, grid; clock) - model = NonhydrostaticModel(; grid, clock, closure = ScalarDiffusivity(ν = κ, κ = κ), biogeochemistry, From 042250808a3e8f2bda38bf104a0f36d4f1d25f2d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 12:49:02 +0100 Subject: [PATCH 083/112] oops --- examples/data_forced.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 3ed220a1a..80f28de37 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -55,13 +55,9 @@ s_function(x, y, z, t) = salinity_itp(mod(t, 364days)) surface_PAR(x, y, t) = PAR_itp(mod(t, 364days)) κₜ(x, y, z, t) = 2e-2 * max(1 - (z + mld_itp(mod(t, 364days)) / 2)^2 / (mld_itp(mod(t, 364days)) / 2)^2, 0) + 1e-4 -clock = Clock(; time = 0.0) - -κ = FunctionField{Center, Center, Center}(κₜ, grid; clock) - nothing #hide -# ## Grid and PAR field +# ## Grid and diffusivity field # Define the grid (in this case a non uniform grid for better resolution near the surface) and an extra Oceananigans field for the PAR to be stored in Nz = 33 Lz = 600meters @@ -74,6 +70,10 @@ z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1) grid = RectilinearGrid(size = (1, 1, Nz), x = (0, 20meters), y = (0, 20meters), z = z_faces) +clock = Clock(; time = 0.0) + +κ = FunctionField{Center, Center, Center}(κₜ, grid; clock) + # ## Biogeochemical and Oceananigans model # Here we instantiate the LOBSTER model with carbonate chemistry and a surface flux of DIC (CO₂) From b133ef0f4741377bcfbf44aa99e3fce2172babcd Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 14:25:03 +0100 Subject: [PATCH 084/112] corrected run criteria on GPU --- test/test_sediments.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 86f47f303..dcb9ab636 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -153,6 +153,8 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin @test false broken = true end + + run = false end if run From 37a21c1e6930c13686e8ce0b8cbde643082f6584 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 14:26:24 +0100 Subject: [PATCH 085/112] corrected run criteria on GPU --- test/test_sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index dcb9ab636..03f011046 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -149,7 +149,7 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) - if isa(sediment_model, SimpleMultiG) && isa(grid, LatitudeLongitudeGrid) && isa(architecture, GPU) + if run == true && isa(sediment_model, SimpleMultiG) && (isa(grid, LatitudeLongitudeGrid) || isa(grid, ImmersedBoundaryGrid)) && isa(architecture, GPU) @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin @test false broken = true end From 16094c9185332aabc65ca9d91e19e7da9edf9eff Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 17:02:38 +0100 Subject: [PATCH 086/112] fix box model --- src/BoxModel/boxmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BoxModel/boxmodel.jl b/src/BoxModel/boxmodel.jl index 935e944db..25dc71a6e 100644 --- a/src/BoxModel/boxmodel.jl +++ b/src/BoxModel/boxmodel.jl @@ -77,7 +77,7 @@ function calculate_tendencies!(model) for tracer in required_biogeochemical_tracers(model.biogeochemistry) if !(tracer == :T) - @inbounds getproperty(Gⁿ, tracer) .= model.biogeochemistry(Val(tracer), 0.0, 0.0, 0.0, model.clock.time, model.values[1]...) + model.forcing[tracer](model.clock.time, model.values[1]...) + @inbounds getproperty(Gⁿ, tracer) .= model.biogeochemistry.underlying_biogeochemistry(Val(tracer), 0.0, 0.0, 0.0, model.clock.time, model.values[1]...) + model.forcing[tracer](model.clock.time, model.values[1]...) end end From f4d17ed998b29a2d254ed3fb55e66fbd55b26975 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 17:04:54 +0100 Subject: [PATCH 087/112] fix box model 2 --- src/BoxModel/boxmodel.jl | 2 +- src/OceanBioME.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/BoxModel/boxmodel.jl b/src/BoxModel/boxmodel.jl index 25dc71a6e..935e944db 100644 --- a/src/BoxModel/boxmodel.jl +++ b/src/BoxModel/boxmodel.jl @@ -77,7 +77,7 @@ function calculate_tendencies!(model) for tracer in required_biogeochemical_tracers(model.biogeochemistry) if !(tracer == :T) - @inbounds getproperty(Gⁿ, tracer) .= model.biogeochemistry.underlying_biogeochemistry(Val(tracer), 0.0, 0.0, 0.0, model.clock.time, model.values[1]...) + model.forcing[tracer](model.clock.time, model.values[1]...) + @inbounds getproperty(Gⁿ, tracer) .= model.biogeochemistry(Val(tracer), 0.0, 0.0, 0.0, model.clock.time, model.values[1]...) + model.forcing[tracer](model.clock.time, model.values[1]...) end end diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 92ade5aed..c6d106e6a 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -120,6 +120,8 @@ update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modi @inline biogeochemical_transition(i, j, k, grid, bgc::Biogeochemistry, val_tracer_name, clock, fields) = biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) +@inoine (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) + function update_biogeochemical_state!(bgc::Biogeochemistry, model) update_biogeochemical_state!(model, bgc.modifiers) update_biogeochemical_state!(model, bgc.light_attenuation) From 0cb434ee21abcb58902fdde4ae6b9b2865043f98 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 17:05:42 +0100 Subject: [PATCH 088/112] typo --- src/OceanBioME.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index c6d106e6a..ac61ed17c 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -120,7 +120,7 @@ update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modi @inline biogeochemical_transition(i, j, k, grid, bgc::Biogeochemistry, val_tracer_name, clock, fields) = biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) -@inoine (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) +@inline (bgc::Biogeochemistry)(args...) = bgc.underlying_biogeochemistry(args...) function update_biogeochemical_state!(bgc::Biogeochemistry, model) update_biogeochemical_state!(model, bgc.modifiers) From ba39a35d1fb83b3d0a6eff57d46c310d4632cf4c Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 17:08:45 +0100 Subject: [PATCH 089/112] hopefully reduced sediment kernel parameter size --- src/Boundaries/Sediments/Sediments.jl | 7 +- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 5 +- src/Models/AdvectedPopulations/NPZD.jl | 3 +- src/TimeSteppers/TimeSteppers.jl | 7 + src/TimeSteppers/positive_conservative_AB2.jl | 162 ++++++++++++++++++ 5 files changed, 181 insertions(+), 3 deletions(-) create mode 100644 src/TimeSteppers/TimeSteppers.jl create mode 100644 src/TimeSteppers/positive_conservative_AB2.jl diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 43fbd1314..733a95aaf 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -37,11 +37,12 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - sediment, bgc.underlying_biogeochemistry, model.grid, model.advection, model.tracers, model.timestepper.Gⁿ, sediment.tendencies.Gⁿ) + sediment, bgc.underlying_biogeochemistry, model.grid, sinking_advection(bgc.underlying_biogeochemistry, model.advection), model.tracers, model.timestepper.Gⁿ, sediment.tendencies.Gⁿ) return nothing end + @kernel function store_flat_tendencies!(G⁻, G⁰) i, j = @index(Global, NTuple) @inbounds G⁻[i, j, 1] = G⁰[i, j, 1] @@ -50,6 +51,10 @@ end @inline nitrogen_flux() = 0 @inline carbon_flux() = 0 @inline remineralisation_receiver() = nothing +@inline sinking_tracers() = nothing + +@inline sinking_advection(bgc, advection) = advection +@inline sinking_advection(bgc, advection::NamedTuple) = advection[sinking_tracers(bgc)] @inline nitrogen_flux(i, j, k, grid, advection, bgc::Biogeochemistry, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc.underlying_biogeochemistry, tracers) @inline carbon_flux(i, j, k, grid, advection, bgc::Biogeochemistry, tracers) = carbon_flux(i, j, k, grid, advection, bgc.underlying_biogeochemistry, tracers) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 31fbee405..905afbaa8 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -65,7 +65,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 +import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry phytoplankton_preference :: FT @@ -489,4 +489,7 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true @inline conserved_tracers(::LOBSTER) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) @inline conserved_tracers(::lobster_variable_redfield) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) + +@inline sinking_tracers(::LOBSTER) = (:sPOM, :bPOM) +@inline sinking_tracers(::lobster_variable_redfield) = (:sPON, :bPON, :sPOC, :bPOC) end # module diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 5d54b43f9..d605c5b13 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -37,7 +37,7 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, import OceanBioME: maximum_sinking_velocity -import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver +import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers import Adapt: adapt_structure, adapt @@ -332,4 +332,5 @@ adapt_structure(to, npzd::NPZD) = @inline remineralisation_receiver(::NPZD) = :N @inline conserved_tracers(::NPZD) = (:N, :P, :Z, :D) +@inline sinking_tracers(::NPZD) = (:P, :D) end # module diff --git a/src/TimeSteppers/TimeSteppers.jl b/src/TimeSteppers/TimeSteppers.jl new file mode 100644 index 000000000..eb7da2f92 --- /dev/null +++ b/src/TimeSteppers/TimeSteppers.jl @@ -0,0 +1,7 @@ +module TimeSteppers + +using Oceananigans.TimeSteppers: AbstractTimeStepper + +include("positive_conservative_AB2.jl") + +end \ No newline at end of file diff --git a/src/TimeSteppers/positive_conservative_AB2.jl b/src/TimeSteppers/positive_conservative_AB2.jl new file mode 100644 index 000000000..b4260b7ee --- /dev/null +++ b/src/TimeSteppers/positive_conservative_AB2.jl @@ -0,0 +1,162 @@ +using Oceananigans.Fields: FunctionField, location +using Oceananigans.TurbulenceClosures: implicit_step! +using Oceananigans.Utils: @apply_regionally, apply_regionally! + +mutable struct PositiveQuasiAdamsBashforth2TimeStepper{FT, GT, IT} <: AbstractTimeStepper + χ :: FT + previous_Δt :: FT + Gⁿ :: GT + G⁻ :: GT + implicit_solver :: IT + intermediate_step :: GT +end + +""" +PositiveQuasiAdamsBashforth2TimeStepper(grid, tracers, + χ = 0.1; + implicit_solver = nothing, + Gⁿ = TendencyFields(grid, tracers), + G⁻ = TendencyFields(grid, tracers)) + +Return a 2nd-order quasi Adams-Bashforth (AB2) time stepper (`QuasiAdamsBashforth2TimeStepper`) +on `grid`, with `tracers`, and AB2 parameter `χ`. The tendency fields `Gⁿ` and `G⁻` can be +specified via optional `kwargs`. + +The 2nd-order quasi Adams-Bashforth timestepper steps forward the state `Uⁿ` by `Δt` via + +```julia +Uⁿ⁺¹ = Uⁿ + Δt * [(3/2 + χ) * Gⁿ - (1/2 + χ) * Gⁿ⁻¹] +``` + +where `Uⁿ` is the state at the ``n``-th timestep, `Gⁿ` is the tendency +at the ``n``-th timestep, and `Gⁿ⁻¹` is the tendency at the previous +timestep (`G⁻`). + +!!! note "First timestep" + For the first timestep, since there are no saved tendencies from the previous timestep, + the `QuasiAdamsBashforth2TimeStepper` performs an Euler timestep: + + ```julia + Uⁿ⁺¹ = Uⁿ + Δt * Gⁿ + ``` +""" +function PositiveQuasiAdamsBashforth2TimeStepper(grid, tracers, + χ = 0.1; + implicit_solver::IT = nothing, + Gⁿ = TendencyFields(grid, tracers), + G⁻ = TendencyFields(grid, tracers), + Cⁿ⁺¹ = TendencyFields(grid, tracers)) where IT + + FT = eltype(grid) + GT = typeof(Gⁿ) + + return PositiveQuasiAdamsBashforth2TimeStepper{FT, GT, IT}(χ, Inf, Gⁿ, G⁻, implicit_solver, Cⁿ⁺¹) +end + +function reset!(timestepper::PositiveQuasiAdamsBashforth2TimeStepper) + timestepper.previous_Δt = Inf + return nothing +end + +##### +##### Time steppping +##### + +""" + time_step!(model::AbstractModel{<:PositiveQuasiAdamsBashforth2TimeStepper}, Δt; euler=false) + +Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and +pressure-correction substep. Setting `euler=true` will take a forward Euler time step. +""" +function time_step!(model::AbstractModel{<:PositiveQuasiAdamsBashforth2TimeStepper}, Δt; callbacks = [], euler=false) + Δt == 0 && @warn "Δt == 0 may cause model blowup!" + + # Shenanigans for properly starting the AB2 loop with an Euler step + euler = euler || (Δt != model.timestepper.previous_Δt) + + χ = ifelse(euler, convert(eltype(model.grid), -0.5), model.timestepper.χ) + + if euler + @debug "Taking a forward Euler step." + # Ensure zeroing out all previous tendency fields to avoid errors in + # case G⁻ includes NaNs. See https://github.com/CliMA/Oceananigans.jl/issues/2259 + for field in model.timestepper.G⁻ + !isnothing(field) && @apply_regionally fill!(field, 0) + end + end + + model.timestepper.previous_Δt = Δt + + # Be paranoid and update state at iteration 0 + model.clock.iteration == 0 && update_state!(model, callbacks) + + @apply_regionally calculate_tendencies!(model, callbacks) + + ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities. + calculate_pressure_correction!(model, Δt) + + @apply_regionally correct_velocities_and_store_tendecies!(model, Δt) + + tick!(model.clock, Δt) + update_state!(model, callbacks) + step_lagrangian_particles!(model, Δt) + + return nothing +end + +function correct_velocities_and_store_tendecies!(model, Δt) + pressure_correct_velocities!(model, Δt) + store_tendencies!(model) + return nothing +end + +##### +##### Time stepping in each step +##### + +""" Generic implementation. """ +function ab2_step!(model, Δt, χ) + + workgroup, worksize = work_layout(model.grid, :xyz) + arch = model.architecture + step_field_kernel! = ab2_step_field!(device(arch), workgroup, worksize) + model_fields = prognostic_fields(model) + + for (i, field) in enumerate(model_fields) + + step_field_kernel!(field, Δt, χ, + model.timestepper.Gⁿ[i], + model.timestepper.G⁻[i]) + + # TODO: function tracer_index(model, field_index) = field_index - 3, etc... + tracer_index = Val(i - 3) # assumption + + implicit_step!(field, + model.timestepper.implicit_solver, + model.closure, + model.diffusivity_fields, + tracer_index, + model.clock, + Δt) + end + + return nothing +end + +""" +Time step velocity fields via the 2nd-order quasi Adams-Bashforth method + + `U^{n+1} = U^n + Δt ((3/2 + χ) * G^{n} - (1/2 + χ) G^{n-1})` + +""" +@kernel function ab2_step_field!(u, Δt, χ, Gⁿ, G⁻) + i, j, k = @index(Global, NTuple) + + T = eltype(u) + one_point_five = convert(T, 1.5) + oh_point_five = convert(T, 0.5) + + @inbounds u[i, j, k] += Δt * ((one_point_five + χ) * Gⁿ[i, j, k] - (oh_point_five + χ) * G⁻[i, j, k]) +end + +@kernel ab2_step_field!(::FunctionField, Δt, χ, Gⁿ, G⁻) = nothing From 5abb4217ab984f2cd32f2bee8742201c68f7e59b Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 19 Sep 2023 17:10:05 +0100 Subject: [PATCH 090/112] renamed lobster variable redfield --- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 22 +++++++++---------- test/test_sediments.jl | 6 ++--- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 905afbaa8..11c04d25b 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -450,7 +450,7 @@ include("carbonate_chemistry.jl") include("oxygen_chemistry.jl") include("variable_redfield.jl") -const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true)}, <:Any}, +const VariableRedfieldLobster = Union{LOBSTER{<:Any, <:Val{(false, false, true)}, <:Any}, LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any}, LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}, LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}} @@ -463,13 +463,13 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true @inline redfield(::Union{Val{:sPOM}, Val{:bPOM}, Val{:DOM}}, bgc::LOBSTER) = bgc.organic_redfield @inline redfield(::Union{Val{:sPOC}, Val{:bPOC}, Val{:DOC}, Val{:DIC}}, bgc::LOBSTER) = 1 -@inline redfield(i, j, k, ::Val{:sPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.sPOC[i, j, k] / tracers.sPON[i, j, k] -@inline redfield(i, j, k, ::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.bPOC[i, j, k] / tracers.bPON[i, j, k] -@inline redfield(i, j, k, ::Val{:DON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.DOC[i, j, k] / tracers.DON[i, j, k] +@inline redfield(i, j, k, ::Val{:sPON}, bgc::VariableRedfieldLobster, tracers) = @inbounds tracers.sPOC[i, j, k] / tracers.sPON[i, j, k] +@inline redfield(i, j, k, ::Val{:bPON}, bgc::VariableRedfieldLobster, tracers) = @inbounds tracers.bPOC[i, j, k] / tracers.bPON[i, j, k] +@inline redfield(i, j, k, ::Val{:DON}, bgc::VariableRedfieldLobster, tracers) = @inbounds tracers.DOC[i, j, k] / tracers.DON[i, j, k] -@inline redfield(::Val{:sPON}, bgc::lobster_variable_redfield, tracers) = tracers.sPOC / tracers.sPON -@inline redfield(::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = tracers.bPOC / tracers.bPON -@inline redfield(::Val{:DON}, bgc::lobster_variable_redfield, tracers) = tracers.DOC / tracers.DON +@inline redfield(::Val{:sPON}, bgc::VariableRedfieldLobster, tracers) = tracers.sPOC / tracers.sPON +@inline redfield(::Val{:bPON}, bgc::VariableRedfieldLobster, tracers) = tracers.bPOC / tracers.bPON +@inline redfield(::Val{:DON}, bgc::VariableRedfieldLobster, tracers) = tracers.DOC / tracers.DON @inline nitrogen_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = sinking_flux(i, j, k, grid, advection, Val(:sPOM), bgc, tracers) + @@ -477,19 +477,19 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true @inline carbon_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * redfield(Val(:sPOM), bgc) -@inline nitrogen_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) = +@inline nitrogen_flux(i, j, k, grid, advection, bgc::VariableRedfieldLobster, tracers) = sinking_flux(i, j, k, grid, advection, Val(:sPON), bgc, tracers) + sinking_flux(i, j, k, grid, advection, Val(:bPON), bgc, tracers) -@inline carbon_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) = +@inline carbon_flux(i, j, k, grid, advection, bgc::VariableRedfieldLobster, tracers) = 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) -@inline conserved_tracers(::lobster_variable_redfield) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) +@inline conserved_tracers(::VariableRedfieldLobster) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) @inline sinking_tracers(::LOBSTER) = (:sPOM, :bPOM) -@inline sinking_tracers(::lobster_variable_redfield) = (:sPON, :bPON, :sPOC, :bPOC) +@inline sinking_tracers(::VariableRedfieldLobster) = (:sPON, :bPON, :sPOC, :bPOC) end # module diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 03f011046..c372b365a 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -9,7 +9,7 @@ using Oceananigans.Fields: TracerFields using Oceananigans.Operators: volume, Azᶠᶜᶜ -using OceanBioME.LOBSTERModel: lobster_variable_redfield +using OceanBioME.LOBSTERModel: VariableRedfieldLobster function intercept_tendencies!(model, intercepted_tendencies) for tracer in keys(model.tracers) @@ -34,7 +34,7 @@ set_defaults!(::LOBSTER, model) = sPOM = 0.2299, bPOM = 0.0103) -set_defaults!(::lobster_variable_redfield, model) = +set_defaults!(::VariableRedfieldLobster, model) = set!(model, P = 0.4686, Z = 0.5363, NO₃ = 2.3103, NH₄ = 0.0010, DIC = 2106.9, Alk = 2408.9, @@ -59,7 +59,7 @@ total_nitrogen(::LOBSTER, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.DOM))) + sum(Array(interior(model.tracers.sPOM))) + sum(Array(interior(model.tracers.bPOM))) -total_nitrogen(::lobster_variable_redfield, model) = sum(Array(interior(model.tracers.NO₃))) + +total_nitrogen(::VariableRedfieldLobster, model) = sum(Array(interior(model.tracers.NO₃))) + sum(Array(interior(model.tracers.NH₄))) + sum(Array(interior(model.tracers.P))) + sum(Array(interior(model.tracers.Z))) + From 474647d74674c7314af2721e5118b9e810c1d9b7 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 13:54:01 +0100 Subject: [PATCH 091/112] error in kelp advection --- src/Models/Individuals/SLatissima.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 839a3f6ce..4aac9bd1e 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -30,6 +30,7 @@ using KernelAbstractions.Extras.LoopInfo: @unroll using Oceananigans.Operators: volume using Oceananigans.Grids: AbstractGrid using Oceananigans.Fields: fractional_indices, _interpolate, datatuple +using Oceananigans.Models: total_velocities import Adapt: adapt_structure, adapt import Oceananigans.Biogeochemistry: update_tendencies! @@ -342,13 +343,13 @@ function update_lagrangian_particle_properties!(particles::SLatissima, model, bg advect_particles_kernel!((x = particles.x, y = particles.y, z = particles.z), 1.0, model.grid, Δt, - datatuple(model.velocities)) + total_velocities(model)) update_particle_properties_kernel! = _update_lagrangian_particle_properties!(device(arch), workgroup, worksize) update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying_biogeochemistry, model.grid, - model.velocities, model.tracers, model.clock, Δt) + total_velocities(model), model.tracers, model.clock, Δt) particles.custom_dynamics(particles, model, bgc, Δt) end From 28559362ac039b6c32ab7863d51abde294a01813 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 16:36:30 +0100 Subject: [PATCH 092/112] changed params from named tuple to tuple --- src/Boundaries/Sediments/simple_multi_G.jl | 94 +++++++++------------- 1 file changed, 40 insertions(+), 54 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 3d9f228f7..a3d7bddde 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -55,11 +55,10 @@ end fast_fraction = 0.74, slow_fraction = 0.26, refactory_fraction = 0.1, - nitrate_oxidation_params = (A = - 1.9785, B = 0.2261, C = -0.0615, D = -0.0289, E = - 0.36109, F = - 0.0232), - denitrification_params = (A = - 3.0790, B = 1.7509, C = 0.0593, D = - 0.1923, E = 0.0604, F = 0.0662), - anoxic_params = (A = - 3.9476, B = 2.6269, C = - 0.2426, D = -1.3349, E = 0.1826, F = - 0.0143), - depth = abs(znodes(grid, Face())[1]), - solid_dep_params = (A = 0.233, B = 0.336, C = 982, D = - 1.548, depth = depth)) + nitrate_oxidation_params = (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232), + denitrification_params = (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662), + anoxic_params = (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143), + solid_dep_params = (0.233, 0.336, 982.0, - 1.548)) Return a single-layer "multi G" sediment model (`SimpleMultiG`) on `grid`, where parameters can be optionally specified. @@ -97,30 +96,10 @@ function SimpleMultiG(; grid, fast_fraction = 0.74, slow_fraction = 0.26, refactory_fraction = 0.1, - nitrate_oxidation_params = (A = - 1.9785, - B = 0.2261, - C = -0.0615, - D = -0.0289, - E = - 0.36109, - F = - 0.0232), - denitrification_params = (A = - 3.0790, - B = 1.7509, - C = 0.0593, - D = - 0.1923, - E = 0.0604, - F = 0.0662), - anoxic_params = (A = - 3.9476, - B = 2.6269, - C = - 0.2426, - D = -1.3349, - E = 0.1826, - F = - 0.0143), - depth = abs(znodes(grid, Face())[1]), - solid_dep_params = (A = 0.233, - B = 0.336, - C = 982, - D = - 1.548, - depth = depth)) + nitrate_oxidation_params = (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232), + denitrification_params = (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662), + anoxic_params = (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143), + solid_dep_params = (0.233, 0.336, 982.0, - 1.548)) @warn "Sediment models are an experimental feature and have not yet been validated." @info "This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC." @@ -147,17 +126,17 @@ function SimpleMultiG(; grid, end adapt_structure(to, sediment::SimpleMultiG) = - SimpleMultiG(sediment.fast_decay_rate, - sediment.slow_decay_rate, - sediment.fast_redfield, - sediment.slow_redfield, - sediment.fast_fraction, - sediment.slow_fraction, - sediment.refactory_fraction, - sediment.nitrate_oxidation_params, - sediment.denitrification_params, - sediment.anoxic_params, - sediment.solid_dep_params, + SimpleMultiG(adapt(to, sediment.fast_decay_rate), + adapt(to, sediment.slow_decay_rate), + adapt(to, sediment.fast_redfield), + adapt(to, sediment.slow_redfield), + adapt(to, sediment.fast_fraction), + adapt(to, sediment.slow_fraction), + adapt(to, sediment.refactory_fraction), + adapt(to, sediment.nitrate_oxidation_params), + adapt(to, sediment.denitrification_params), + adapt(to, sediment.anoxic_params), + adapt(to, sediment.solid_dep_params), adapt(to, sediment.fields), nothing, adapt(to, sediment.bottom_indices)) @@ -176,6 +155,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) + depth = @inbounds znodes(grid, Center(), Center(), Center())[k] Δz = zspacing(i, j, k, grid, Center(), Center(), Center()) @@ -209,13 +189,15 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, O₂ = tracers.O₂[i, j, k] NO₃ = tracers.NO₃[i, j, k] NH₄ = tracers.NH₄[i, j, k] + + A, B, C, D, E, F = sediment.nitrate_oxidation_params - pₙᵢₜ = exp(sediment.nitrate_oxidation_params.A + - sediment.nitrate_oxidation_params.B * log(Cᵐⁱⁿ * day) * log(O₂) + - sediment.nitrate_oxidation_params.C * log(Cᵐⁱⁿ * day) ^ 2 + - sediment.nitrate_oxidation_params.D * log(reactivity) * log(NH₄) + - sediment.nitrate_oxidation_params.E * log(Cᵐⁱⁿ * day) + - sediment.nitrate_oxidation_params.F * log(Cᵐⁱⁿ * day) * log(NH₄)) / (Nᵐⁱⁿ * day) + pₙᵢₜ = exp(A + + B * log(Cᵐⁱⁿ * day) * log(O₂) + + C * log(Cᵐⁱⁿ * day) ^ 2 + + D * log(reactivity) * log(NH₄) + + E * log(Cᵐⁱⁿ * day) + + F * log(Cᵐⁱⁿ * day) * log(NH₄)) / (Nᵐⁱⁿ * day) #= pᵈᵉⁿⁱᵗ = exp(sediment.denitrification_params.A + @@ -225,14 +207,18 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, sediment.denitrification_params.E * log(reactivity) ^ 2 + sediment.denitrification_params.F * log(O₂) * log(reactivity)) / (Cᵐⁱⁿ * day) =# - pₐₙₒₓ = exp(sediment.anoxic_params.A + - sediment.anoxic_params.B * log(Cᵐⁱⁿ * day) + - sediment.anoxic_params.C * log(Cᵐⁱⁿ * day) ^ 2 + - sediment.anoxic_params.D * log(reactivity) + - sediment.anoxic_params.E * log(O₂) * log(reactivity) + - sediment.anoxic_params.F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) - - pₛₒₗᵢ = sediment.solid_dep_params.A * (sediment.solid_dep_params.C * sediment.solid_dep_params.depth ^ sediment.solid_dep_params.D) ^ sediment.solid_dep_params.B + + A, B, C, D, E, F = sediment.anoxic_params + + pₐₙₒₓ = exp(A + + B * log(Cᵐⁱⁿ * day) + + C * log(Cᵐⁱⁿ * day) ^ 2 + + D * log(reactivity) + + E * log(O₂) * log(reactivity) + + F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) + + A, B, C, D = sediment.solid_dep_params + pₛₒₗᵢ = A * (C * depth ^ D) ^ B tendencies.NH₄[i, j, k] += Nᵐⁱⁿ * (1 - pₙᵢₜ) / Δz tendencies.NO₃[i, j, k] += Nᵐⁱⁿ * pₙᵢₜ / Δz From b9b42b3174c139ab1bb48cc2e1d0e567efc4ee00 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 16:36:58 +0100 Subject: [PATCH 093/112] improved generality of `adapt_structure` methods --- .../Sediments/instant_remineralization.jl | 6 +- src/Boundaries/gasexchange.jl | 7 +- src/Light/2band.jl | 16 +-- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 60 +++++------ src/Models/AdvectedPopulations/NPZD.jl | 26 ++--- src/Models/Individuals/SLatissima.jl | 102 +++++++++--------- 6 files changed, 110 insertions(+), 107 deletions(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 7270160b5..8e391f70f 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -74,9 +74,9 @@ function InstantRemineralisation(; grid, end adapt_structure(to, sediment::InstantRemineralisation) = - InstantRemineralisation(sediment.burial_efficiency_constant1, - sediment.burial_efficiency_constant2, - sediment.burial_efficiency_half_saturaiton, + InstantRemineralisation(adapt(to, sediment.burial_efficiency_constant1), + adapt(to, sediment.burial_efficiency_constant2), + adapt(to, sediment.burial_efficiency_half_saturaiton), adapt(to, sediment.fields), nothing, adapt(to, sediment.bottom_indices)) diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index dbb0ad06c..53e13c963 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -43,8 +43,11 @@ adapt_structure(to, pCO₂_model::pCO₂) = pCO₂(adapt(to, pCO₂_model.solubi adapt(to, pCO₂_model.carbonate_dissociation), adapt(to, pCO₂_model.boric_acid_dissociation), adapt(to, pCO₂_model.water_dissociaiton), - pCO₂_model.lower_pH_bound, pCO₂_model.upper_pH_bound, - pCO₂_model.boron_ratio, pCO₂_model.thermal_expansion, pCO₂_model.haline_contraction) + adapt(to, pCO₂_model.lower_pH_bound), + adapt(to, pCO₂_model.upper_pH_bound), + adapt(to, pCO₂_model.boron_ratio), + adapt(to, pCO₂_model.thermal_expansion), + adapt(to, pCO₂_model.haline_contraction)) @inline function titrate_alkalinity(H, p) return p.DIC * (p.k¹ * H + 2 * p.k¹ * p.k²) / (H ^ 2 + p.k¹ * H + p.k¹ * p.k²) - diff --git a/src/Light/2band.jl b/src/Light/2band.jl index 35ade4e06..37562902b 100644 --- a/src/Light/2band.jl +++ b/src/Light/2band.jl @@ -133,13 +133,13 @@ show(io::IO, model::TwoBandPhotosyntheticallyActiveRadiation{FT}) where {FT} = p biogeochemical_auxiliary_fields(par::TwoBandPhotosyntheticallyActiveRadiation) = (PAR = par.field, ) adapt_structure(to, par::TwoBandPhotosyntheticallyActiveRadiation) = - TwoBandPhotosyntheticallyActiveRadiation(par.water_red_attenuation, - par.water_blue_attenuation, - par.chlorophyll_red_attenuation, - par.chlorophyll_blue_attenuation, - par.chlorophyll_red_exponent, - par.chlorophyll_blue_exponent, - par.pigment_ratio, - par.phytoplankton_chlorophyll_ratio, + TwoBandPhotosyntheticallyActiveRadiation(adapt(to, par.water_red_attenuation), + adapt(to, par.water_blue_attenuation), + adapt(to, par.chlorophyll_red_attenuation), + adapt(to, par.chlorophyll_blue_attenuation), + adapt(to, par.chlorophyll_red_exponent), + adapt(to, par.chlorophyll_blue_exponent), + adapt(to, par.pigment_ratio), + adapt(to, par.phytoplankton_chlorophyll_ratio), adapt(to, par.field), adapt(to, par.surface_PAR)) diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 11c04d25b..c6a5f678e 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -400,36 +400,36 @@ function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:LOBSTER}, <: end adapt_structure(to, lobster::LOBSTER) = - LOBSTER(lobster.phytoplankton_preference, - lobster.maximum_grazing_rate, - lobster.grazing_half_saturation, - lobster.light_half_saturation, - lobster.nitrate_ammonia_inhibition, - lobster.nitrate_half_saturation, - lobster.ammonia_half_saturation, - lobster.maximum_phytoplankton_growthrate, - lobster.zooplankton_assimilation_fraction, - lobster.zooplankton_mortality, - lobster.zooplankton_excretion_rate, - lobster.phytoplankton_mortality, - lobster.small_detritus_remineralisation_rate, - lobster.large_detritus_remineralisation_rate, - lobster.phytoplankton_exudation_fraction, - lobster.nitrifcaiton_rate, - lobster.ammonia_fraction_of_exudate, - lobster.ammonia_fraction_of_excriment, - lobster.ammonia_fraction_of_detritus, - lobster.phytoplankton_redfield, - lobster.organic_redfield, - lobster.phytoplankton_chlorophyll_ratio, - lobster.organic_carbon_calcate_ratio, - lobster.respiraiton_oxygen_nitrogen_ratio, - lobster.nitrifcation_oxygen_nitrogen_ratio, - lobster.slow_sinking_mortality_fraction, - lobster.fast_sinking_mortality_fraction, - lobster.disolved_organic_breakdown_rate, - lobster.zooplankton_calcite_dissolution, - lobster.optionals, + LOBSTER(adapt(to, lobster.phytoplankton_preference), + adapt(to, lobster.maximum_grazing_rate), + adapt(to, lobster.grazing_half_saturation), + adapt(to, lobster.light_half_saturation), + adapt(to, lobster.nitrate_ammonia_inhibition), + adapt(to, lobster.nitrate_half_saturation), + adapt(to, lobster.ammonia_half_saturation), + adapt(to, lobster.maximum_phytoplankton_growthrate), + adapt(to, lobster.zooplankton_assimilation_fraction), + adapt(to, lobster.zooplankton_mortality), + adapt(to, lobster.zooplankton_excretion_rate), + adapt(to, lobster.phytoplankton_mortality), + adapt(to, lobster.small_detritus_remineralisation_rate), + adapt(to, lobster.large_detritus_remineralisation_rate), + adapt(to, lobster.phytoplankton_exudation_fraction), + adapt(to, lobster.nitrifcaiton_rate), + adapt(to, lobster.ammonia_fraction_of_exudate), + adapt(to, lobster.ammonia_fraction_of_excriment), + adapt(to, lobster.ammonia_fraction_of_detritus), + adapt(to, lobster.phytoplankton_redfield), + adapt(to, lobster.organic_redfield), + adapt(to, lobster.phytoplankton_chlorophyll_ratio), + adapt(to, lobster.organic_carbon_calcate_ratio), + adapt(to, lobster.respiraiton_oxygen_nitrogen_ratio), + adapt(to, lobster.nitrifcation_oxygen_nitrogen_ratio), + adapt(to, lobster.slow_sinking_mortality_fraction), + adapt(to, lobster.fast_sinking_mortality_fraction), + adapt(to, lobster.disolved_organic_breakdown_rate), + adapt(to, lobster.zooplankton_calcite_dissolution), + adapt(to, lobster.optionals), adapt(to, lobster.sinking_velocities)) summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index d605c5b13..31b16f206 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -304,19 +304,19 @@ show(io::IO, model::NPZD) = string(summary(model), " \n", @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w) adapt_structure(to, npzd::NPZD) = - NutrientPhytoplanktonZooplanktonDetritus(npzd.initial_photosynthetic_slope, - npzd.base_maximum_growth, - npzd.nutrient_half_saturation, - npzd.base_respiration_rate, - npzd.phyto_base_mortality_rate, - - npzd.maximum_grazing_rate, - npzd.grazing_half_saturation, - npzd.assimulation_efficiency, - npzd.base_excretion_rate, - npzd.zoo_base_mortality_rate, - - npzd.remineralization_rate, + NutrientPhytoplanktonZooplanktonDetritus(adapt(to, npzd.initial_photosynthetic_slope), + adapt(to, npzd.base_maximum_growth), + adapt(to, npzd.nutrient_half_saturation), + adapt(to, npzd.base_respiration_rate), + adapt(to, npzd.phyto_base_mortality_rate, + + adapt(to, npzd.maximum_grazing_rate), + adapt(to, npzd.grazing_half_saturation), + adapt(to, npzd.assimulation_efficiency), + adapt(to, npzd.base_excretion_rate), + adapt(to, npzd.zoo_base_mortality_rate), + + adapt(to, npzd.remineralization_rate), adapt(to, npzd.sinking_velocities)) diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 4aac9bd1e..122337337 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -201,54 +201,54 @@ Base.@kwdef struct SLatissima{AR, FT, U, T, S, P, F} <: BiogeochemicalParticles latitude :: FT = 57.5 end -adapt_structure(to, kelp::SLatissima) = SLatissima(kelp.architecture, - kelp.growth_rate_adjustement, - kelp.photosynthetic_efficiency, - kelp.minimum_carbon_reserve, - kelp.structural_carbon, - kelp.exudation, - kelp.erosion, - kelp.saturation_irradiance, - kelp.structural_dry_weight_per_area, - kelp.structural_dry_to_wet_weight, - kelp.carbon_reserve_per_carbon, - kelp.nitrogen_reserve_per_nitrogen, - kelp.minimum_nitrogen_reserve, - kelp.maximum_nitrogen_reserve, - kelp.growth_adjustement_2, - kelp.growth_adjustement_1, - kelp.maximum_specific_growth_rate, - kelp.structural_nitrogen, - kelp.photosynthesis_at_ref_temp_1, - kelp.photosynthesis_at_ref_temp_2, - kelp.photosynthesis_ref_temp_1, - kelp.photosynthesis_ref_temp_2, - kelp.photoperiod_1, - kelp.photoperiod_2, - kelp.respiration_at_ref_temp_1, - kelp.respiration_at_ref_temp_2, - kelp.respiration_ref_temp_1, - kelp.respiration_ref_temp_2, - kelp.photosynthesis_arrhenius_temp, - kelp.photosynthesis_low_temp, - kelp.photosynthesis_high_temp, - kelp.photosynthesis_high_arrhenius_temp, - kelp.photosynthesis_low_arrhenius_temp, - kelp.respiration_arrhenius_temp, - kelp.current_speed_for_0p65_uptake, - kelp.nitrate_half_saturation, - kelp.ammonia_half_saturation, - kelp.maximum_nitrate_uptake, - kelp.maximum_ammonia_uptake, - kelp.current_1, - kelp.current_2, - kelp.current_3, - kelp.respiration_reference_A, - kelp.respiration_reference_B, - kelp.exudation_redfield_ratio, - kelp.pescribed_velocity, - kelp.pescribed_temperature, - kelp.pescribed_salinity, +adapt_structure(to, kelp::SLatissima) = SLatissima(adapt(to, kelp.architecture), + adapt(to, kelp.growth_rate_adjustement), + adapt(to, kelp.photosynthetic_efficiency), + adapt(to, kelp.minimum_carbon_reserve), + adapt(to, kelp.structural_carbon), + adapt(to, kelp.exudation), + adapt(to, kelp.erosion), + adapt(to, kelp.saturation_irradiance), + adapt(to, kelp.structural_dry_weight_per_area), + adapt(to, kelp.structural_dry_to_wet_weight), + adapt(to, kelp.carbon_reserve_per_carbon), + adapt(to, kelp.nitrogen_reserve_per_nitrogen), + adapt(to, kelp.minimum_nitrogen_reserve), + adapt(to, kelp.maximum_nitrogen_reserve), + adapt(to, kelp.growth_adjustement_2), + adapt(to, kelp.growth_adjustement_1), + adapt(to, kelp.maximum_specific_growth_rate), + adapt(to, kelp.structural_nitrogen), + adapt(to, kelp.photosynthesis_at_ref_temp_1), + adapt(to, kelp.photosynthesis_at_ref_temp_2), + adapt(to, kelp.photosynthesis_ref_temp_1), + adapt(to, kelp.photosynthesis_ref_temp_2), + adapt(to, kelp.photoperiod_1), + adapt(to, kelp.photoperiod_2), + adapt(to, kelp.respiration_at_ref_temp_1), + adapt(to, kelp.respiration_at_ref_temp_2), + adapt(to, kelp.respiration_ref_temp_1), + adapt(to, kelp.respiration_ref_temp_2), + adapt(to, kelp.photosynthesis_arrhenius_temp), + adapt(to, kelp.photosynthesis_low_temp), + adapt(to, kelp.photosynthesis_high_temp), + adapt(to, kelp.photosynthesis_high_arrhenius_temp), + adapt(to, kelp.photosynthesis_low_arrhenius_temp), + adapt(to, kelp.respiration_arrhenius_temp), + adapt(to, kelp.current_speed_for_0p65_uptake), + adapt(to, kelp.nitrate_half_saturation), + adapt(to, kelp.ammonia_half_saturation), + adapt(to, kelp.maximum_nitrate_uptake), + adapt(to, kelp.maximum_ammonia_uptake), + adapt(to, kelp.current_1), + adapt(to, kelp.current_2), + adapt(to, kelp.current_3), + adapt(to, kelp.respiration_reference_A), + adapt(to, kelp.respiration_reference_B), + adapt(to, kelp.exudation_redfield_ratio), + adapt(to, kelp.pescribed_velocity), + adapt(to, kelp.pescribed_temperature), + adapt(to, kelp.pescribed_salinity), adapt(to, kelp.x), adapt(to, kelp.y), adapt(to, kelp.z), @@ -261,9 +261,9 @@ adapt_structure(to, kelp::SLatissima) = SLatissima(kelp.architecture, adapt(to, kelp.frond_exudation), adapt(to, kelp.nitrogen_erosion), adapt(to, kelp.carbon_erosion), - kelp.custom_dynamics, - kelp.scalefactor, - kelp.latitude) + adapt(to, kelp.custom_dynamics), + adapt(to, kelp.scalefactor), + adapt(to, kelp.latitude)) function update_tendencies!(bgc, particles::SLatissima, model) num_particles = length(particles) From 21622ff4f232f89721bdecbbcbf3fc0007a78eb5 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 17:28:35 +0100 Subject: [PATCH 094/112] updated broken tests --- test/test_sediments.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index c372b365a..4f43428d9 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -149,14 +149,6 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) - if run == true && isa(sediment_model, SimpleMultiG) && (isa(grid, LatitudeLongitudeGrid) || isa(grid, ImmersedBoundaryGrid)) && isa(architecture, GPU) - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" begin - @test false broken = true - end - - run = false - end - if run @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry)) with $(display_name(grid))" @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry)), $(display_name(grid))" test_flat_sediment(grid, biogeochemistry, model; timestepper) From d3f08828adefb5ad95e1b29b55b894c3a1c5a801 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 17:38:20 +0100 Subject: [PATCH 095/112] typo --- 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 31b16f206..95071cf80 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -308,7 +308,7 @@ adapt_structure(to, npzd::NPZD) = adapt(to, npzd.base_maximum_growth), adapt(to, npzd.nutrient_half_saturation), adapt(to, npzd.base_respiration_rate), - adapt(to, npzd.phyto_base_mortality_rate, + adapt(to, npzd.phyto_base_mortality_rate), adapt(to, npzd.maximum_grazing_rate), adapt(to, npzd.grazing_half_saturation), From 2351d64d4661814b4ad57156c8f40ee6ff894af8 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 18:59:19 +0100 Subject: [PATCH 096/112] sign error --- src/Boundaries/Sediments/simple_multi_G.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index a3d7bddde..9f849787b 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -155,7 +155,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) - depth = @inbounds znodes(grid, Center(), Center(), Center())[k] + depth = @inbounds -znodes(grid, Center(), Center(), Center())[k] Δz = zspacing(i, j, k, grid, Center(), Center(), Center()) From b8e1af38c0e3205639c4f21356dc28534d644f6d Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 19:06:09 +0100 Subject: [PATCH 097/112] fixed sediment docs --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index bccc8f06f..0e397e3e5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -53,7 +53,7 @@ model_parameters = (LOBSTER(; grid = BoxModelGrid()), NutrientPhytoplanktonZooplanktonDetritus(; grid = BoxModelGrid()), SLatissima(), TwoBandPhotosyntheticallyActiveRadiation(; grid = BoxModelGrid()), - SimpleMultiG(; grid = BoxModelGrid(), depth = 1000), + SimpleMultiG(; grid = BoxModelGrid()), InstantRemineralisation(; grid = BoxModelGrid()), OCMIP_default, GasExchange(; gas = :CO₂).condition.parameters, From 2cdee7bd4465ddff66bf5c76d9b664fbdb35b4ed Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 21:35:53 +0100 Subject: [PATCH 098/112] removed unneccesary `Array(interior(...))` --- test/test_sediments.jl | 50 ++++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 4f43428d9..1a2214bf4 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -46,30 +46,32 @@ set_defaults!(::VariableRedfieldLobster, model) = set_defaults!(::NutrientPhytoplanktonZooplanktonDetritus, model) = set!(model, N = 2.3, P = 0.4, Z = 0.5, D = 0.2) -total_nitrogen(sed::SimpleMultiG) = sum(Array(interior(sed.fields.N_fast))) + - sum(Array(interior(sed.fields.N_slow))) + - sum(Array(interior(sed.fields.N_ref))) - -total_nitrogen(sed::InstantRemineralisation) = sum(Array(interior(sed.fields.N_storage))) - -total_nitrogen(::LOBSTER, model) = sum(Array(interior(model.tracers.NO₃))) + - sum(Array(interior(model.tracers.NH₄))) + - sum(Array(interior(model.tracers.P))) + - sum(Array(interior(model.tracers.Z))) + - sum(Array(interior(model.tracers.DOM))) + - sum(Array(interior(model.tracers.sPOM))) + - sum(Array(interior(model.tracers.bPOM))) -total_nitrogen(::VariableRedfieldLobster, model) = sum(Array(interior(model.tracers.NO₃))) + - sum(Array(interior(model.tracers.NH₄))) + - sum(Array(interior(model.tracers.P))) + - sum(Array(interior(model.tracers.Z))) + - sum(Array(interior(model.tracers.DON))) + - sum(Array(interior(model.tracers.sPON))) + - sum(Array(interior(model.tracers.bPON))) -total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(Array(interior(model.tracers.N))) + - sum(Array(interior(model.tracers.P))) + - sum(Array(interior(model.tracers.Z))) + - sum(Array(interior(model.tracers.D))) +total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) + + sum(sed.fields.N_slow) + + sum(sed.fields.N_ref) + +total_nitrogen(sed::InstantRemineralisation) = sum(sed.fields.N_storage) + +total_nitrogen(::LOBSTER, model) = sum(model.tracers.NO₃) + + sum(model.tracers.NH₄) + + sum(model.tracers.P) + + sum(model.tracers.Z) + + sum(model.tracers.DOM) + + sum(model.tracers.sPOM) + + sum(model.tracers.bPOM) + +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) + +total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(model.tracers.N) + + sum(model.tracers.P) + + sum(model.tracers.Z) + + sum(model.tracers.D) function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAdamsBashforth2) model = isa(model, NonhydrostaticModel) ? model(; grid, From bc6bb3ea38c8e02eb4c7286e43f3b5f169f413d6 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:02:58 +0100 Subject: [PATCH 099/112] flattened sinking velocities --- src/Boundaries/Sediments/simple_multi_G.jl | 30 +++++++++---------- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 2 +- src/Models/AdvectedPopulations/NPZD.jl | 8 ----- src/Utils/sinking_velocity_fields.jl | 3 +- 4 files changed, 17 insertions(+), 26 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 9f849787b..4ad9be129 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -55,10 +55,10 @@ end fast_fraction = 0.74, slow_fraction = 0.26, refactory_fraction = 0.1, - nitrate_oxidation_params = (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232), - denitrification_params = (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662), - anoxic_params = (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143), - solid_dep_params = (0.233, 0.336, 982.0, - 1.548)) + nitrate_oxidation_params = arch_array(architecture(grid), (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232)), + denitrification_params = arch_array(architecture(grid), (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662)), + anoxic_params = arch_array(architecture(grid), (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143)), + solid_dep_params = arch_array(architecture(grid), (0.233, 0.336, 982.0, - 1.548))) Return a single-layer "multi G" sediment model (`SimpleMultiG`) on `grid`, where parameters can be optionally specified. @@ -89,17 +89,17 @@ Single-layer multi-G sediment model (Float64) ``` """ function SimpleMultiG(; grid, - fast_decay_rate = 2/day, - slow_decay_rate = 0.2/day, - fast_redfield = 0.1509, - slow_redfield = 0.13, - fast_fraction = 0.74, - slow_fraction = 0.26, - refactory_fraction = 0.1, - nitrate_oxidation_params = (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232), - denitrification_params = (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662), - anoxic_params = (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143), - solid_dep_params = (0.233, 0.336, 982.0, - 1.548)) + fast_decay_rate = 2/day, + slow_decay_rate = 0.2/day, + fast_redfield = 0.1509, + slow_redfield = 0.13, + fast_fraction = 0.74, + slow_fraction = 0.26, + refactory_fraction = 0.1, + nitrate_oxidation_params = arch_array(architecture(grid), (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232)), + denitrification_params = arch_array(architecture(grid), (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662)), + anoxic_params = arch_array(architecture(grid), (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143)), + solid_dep_params = arch_array(architecture(grid), (0.233, 0.336, 982.0, - 1.548))) @warn "Sediment models are an experimental feature and have not yet been validated." @info "This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC." diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index c6a5f678e..422bc2b7e 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -389,7 +389,7 @@ const DOM = Union{Val{:DOM}, Val{:DON}} @inline function biogeochemical_drift_velocity(bgc::LOBSTER, ::Val{tracer_name}) where tracer_name if tracer_name in keys(bgc.sinking_velocities) - return bgc.sinking_velocities[tracer_name] + return (u = ZeroField(), v = ZeroField(), w = bgc.sinking_velocities[tracer_name]) else return (u = ZeroField(), v = ZeroField(), w = ZeroField()) end diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 95071cf80..74febd0cc 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -284,14 +284,6 @@ end return phytoplankton_mortality_loss + zooplankton_assimilation_loss + zooplankton_mortality_loss - remineralization end -@inline function biogeochemical_drift_velocity(bgc::NPZD, ::Val{tracer_name}) where tracer_name - if tracer_name in keys(bgc.sinking_velocities) - return bgc.sinking_velocities[tracer_name] - else - return (u = ZeroField(), v = ZeroField(), w = ZeroField()) - end -end - function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:NPZD}, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) diff --git a/src/Utils/sinking_velocity_fields.jl b/src/Utils/sinking_velocity_fields.jl index 2a2805814..90208a206 100644 --- a/src/Utils/sinking_velocity_fields.jl +++ b/src/Utils/sinking_velocity_fields.jl @@ -7,7 +7,6 @@ import Adapt: adapt_structure, adapt function setup_velocity_fields(drift_speeds, grid::AbstractGrid, open_bottom; smoothing_distance = 2) drift_velocities = [] for w in values(drift_speeds) - u, v = maybe_constant_field.((0, 0)) if isa(values(w), Number) w_field = ZFaceField(grid) for k=1:grid.Nz @@ -18,7 +17,7 @@ function setup_velocity_fields(drift_speeds, grid::AbstractGrid, open_bottom; sm error("Provided sinking speeds are an unsuitable value, must be a `NamedTuple` of scalar values (positive = sinking) with keys of tracer names, or `NamedTuple` of `VelocityFields`") end - push!(drift_velocities, (; u, v, w)) + push!(drift_velocities, w) end return NamedTuple{keys(drift_speeds)}(drift_velocities) From 9ccda063622c87a6b4a91e73b7bdcc79ad37d093 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:22:42 +0100 Subject: [PATCH 100/112] reduced complexity of sediment models further --- src/Boundaries/Sediments/Sediments.jl | 9 ++++++++- src/Boundaries/Sediments/instant_remineralization.jl | 3 +++ src/Boundaries/Sediments/simple_multi_G.jl | 3 +++ src/Models/AdvectedPopulations/NPZD.jl | 8 ++++++++ 4 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 733a95aaf..2cba78810 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -35,9 +35,14 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, store_flat_tendencies!, sediment.tendencies.G⁻[i], sediment.tendencies.Gⁿ[i]) end + Biogeochemistry = bgc.underlying_biogeochemistry + launch!(arch, model.grid, :xy, _calculate_tendencies!, - sediment, bgc.underlying_biogeochemistry, model.grid, sinking_advection(bgc.underlying_biogeochemistry, model.advection), model.tracers, model.timestepper.Gⁿ, sediment.tendencies.Gⁿ) + sediment, biogeochemistry, model.grid, + sinking_advection(biogeochemistry, model.advection), + required_tracers(sediment, biogeochemistry, model.tracers), + required_tendencies(sediment, biogeochemistry, model.timestepper.Gⁿ), sediment.tendencies.Gⁿ) return nothing end @@ -52,6 +57,8 @@ end @inline carbon_flux() = 0 @inline remineralisation_receiver() = nothing @inline sinking_tracers() = nothing +@inline required_tracers() = nothing +@inline required_tendencies() = nothing @inline sinking_advection(bgc, advection) = advection @inline sinking_advection(bgc, advection::NamedTuple) = advection[sinking_tracers(bgc)] diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 8e391f70f..774dbfb1f 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -84,6 +84,9 @@ adapt_structure(to, sediment::InstantRemineralisation) = sediment_tracers(::InstantRemineralisation) = (:N_storage, ) sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, ) +@inline required_tracers(::InstantRemineralisation, bgc, tracers) = tracers[(sinking_tracers(bgc)..., remineralisation_receiver(bgc)...)] +@inline required_tendencies(::InstantRemineralisation, bgc, tracers) = tracers[remineralisation_receiver(bgc)] + @inline bottom_index_array(sediment::InstantRemineralisation) = sediment.bottom_indices @kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, tendencies, sediment_tendencies) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 4ad9be129..ee531cbcc 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -149,6 +149,9 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, C_ref = model.fields.C_ref, N_ref = model.fields.N_ref) +@inline required_tracers(::SimpleMultiG, bgc, tracers) = tracers[(:NO₃, :NH₄, :O₂, sinking_tracers(bgc)...)] +@inline required_tendencies(::SimpleMultiG, bgc, tracers) = tracers[(:NO₃, :NH₄, :O₂,, :DIC)] + @inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices @kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, tendencies, sediment_tendencies) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 74febd0cc..35f01465c 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -284,6 +284,14 @@ end return phytoplankton_mortality_loss + zooplankton_assimilation_loss + zooplankton_mortality_loss - remineralization end +@inline function biogeochemical_drift_velocity(bgc::NPZD, ::Val{tracer_name}) where tracer_name + if tracer_name in keys(bgc.sinking_velocities) + return (u = ZeroField(), v = ZeroField(), w = bgc.sinking_velocities[tracer_name]) + else + return (u = ZeroField(), v = ZeroField(), w = ZeroField()) + end +end + function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:NPZD}, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) From e9f9206afb1fe87187d8d00d663a6c24e59fb753 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:24:47 +0100 Subject: [PATCH 101/112] typo --- src/Boundaries/Sediments/simple_multi_G.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index ee531cbcc..d7e3b9697 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -150,7 +150,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, N_ref = model.fields.N_ref) @inline required_tracers(::SimpleMultiG, bgc, tracers) = tracers[(:NO₃, :NH₄, :O₂, sinking_tracers(bgc)...)] -@inline required_tendencies(::SimpleMultiG, bgc, tracers) = tracers[(:NO₃, :NH₄, :O₂,, :DIC)] +@inline required_tendencies(::SimpleMultiG, bgc, tracers) = tracers[(:NO₃, :NH₄, :O₂, :DIC)] @inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices From 36d274163b44eaa0397325ffac8b604375eaf068 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:30:35 +0100 Subject: [PATCH 102/112] another mistake [skip ci] --- src/Boundaries/Sediments/simple_multi_G.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index d7e3b9697..c30192a9b 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -55,10 +55,10 @@ end fast_fraction = 0.74, slow_fraction = 0.26, refactory_fraction = 0.1, - nitrate_oxidation_params = arch_array(architecture(grid), (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232)), - denitrification_params = arch_array(architecture(grid), (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662)), - anoxic_params = arch_array(architecture(grid), (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143)), - solid_dep_params = arch_array(architecture(grid), (0.233, 0.336, 982.0, - 1.548))) + nitrate_oxidation_params = arch_array(architecture(grid), [- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232]), + denitrification_params = arch_array(architecture(grid), [- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662]), + anoxic_params = arch_array(architecture(grid), [- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143]), + solid_dep_params = arch_array(architecture(grid), [0.233, 0.336, 982.0, - 1.548])) Return a single-layer "multi G" sediment model (`SimpleMultiG`) on `grid`, where parameters can be optionally specified. @@ -96,10 +96,10 @@ function SimpleMultiG(; grid, fast_fraction = 0.74, slow_fraction = 0.26, refactory_fraction = 0.1, - nitrate_oxidation_params = arch_array(architecture(grid), (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232)), - denitrification_params = arch_array(architecture(grid), (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662)), - anoxic_params = arch_array(architecture(grid), (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143)), - solid_dep_params = arch_array(architecture(grid), (0.233, 0.336, 982.0, - 1.548))) + nitrate_oxidation_params = arch_array(architecture(grid), [- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232]), + denitrification_params = arch_array(architecture(grid), [- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662]), + anoxic_params = arch_array(architecture(grid), [- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143]), + solid_dep_params = arch_array(architecture(grid), [0.233, 0.336, 982.0, - 1.548])) @warn "Sediment models are an experimental feature and have not yet been validated." @info "This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC." From 70cb036602e9927f02258512ac75a79ba9a75743 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:36:30 +0100 Subject: [PATCH 103/112] typo [skip ci] --- src/Boundaries/Sediments/Sediments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 2cba78810..0351fce90 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -35,7 +35,7 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, store_flat_tendencies!, sediment.tendencies.G⁻[i], sediment.tendencies.Gⁿ[i]) end - Biogeochemistry = bgc.underlying_biogeochemistry + biogeochemistry = bgc.underlying_biogeochemistry launch!(arch, model.grid, :xy, _calculate_tendencies!, From 13314f060b9ce011e4667d6851aac76c6d29a3a7 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:44:01 +0100 Subject: [PATCH 104/112] tidy up --- src/Boundaries/Sediments/Sediments.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 0351fce90..51d8be670 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -42,7 +42,8 @@ function update_tendencies!(bgc, sediment::FlatSediment, model) sediment, biogeochemistry, model.grid, sinking_advection(biogeochemistry, model.advection), required_tracers(sediment, biogeochemistry, model.tracers), - required_tendencies(sediment, biogeochemistry, model.timestepper.Gⁿ), sediment.tendencies.Gⁿ) + required_tendencies(sediment, biogeochemistry, model.timestepper.Gⁿ), + sediment.tendencies.Gⁿ) return nothing end From aecd5b716afc8b4296475442cc0a05b2da9cfec2 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:46:55 +0100 Subject: [PATCH 105/112] minor issue with instant remin --- src/Boundaries/Sediments/instant_remineralization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 774dbfb1f..c41ddcb71 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -84,7 +84,7 @@ adapt_structure(to, sediment::InstantRemineralisation) = sediment_tracers(::InstantRemineralisation) = (:N_storage, ) sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, ) -@inline required_tracers(::InstantRemineralisation, bgc, tracers) = tracers[(sinking_tracers(bgc)..., remineralisation_receiver(bgc)...)] +@inline required_tracers(::InstantRemineralisation, bgc, tracers) = tracers[(sinking_tracers(bgc)..., remineralisation_receiver(bgc))] @inline required_tendencies(::InstantRemineralisation, bgc, tracers) = tracers[remineralisation_receiver(bgc)] @inline bottom_index_array(sediment::InstantRemineralisation) = sediment.bottom_indices From 3b953dd63ed5ed2cdfa149cf0131c68e68aab1f6 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Wed, 20 Sep 2023 22:56:52 +0100 Subject: [PATCH 106/112] minor issue with instant remin --- src/Boundaries/Sediments/instant_remineralization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index c41ddcb71..3828b1629 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -103,7 +103,7 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st # sediment evolution @inbounds sediment_tendencies.N_storage[i, j, 1] = burial_efficiency * flux - @inbounds tendencies[remineralisation_receiver(bgc)][i, j, k] += flux * (1 - burial_efficiency) / Δz + @inbounds tendencies[i, j, k] += flux * (1 - burial_efficiency) / Δz end summary(::InstantRemineralisation{FT}) where {FT} = string("Single-layer instant remineralisaiton ($FT)") From 2dfb211aa9b6287df50f33ca6c0c37305e0b4775 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 21 Sep 2023 11:46:48 +0100 Subject: [PATCH 107/112] updated model implementation page --- docs/src/model_implementation.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 59fcd5185..32933dfc9 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -213,7 +213,7 @@ Another aspect that OceanBioME includes is sediment models. Doing this varies be ```@example implementing using OceanBioME.Boundaries.Sediments: sinking_flux -import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver +import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers @inline nitrogen_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) @@ -221,6 +221,8 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa @inline carbon_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 @inline remineralisation_receiver(::NutrientPhytoplankton) = :N + +@inline siniking_tracers(::NutrientPhytoplankton) = (:P, ) ``` ### Putting it together From 3ddd9470159d197b7a4b81dd4737fe9e5c0fe9c3 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Thu, 21 Sep 2023 14:21:13 +0100 Subject: [PATCH 108/112] typo --- 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 32933dfc9..c8e8d1c51 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -222,7 +222,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa @inline remineralisation_receiver(::NutrientPhytoplankton) = :N -@inline siniking_tracers(::NutrientPhytoplankton) = (:P, ) +@inline sinking_tracers(::NutrientPhytoplankton) = (:P, ) ``` ### Putting it together From 6b7b3fb50c73e63f020fcef5cf48b1781bb578bf Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 23 Sep 2023 17:16:53 +0100 Subject: [PATCH 109/112] Remove redundant fallbacks --- src/Boundaries/Sediments/Sediments.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 51d8be670..77d898ade 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -64,10 +64,6 @@ end @inline sinking_advection(bgc, advection) = advection @inline sinking_advection(bgc, advection::NamedTuple) = advection[sinking_tracers(bgc)] -@inline nitrogen_flux(i, j, k, grid, advection, bgc::Biogeochemistry, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc.underlying_biogeochemistry, tracers) -@inline carbon_flux(i, j, k, grid, advection, bgc::Biogeochemistry, tracers) = carbon_flux(i, j, k, grid, advection, bgc.underlying_biogeochemistry, tracers) -@inline remineralisation_receiver(bgc::Biogeochemistry) = remineralisation_receiver(bgc.underlying_biogeochemistry) - @inline advection_scheme(advection, val_tracer) = advection @inline advection_scheme(advection::NamedTuple, val_tracer::Val{T}) where T = advection[T] From cd992995b7d585d034a4d54cc925cff2eaa13e4a Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Sat, 23 Sep 2023 17:17:21 +0100 Subject: [PATCH 110/112] generalised sinking tracers for NPZD --- 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 35f01465c..1134967c4 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -332,5 +332,5 @@ adapt_structure(to, npzd::NPZD) = @inline remineralisation_receiver(::NPZD) = :N @inline conserved_tracers(::NPZD) = (:N, :P, :Z, :D) -@inline sinking_tracers(::NPZD) = (:P, :D) +@inline sinking_tracers(bgc::NPZD) = keys(bgc.sinking_velocities) end # module From 745501e43beb69bd347ef3f9229465c272a5f81a Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Mon, 25 Sep 2023 12:11:20 +0100 Subject: [PATCH 111/112] bump oceananigans compat entry --- Manifest.toml | 6 +++--- Project.toml | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 3a78e27c1..b4638efb1 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.9.3" manifest_format = "2.0" -project_hash = "b8336a615981d764881d88f0d8193dd7f4667beb" +project_hash = "d35eb03107a4f7f3d176d38b0533075a0ad26402" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -500,9 +500,9 @@ 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 = "8746c619e15950de59160b2d99e61961cfa57fcc" +git-tree-sha1 = "1821dce5269caa6f9dc3df3e1b97928b39f233b5" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.87.2" +version = "0.88.0" [[deps.OffsetArrays]] deps = ["Adapt"] diff --git a/Project.toml b/Project.toml index cbb241cdf..0f4b34186 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,7 @@ StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" Adapt = "3" JLD2 = "0.4" KernelAbstractions = "0.9" -Oceananigans = "0.84.1, 0.85, 0.86, 0.87" +Oceananigans = "0.84.1, 0.85, 0.86, 0.87, 0.88" Roots = "2" StructArrays = "0.4, 0.5, 0.6" julia = "1.9" From 103c9561b785d02356f0441e5393b8e618a025e4 Mon Sep 17 00:00:00 2001 From: Jago Stong-Wright Date: Tue, 3 Oct 2023 13:46:30 +0100 Subject: [PATCH 112/112] removed accidentally committed new files --- src/TimeSteppers/TimeSteppers.jl | 7 - src/TimeSteppers/positive_conservative_AB2.jl | 162 ------------------ 2 files changed, 169 deletions(-) delete mode 100644 src/TimeSteppers/TimeSteppers.jl delete mode 100644 src/TimeSteppers/positive_conservative_AB2.jl diff --git a/src/TimeSteppers/TimeSteppers.jl b/src/TimeSteppers/TimeSteppers.jl deleted file mode 100644 index eb7da2f92..000000000 --- a/src/TimeSteppers/TimeSteppers.jl +++ /dev/null @@ -1,7 +0,0 @@ -module TimeSteppers - -using Oceananigans.TimeSteppers: AbstractTimeStepper - -include("positive_conservative_AB2.jl") - -end \ No newline at end of file diff --git a/src/TimeSteppers/positive_conservative_AB2.jl b/src/TimeSteppers/positive_conservative_AB2.jl deleted file mode 100644 index b4260b7ee..000000000 --- a/src/TimeSteppers/positive_conservative_AB2.jl +++ /dev/null @@ -1,162 +0,0 @@ -using Oceananigans.Fields: FunctionField, location -using Oceananigans.TurbulenceClosures: implicit_step! -using Oceananigans.Utils: @apply_regionally, apply_regionally! - -mutable struct PositiveQuasiAdamsBashforth2TimeStepper{FT, GT, IT} <: AbstractTimeStepper - χ :: FT - previous_Δt :: FT - Gⁿ :: GT - G⁻ :: GT - implicit_solver :: IT - intermediate_step :: GT -end - -""" -PositiveQuasiAdamsBashforth2TimeStepper(grid, tracers, - χ = 0.1; - implicit_solver = nothing, - Gⁿ = TendencyFields(grid, tracers), - G⁻ = TendencyFields(grid, tracers)) - -Return a 2nd-order quasi Adams-Bashforth (AB2) time stepper (`QuasiAdamsBashforth2TimeStepper`) -on `grid`, with `tracers`, and AB2 parameter `χ`. The tendency fields `Gⁿ` and `G⁻` can be -specified via optional `kwargs`. - -The 2nd-order quasi Adams-Bashforth timestepper steps forward the state `Uⁿ` by `Δt` via - -```julia -Uⁿ⁺¹ = Uⁿ + Δt * [(3/2 + χ) * Gⁿ - (1/2 + χ) * Gⁿ⁻¹] -``` - -where `Uⁿ` is the state at the ``n``-th timestep, `Gⁿ` is the tendency -at the ``n``-th timestep, and `Gⁿ⁻¹` is the tendency at the previous -timestep (`G⁻`). - -!!! note "First timestep" - For the first timestep, since there are no saved tendencies from the previous timestep, - the `QuasiAdamsBashforth2TimeStepper` performs an Euler timestep: - - ```julia - Uⁿ⁺¹ = Uⁿ + Δt * Gⁿ - ``` -""" -function PositiveQuasiAdamsBashforth2TimeStepper(grid, tracers, - χ = 0.1; - implicit_solver::IT = nothing, - Gⁿ = TendencyFields(grid, tracers), - G⁻ = TendencyFields(grid, tracers), - Cⁿ⁺¹ = TendencyFields(grid, tracers)) where IT - - FT = eltype(grid) - GT = typeof(Gⁿ) - - return PositiveQuasiAdamsBashforth2TimeStepper{FT, GT, IT}(χ, Inf, Gⁿ, G⁻, implicit_solver, Cⁿ⁺¹) -end - -function reset!(timestepper::PositiveQuasiAdamsBashforth2TimeStepper) - timestepper.previous_Δt = Inf - return nothing -end - -##### -##### Time steppping -##### - -""" - time_step!(model::AbstractModel{<:PositiveQuasiAdamsBashforth2TimeStepper}, Δt; euler=false) - -Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and -pressure-correction substep. Setting `euler=true` will take a forward Euler time step. -""" -function time_step!(model::AbstractModel{<:PositiveQuasiAdamsBashforth2TimeStepper}, Δt; callbacks = [], euler=false) - Δt == 0 && @warn "Δt == 0 may cause model blowup!" - - # Shenanigans for properly starting the AB2 loop with an Euler step - euler = euler || (Δt != model.timestepper.previous_Δt) - - χ = ifelse(euler, convert(eltype(model.grid), -0.5), model.timestepper.χ) - - if euler - @debug "Taking a forward Euler step." - # Ensure zeroing out all previous tendency fields to avoid errors in - # case G⁻ includes NaNs. See https://github.com/CliMA/Oceananigans.jl/issues/2259 - for field in model.timestepper.G⁻ - !isnothing(field) && @apply_regionally fill!(field, 0) - end - end - - model.timestepper.previous_Δt = Δt - - # Be paranoid and update state at iteration 0 - model.clock.iteration == 0 && update_state!(model, callbacks) - - @apply_regionally calculate_tendencies!(model, callbacks) - - ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities. - calculate_pressure_correction!(model, Δt) - - @apply_regionally correct_velocities_and_store_tendecies!(model, Δt) - - tick!(model.clock, Δt) - update_state!(model, callbacks) - step_lagrangian_particles!(model, Δt) - - return nothing -end - -function correct_velocities_and_store_tendecies!(model, Δt) - pressure_correct_velocities!(model, Δt) - store_tendencies!(model) - return nothing -end - -##### -##### Time stepping in each step -##### - -""" Generic implementation. """ -function ab2_step!(model, Δt, χ) - - workgroup, worksize = work_layout(model.grid, :xyz) - arch = model.architecture - step_field_kernel! = ab2_step_field!(device(arch), workgroup, worksize) - model_fields = prognostic_fields(model) - - for (i, field) in enumerate(model_fields) - - step_field_kernel!(field, Δt, χ, - model.timestepper.Gⁿ[i], - model.timestepper.G⁻[i]) - - # TODO: function tracer_index(model, field_index) = field_index - 3, etc... - tracer_index = Val(i - 3) # assumption - - implicit_step!(field, - model.timestepper.implicit_solver, - model.closure, - model.diffusivity_fields, - tracer_index, - model.clock, - Δt) - end - - return nothing -end - -""" -Time step velocity fields via the 2nd-order quasi Adams-Bashforth method - - `U^{n+1} = U^n + Δt ((3/2 + χ) * G^{n} - (1/2 + χ) G^{n-1})` - -""" -@kernel function ab2_step_field!(u, Δt, χ, Gⁿ, G⁻) - i, j, k = @index(Global, NTuple) - - T = eltype(u) - one_point_five = convert(T, 1.5) - oh_point_five = convert(T, 0.5) - - @inbounds u[i, j, k] += Δt * ((one_point_five + χ) * Gⁿ[i, j, k] - (oh_point_five + χ) * G⁻[i, j, k]) -end - -@kernel ab2_step_field!(::FunctionField, Δt, χ, Gⁿ, G⁻) = nothing