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 4761ddc56..0f4b34186 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" @@ -15,10 +15,11 @@ 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" +Documenter = "0.27" [extras] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" 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" 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, 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/docs/src/model_implementation.md b/docs/src/model_implementation.md index 59fcd5185..c8e8d1c51 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 sinking_tracers(::NutrientPhytoplankton) = (:P, ) ``` ### Putting it together diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 506befc75..80f28de37 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -20,6 +20,12 @@ 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(κ))) +maximum_numeric_diffusivity(κ::FunctionField) = maximum(κ) const year = years = 365days # just for these idealised cases nothing #hide @@ -48,9 +54,10 @@ 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 + 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 @@ -63,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₂) @@ -73,8 +84,8 @@ biogeochemistry = LOBSTER(; grid, CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = s_function) -model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), +model = NonhydrostaticModel(; grid, clock, + closure = ScalarDiffusivity(ν = κ, κ = κ), biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),)) @@ -106,7 +117,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/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index d1fcc9767..77d898ade 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -7,7 +7,7 @@ using KernelAbstractions using OceanBioME: Biogeochemistry, 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 @@ -35,13 +35,20 @@ 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, model.grid, model.advection, model.tracers, model.timestepper) + 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 + @kernel function store_flat_tendencies!(G⁻, G⁰) i, j = @index(Global, NTuple) @inbounds G⁻[i, j, 1] = G⁰[i, j, 1] @@ -50,10 +57,12 @@ end @inline nitrogen_flux() = 0 @inline carbon_flux() = 0 @inline remineralisation_receiver() = nothing +@inline sinking_tracers() = nothing +@inline required_tracers() = nothing +@inline required_tendencies() = nothing -@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 sinking_advection(bgc, advection) = advection +@inline sinking_advection(bgc, advection::NamedTuple) = advection[sinking_tracers(bgc)] @inline advection_scheme(advection, val_tracer) = advection @inline advection_scheme(advection::NamedTuple, val_tracer::Val{T}) where T = advection[T] @@ -81,7 +90,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) diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 62046ccba..3828b1629 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, @@ -74,18 +74,22 @@ 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(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)) 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, timestepper) +@kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, tendencies, sediment_tendencies) i, j = @index(Global, NTuple) k = bottom_index(i, j, sediment) @@ -97,9 +101,9 @@ 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 timestepper.Gⁿ[remineralisation_receiver(bgc)][i, j, 1] += 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)") diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index 77e804910..c30192a9b 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -25,22 +25,40 @@ 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 """ 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), - depth = abs(znodes(grid, Face())[1]), - solid_dep_params::P4 = (A = 0.233, B = 0.336, C = 982, D = - 1.548, depth = depth)) + 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])) Return a single-layer "multi G" sediment model (`SimpleMultiG`) on `grid`, where parameters can be optionally specified. @@ -71,37 +89,17 @@ 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, - 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), - depth = abs(znodes(grid, Face())[1]), - solid_dep_params::P4 = (A = 0.233, - B = 0.336, - C = 982, - D = - 1.548, - depth = depth)) where {FT, P1, P2, P3, P4} + 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." @@ -113,38 +111,35 @@ 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) - - 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) + bottom_indices = arch_array(architecture(grid), calculate_bottom_indices(grid)) + + 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) = - 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), - adapt(to, sediment.tendencies)) + nothing, + 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, @@ -154,17 +149,20 @@ 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, timestepper) +@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, tendencies, sediment_tendencies) 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()) @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 @@ -179,57 +177,56 @@ 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 - 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, 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] + + 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(k) * 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 + 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.F * log(NO₃) ^ 2) / (Cᵐⁱⁿ * day) - - if isnan(pₐₙₒₓ) - 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] - - 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 + + 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 + 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 (I think) end end diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index badb7610a..53e13c963 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: fractional_indices, indices + struct pCO₂{P0, P1, P2, PB, PW, FT} solubility :: P0 bicarbonate_dissociation :: P1 @@ -41,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²) - @@ -257,9 +262,23 @@ 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 + + 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) -# 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 + _, _, ks = indices(conc) + + 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/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 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 69fdbba6c..422bc2b7e 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 @@ -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 @@ -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 @@ -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)") @@ -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,16 +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(::VariableRedfieldLobster) = (:sPON, :bPON, :sPOC, :bPOC) end # module diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index ed892a72e..1134967c4 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 @@ -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 @@ -286,7 +286,7 @@ 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] + return (u = ZeroField(), v = ZeroField(), w = bgc.sinking_velocities[tracer_name]) else return (u = ZeroField(), v = ZeroField(), w = ZeroField()) end @@ -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, + 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), - npzd.maximum_grazing_rate, - npzd.grazing_half_saturation, - npzd.assimulation_efficiency, - npzd.base_excretion_rate, - npzd.zoo_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), - npzd.remineralization_rate, + adapt(to, npzd.remineralization_rate), adapt(to, npzd.sinking_velocities)) @@ -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(bgc::NPZD) = keys(bgc.sinking_velocities) end # module diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 1d8ead7c3..122337337 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! @@ -200,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), @@ -260,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) @@ -342,18 +343,18 @@ 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, model.grid, - model.velocities, model.tracers, model.clock, Δt) + update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying_biogeochemistry, model.grid, + total_velocities(model), 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 +370,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(light_attenuation).PAR) photo = photosynthesis(T, PAR, p) e = exudation(C, p) @@ -510,10 +511,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, PAR_field) bgc_tracers = required_biogeochemical_tracers(bgc) - auxiliary_fields = biogeochemical_auxiliary_fields(bgc) i, j, k = fractional_indices(x, y, z, (Center(), Center(), Center()), grid) @@ -524,7 +523,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)) diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index a7b872642..ac61ed17c 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 @@ -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, @@ -102,24 +104,24 @@ 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) update_tendencies!(bgc, bgc.particles, model) update_tendencies!(bgc, bgc.modifiers, model) + synchronize(device(architecture(model))) 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) +@inline (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) diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index a49229adc..8cc839ac9 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 @@ -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 @@ -77,42 +77,51 @@ 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 = NamedTuple{tracers}(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 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.") + + t < 0 && (t = NaN) + @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 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) diff --git a/src/Utils/timestep.jl b/src/Utils/timestep.jl index d9d6096e7..f6aa7db6e 100644 --- a/src/Utils/timestep.jl +++ b/src/Utils/timestep.jl @@ -1,21 +1,6 @@ using Oceananigans.Advection: cell_advection_timescale using Oceananigans.Grids: Center, znodes, zspacing, minimum_zspacing -@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/runtests.jl b/test/runtests.jl index 8dff45c8d..6a721a0e4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,6 @@ -using OceanBioME, Documenter, Test +using OceanBioME, Documenter, Test, CUDA, Oceananigans + +architecture = CUDA.has_cuda() ? GPU() : CPU() include("test_utils.jl") include("test_light.jl") @@ -8,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 diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index db25b6d49..df86a1e09 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(architecture)) with carbonates $(carbonates ? :✅ : :❌), oxygen $(oxygen ? :✅ : :❌), variable redfield $(variable_redfield ? :✅ : :❌), sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" + @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 end \ No newline at end of file diff --git a/test/test_NPZD.jl b/test/test_NPZD.jl index 465d41266..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,27 +34,32 @@ 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 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(architecture)) with sinking $(sinking ? :✅ : :❌), open bottom $(open_bottom ? :✅ : :❌))" + @testset "NPZD $architecture, $sinking, $open_bottom" begin + test_NPZD(grid, sinking, open_bottom) end end end diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl index 0c3dcad08..08717ecce 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", @@ -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)) @@ -51,12 +53,16 @@ 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) +#Oceananigans.fill_halo_regions!(conc_field) + @testset "Gas exchange coupling" begin - grid = RectilinearGrid(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] + for air_concentration in [413.1, conc_function]#, conc_field] @info "Testing with $(typeof(air_concentration))" test_gas_exchange_model(grid, air_concentration) end diff --git a/test/test_light.jl b/test/test_light.jl index fdcf99169..3b274be3e 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -37,18 +37,16 @@ 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 -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..1a2214bf4 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -9,11 +9,11 @@ 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) - copyto!(intercepted_tendencies[tracer], model.timestepper.Gⁿ[tracer]) + intercepted_tendencies[tracer] = Array(interior(model.timestepper.Gⁿ[tracer])) end end @@ -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, @@ -52,9 +52,26 @@ total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) + 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(::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, @@ -74,29 +91,30 @@ 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) - 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) # 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) - @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) + 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 - @test N₁ ≈ N₀ + rtol = ifelse(isa(architecture, CPU), max(√eps(N₀), √eps(N₁)), 5e-7) + @test isapprox(N₀, N₁; rtol) return nothing end @@ -113,30 +131,29 @@ 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)) 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 diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index ada407cd9..9cd4bdc48 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -2,23 +2,39 @@ 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) - 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, 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) .+ Array(interior(tracers.Z))) .* redfield + @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 - 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) @@ -36,13 +52,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_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_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_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, 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 @@ -50,29 +64,33 @@ 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, 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 @test initial_kelp_N != final_kelp_N @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) - intercepted_tendencies = TracerFields(keys(model.tracers), grid) + # slow but easiest to have this just done on CPU + 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) 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[idx] .!= Array(interior(model.timestepper.Gⁿ[tracer]))) for (idx, tracer) in enumerate((:NO₃, :NH₄, :DIC, :DOC, :bPON, :bPOC))]) end diff --git a/test/test_utils.jl b/test/test_utils.jl index 86586a791..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)) @@ -28,7 +11,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,13 +28,14 @@ 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 - for arch in (CPU(), ) - @test test_column_diffusion_timescale(arch) - @test test_negative_scaling(arch) - @test test_negative_zeroing(arch) - end + @test test_negative_scaling(architecture) + @test test_negative_zeroing(architecture) end