diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6dbdee3a..ff53d3dd 100755 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,7 +22,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -47,7 +47,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -72,7 +72,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -97,7 +97,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -122,7 +122,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -147,7 +147,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -172,7 +172,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -197,7 +197,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -222,7 +222,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -247,7 +247,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: @@ -272,7 +272,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.12' os: - ubuntu-latest arch: diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 04bad146..e92af96a 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -19,7 +19,7 @@ jobs: - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@v2 with: - version: '1.10' + version: '1.12' - name: Install dependencies run: julia --color=yes --project -e 'using Pkg; Pkg.instantiate()' - name: Build and deploy diff --git a/Project.toml b/Project.toml index 85342022..df260ad4 100755 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Oceanostics" uuid = "d0ccf422-c8fb-49b5-a76d-74acdde946ac" +version = "0.16.5" authors = ["tomchor "] -version = "0.16.4" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -11,7 +11,7 @@ SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" [compat] DocStringExtensions = "0.9" -Oceananigans = "^0.100.3" +Oceananigans = "0.101" SeawaterPolynomials = "0.3" julia = "1.9" diff --git a/src/KineticEnergyEquation.jl b/src/KineticEnergyEquation.jl index 016a197d..3683fac7 100644 --- a/src/KineticEnergyEquation.jl +++ b/src/KineticEnergyEquation.jl @@ -87,11 +87,11 @@ KineticEnergy(model; kwargs...) = KineticEnergy(model, model.velocities...; kwar velocities, tracers, auxiliary_fields, - diffusivity_fields, + closure_fields, pHY′, clock, forcings) - common_args = (buoyancy, background_fields, velocities, tracers, auxiliary_fields, diffusivity_fields, pHY′, clock) + common_args = (buoyancy, background_fields, velocities, tracers, auxiliary_fields, closure_fields, pHY′, clock) u∂ₜu = ℑxᶜᵃᵃ(i, j, k, grid, ψf, velocities.u, u_velocity_tendency, advection, coriolis, stokes_drift, closure, u_immersed_bc, common_args..., forcings.u) v∂ₜv = ℑyᵃᶜᵃ(i, j, k, grid, ψf, velocities.v, v_velocity_tendency, advection, coriolis, stokes_drift, closure, v_immersed_bc, common_args..., forcings.v) w∂ₜw = ℑzᵃᵃᶜ(i, j, k, grid, ψf, velocities.w, w_velocity_tendency, advection, coriolis, stokes_drift, closure, w_immersed_bc, common_args..., forcings.w) @@ -138,7 +138,7 @@ function KineticEnergyTendency(model::NonhydrostaticModel; location = (Center, C model.velocities, model.tracers, model.auxiliary_fields, - model.diffusivity_fields, + model.closure_fields, model.pressures.pHY′, model.clock, model.forcing,) @@ -191,14 +191,14 @@ end #+++ KineticEnergyStress @inline function uᵢ∂ⱼ_τᵢⱼᶜᶜᶜ(i, j, k, grid, closure, - diffusivity_fields, + closure_fields, clock, model_fields, buoyancy) - u∂ⱼ_τ₁ⱼ = ℑxᶜᵃᵃ(i, j, k, grid, ψf, model_fields.u, ∂ⱼ_τ₁ⱼ, closure, diffusivity_fields, clock, model_fields, buoyancy) - v∂ⱼ_τ₂ⱼ = ℑyᵃᶜᵃ(i, j, k, grid, ψf, model_fields.v, ∂ⱼ_τ₂ⱼ, closure, diffusivity_fields, clock, model_fields, buoyancy) - w∂ⱼ_τ₃ⱼ = ℑzᵃᵃᶜ(i, j, k, grid, ψf, model_fields.w, ∂ⱼ_τ₃ⱼ, closure, diffusivity_fields, clock, model_fields, buoyancy) + u∂ⱼ_τ₁ⱼ = ℑxᶜᵃᵃ(i, j, k, grid, ψf, model_fields.u, ∂ⱼ_τ₁ⱼ, closure, closure_fields, clock, model_fields, buoyancy) + v∂ⱼ_τ₂ⱼ = ℑyᵃᶜᵃ(i, j, k, grid, ψf, model_fields.v, ∂ⱼ_τ₂ⱼ, closure, closure_fields, clock, model_fields, buoyancy) + w∂ⱼ_τ₃ⱼ = ℑzᵃᵃᶜ(i, j, k, grid, ψf, model_fields.w, ∂ⱼ_τ₃ⱼ, closure, closure_fields, clock, model_fields, buoyancy) return u∂ⱼ_τ₁ⱼ+ v∂ⱼ_τ₂ⱼ + w∂ⱼ_τ₃ⱼ end @@ -240,7 +240,7 @@ function KineticEnergyStress(model; location = (Center, Center, Center)) model_fields = (; model_fields..., w=ZeroField()) end dependencies = (model.closure, - model.diffusivity_fields, + model.closure_fields, model.clock, fields(model), model.buoyancy) @@ -434,18 +434,18 @@ Axᶠᶜᶠ_δwᶠᶜᶠ_F₃₁ᶠᶜᶠ(i, j, k, grid, closure, K_fields, clo, Ayᶜᶠᶠ_δwᶜᶠᶠ_F₃₂ᶜᶠᶠ(i, j, k, grid, closure, K_fields, clo, fields, b) = -Ayᶜᶠᶠ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, grid, fields.w) * viscous_flux_wy(i, j, k, grid, closure, K_fields, clo, fields, b) Azᶜᶜᶜ_δwᶜᶜᶜ_F₃₃ᶜᶜᶜ(i, j, k, grid, closure, K_fields, clo, fields, b) = -Azᶜᶜᶜ(i, j, k, grid) * δzᵃᵃᶜ(i, j, k, grid, fields.w) * viscous_flux_wz(i, j, k, grid, closure, K_fields, clo, fields, b) -@inline viscous_dissipation_rate_ccc(i, j, k, grid, diffusivity_fields, fields, p) = - (Axᶜᶜᶜ_δuᶜᶜᶜ_F₁₁ᶜᶜᶜ(i, j, k, grid, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # C, C, C - ℑxyᶜᶜᵃ(i, j, k, grid, Ayᶠᶠᶜ_δuᶠᶠᶜ_F₁₂ᶠᶠᶜ, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # F, F, C → C, C, C - ℑxzᶜᵃᶜ(i, j, k, grid, Azᶠᶜᶠ_δuᶠᶜᶠ_F₁₃ᶠᶜᶠ, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # F, C, F → C, C, C +@inline viscous_dissipation_rate_ccc(i, j, k, grid, closure_fields, fields, p) = + (Axᶜᶜᶜ_δuᶜᶜᶜ_F₁₁ᶜᶜᶜ(i, j, k, grid, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # C, C, C + ℑxyᶜᶜᵃ(i, j, k, grid, Ayᶠᶠᶜ_δuᶠᶠᶜ_F₁₂ᶠᶠᶜ, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # F, F, C → C, C, C + ℑxzᶜᵃᶜ(i, j, k, grid, Azᶠᶜᶠ_δuᶠᶜᶠ_F₁₃ᶠᶜᶠ, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # F, C, F → C, C, C - ℑxyᶜᶜᵃ(i, j, k, grid, Axᶠᶠᶜ_δvᶠᶠᶜ_F₂₁ᶠᶠᶜ, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # F, F, C → C, C, C - Ayᶜᶜᶜ_δvᶜᶜᶜ_F₂₂ᶜᶜᶜ(i, j, k, grid, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # C, C, C - ℑyzᵃᶜᶜ(i, j, k, grid, Azᶜᶠᶠ_δvᶜᶠᶠ_F₂₃ᶜᶠᶠ, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # C, F, F → C, C, C + ℑxyᶜᶜᵃ(i, j, k, grid, Axᶠᶠᶜ_δvᶠᶠᶜ_F₂₁ᶠᶠᶜ, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # F, F, C → C, C, C + Ayᶜᶜᶜ_δvᶜᶜᶜ_F₂₂ᶜᶜᶜ(i, j, k, grid, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # C, C, C + ℑyzᵃᶜᶜ(i, j, k, grid, Azᶜᶠᶠ_δvᶜᶠᶠ_F₂₃ᶜᶠᶠ, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # C, F, F → C, C, C - ℑxzᶜᵃᶜ(i, j, k, grid, Axᶠᶜᶠ_δwᶠᶜᶠ_F₃₁ᶠᶜᶠ, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # F, C, F → C, C, C - ℑyzᵃᶜᶜ(i, j, k, grid, Ayᶜᶠᶠ_δwᶜᶠᶠ_F₃₂ᶜᶠᶠ, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) + # C, F, F → C, C, C - Azᶜᶜᶜ_δwᶜᶜᶜ_F₃₃ᶜᶜᶜ(i, j, k, grid, p.closure, diffusivity_fields, p.clock, fields, p.buoyancy) # C, C, C + ℑxzᶜᵃᶜ(i, j, k, grid, Axᶠᶜᶠ_δwᶠᶜᶠ_F₃₁ᶠᶜᶠ, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # F, C, F → C, C, C + ℑyzᵃᶜᶜ(i, j, k, grid, Ayᶜᶠᶠ_δwᶜᶠᶠ_F₃₂ᶜᶠᶠ, p.closure, closure_fields, p.clock, fields, p.buoyancy) + # C, F, F → C, C, C + Azᶜᶜᶜ_δwᶜᶜᶜ_F₃₃ᶜᶜᶜ(i, j, k, grid, p.closure, closure_fields, p.clock, fields, p.buoyancy) # C, C, C ) / Vᶜᶜᶜ(i, j, k, grid) # This division by volume, coupled with the call to A*δuᵢ above, ensures a derivative operation const KineticEnergyDissipationRate = CustomKFO{<:typeof(viscous_dissipation_rate_ccc)} @@ -484,7 +484,7 @@ function DissipationRate(model; U=ZeroField(), V=ZeroField(), W=ZeroField(), model.buoyancy) return KernelFunctionOperation{Center, Center, Center}(viscous_dissipation_rate_ccc, model.grid, - model.diffusivity_fields, model_fields, parameters) + model.closure_fields, model_fields, parameters) end #--- @@ -499,7 +499,7 @@ end Σˣᶻ² = ℑxzᶜᵃᶜ(i, j, k, grid, fψ_plus_gφ², ∂zᶠᶜᶠ, u, ∂xᶠᶜᶠ, w) / 4 Σʸᶻ² = ℑyzᵃᶜᶜ(i, j, k, grid, fψ_plus_gφ², ∂zᶜᶠᶠ, v, ∂yᶜᶠᶠ, w) / 4 - ν = _νᶜᶜᶜ(i, j, k, grid, p.closure, p.diffusivity_fields, p.clock) + ν = _νᶜᶜᶜ(i, j, k, grid, p.closure, p.closure_fields, p.clock, p.model_fields) return 2ν * (Σˣˣ² + Σʸʸ² + Σᶻᶻ² + 2 * (Σˣʸ² + Σˣᶻ² + Σʸᶻ²)) end @@ -531,17 +531,17 @@ KernelFunctionOperation at (Center, Center, Center) └── arguments: ("Field", "Field", "Field", "NamedTuple") ``` """ -function KineticEnergyIsotropicDissipationRate(u, v, w, closure, diffusivity_fields, clock; location = (Center, Center, Center)) +function KineticEnergyIsotropicDissipationRate(u, v, w, closure, closure_fields, model_fields, clock; location = (Center, Center, Center)) validate_location(location, "KineticEnergyIsotropicDissipationRate") validate_dissipative_closure(closure) - parameters = (; closure, diffusivity_fields, clock) + parameters = (; closure, closure_fields, clock, model_fields) return KernelFunctionOperation{Center, Center, Center}(isotropic_viscous_dissipation_rate_ccc, u.grid, u, v, w, parameters) end @inline KineticEnergyIsotropicDissipationRate(model; location = (Center, Center, Center)) = - KineticEnergyIsotropicDissipationRate(model.velocities..., model.closure, model.diffusivity_fields, model.clock; location = location) + KineticEnergyIsotropicDissipationRate(model.velocities..., model.closure, model.closure_fields, fields(model), model.clock; location = location) #--- end # module diff --git a/src/ProgressMessengers/cfl.jl b/src/ProgressMessengers/cfl.jl index fea9c493..3f231425 100644 --- a/src/ProgressMessengers/cfl.jl +++ b/src/ProgressMessengers/cfl.jl @@ -41,7 +41,7 @@ tuple_to_op(::Nothing) = nothing tuple_to_op(ν_tuple::Tuple) = sum(ν_tuple) @inline function (maxν::MaxViscosity)(simulation) - ν = tuple_to_op(viscosity(simulation.model.closure, simulation.model.diffusivity_fields)) + ν = tuple_to_op(viscosity(simulation.model.closure, simulation.model.closure_fields)) ν_max = maximum(abs, ν) message = @sprintf("%.2g", ν_max) maxν.with_prefix && (message = "νₘₐₓ = " * message) diff --git a/src/TracerEquation.jl b/src/TracerEquation.jl index 66f8e0ef..980a0e2a 100755 --- a/src/TracerEquation.jl +++ b/src/TracerEquation.jl @@ -11,9 +11,9 @@ export Advection, Diffusion, ImmersedDiffusion, TotalDiffusion, Forcing, TracerAdvection, TracerDiffusion, TracerImmersedDiffusion, TracerTotalDiffusion, TracerForcing # Inline function for total diffusion -@inline total_∇_dot_qᶜ(i, j, k, grid, c, c_immersed_bc, closure, diffusivity_fields, val_tracer_index, clock, model_fields, buoyancy) = - ∇_dot_qᶜ(i, j, k, grid, closure, diffusivity_fields, val_tracer_index, c, clock, model_fields, buoyancy) + - immersed_∇_dot_qᶜ(i, j, k, grid, c, c_immersed_bc, closure, diffusivity_fields, val_tracer_index, clock, model_fields, buoyancy) +@inline total_∇_dot_qᶜ(i, j, k, grid, c, c_immersed_bc, closure, closure_fields, val_tracer_index, clock, model_fields, buoyancy) = + ∇_dot_qᶜ(i, j, k, grid, closure, closure_fields, val_tracer_index, c, clock, model_fields, buoyancy) + + immersed_∇_dot_qᶜ(i, j, k, grid, c, c_immersed_bc, closure, closure_fields, val_tracer_index, clock, model_fields, buoyancy) # Type aliases for major functions const Advection = CustomKFO{<:typeof(div_Uc)} @@ -93,15 +93,15 @@ KernelFunctionOperation at (Center, Center, Center) └── arguments: ("Nothing", "Nothing", "Val", "Field", "Clock", "NamedTuple", "Nothing") ``` """ -function Diffusion(model, val_tracer_index, c, closure, diffusivity_fields, clock, model_fields, buoyancy; location = (Center, Center, Center)) +function Diffusion(model, val_tracer_index, c, closure, closure_fields, clock, model_fields, buoyancy; location = (Center, Center, Center)) validate_location(location, "Diffusion", (Center, Center, Center)) - return KernelFunctionOperation{Center, Center, Center}(∇_dot_qᶜ, model.grid, closure, diffusivity_fields, val_tracer_index, c, clock, model_fields, buoyancy) + return KernelFunctionOperation{Center, Center, Center}(∇_dot_qᶜ, model.grid, closure, closure_fields, val_tracer_index, c, clock, model_fields, buoyancy) end function Diffusion(model, tracer_name; kwargs...) tracer_index = findfirst(x -> x == tracer_name, keys(model.tracers)) @inbounds c = model.tracers[tracer_name] - return Diffusion(model, Val(tracer_index), c, model.closure, model.diffusivity_fields, model.clock, fields(model), model.buoyancy; kwargs...) + return Diffusion(model, Val(tracer_index), c, model.closure, model.closure_fields, model.clock, fields(model), model.buoyancy; kwargs...) end @@ -129,16 +129,16 @@ KernelFunctionOperation at (Center, Center, Center) └── arguments: ("Field", "Nothing", "Nothing", "Nothing", "Val", "Clock", "NamedTuple") ``` """ -function ImmersedDiffusion(model, c, c_immersed_bc, closure, diffusivity_fields, val_tracer_index, clock, model_fields; location = (Center, Center, Center)) +function ImmersedDiffusion(model, c, c_immersed_bc, closure, closure_fields, val_tracer_index, clock, model_fields; location = (Center, Center, Center)) validate_location(location, "ImmersedDiffusion", (Center, Center, Center)) - return KernelFunctionOperation{Center, Center, Center}(immersed_∇_dot_qᶜ, model.grid, c, c_immersed_bc, closure, diffusivity_fields, val_tracer_index, clock, model_fields) + return KernelFunctionOperation{Center, Center, Center}(immersed_∇_dot_qᶜ, model.grid, c, c_immersed_bc, closure, closure_fields, val_tracer_index, clock, model_fields) end function ImmersedDiffusion(model, tracer_name; kwargs...) tracer_index = findfirst(x -> x == tracer_name, keys(model.tracers)) tracer = model.tracers[tracer_name] immersed_bc = tracer.boundary_conditions.immersed - return ImmersedDiffusion(model, tracer, immersed_bc, model.closure, model.diffusivity_fields, Val(tracer_index), model.clock, fields(model); kwargs...) + return ImmersedDiffusion(model, tracer, immersed_bc, model.closure, model.closure_fields, Val(tracer_index), model.clock, fields(model); kwargs...) end """ @@ -165,16 +165,16 @@ KernelFunctionOperation at (Center, Center, Center) └── arguments: ("Field", "Nothing", "Nothing", "Nothing", "Val", "Clock", "NamedTuple", "Nothing") ``` """ -function TotalDiffusion(model, c, c_immersed_bc, closure, diffusivity_fields, val_tracer_index, clock, model_fields, buoyancy; location = (Center, Center, Center)) +function TotalDiffusion(model, c, c_immersed_bc, closure, closure_fields, val_tracer_index, clock, model_fields, buoyancy; location = (Center, Center, Center)) validate_location(location, "TotalDiffusion", (Center, Center, Center)) - return KernelFunctionOperation{Center, Center, Center}(total_∇_dot_qᶜ, model.grid, c, c_immersed_bc, closure, diffusivity_fields, val_tracer_index, clock, model_fields, buoyancy) + return KernelFunctionOperation{Center, Center, Center}(total_∇_dot_qᶜ, model.grid, c, c_immersed_bc, closure, closure_fields, val_tracer_index, clock, model_fields, buoyancy) end function TotalDiffusion(model, tracer_name; kwargs...) tracer_index = findfirst(x -> x == tracer_name, keys(model.tracers)) tracer = model.tracers[tracer_index] immersed_bc = tracer.boundary_conditions.immersed - return TotalDiffusion(model, tracer, immersed_bc, model.closure, model.diffusivity_fields, Val(tracer_index), model.clock, fields(model), model.buoyancy; kwargs...) + return TotalDiffusion(model, tracer, immersed_bc, model.closure, model.closure_fields, Val(tracer_index), model.clock, fields(model), model.buoyancy; kwargs...) end #--- diff --git a/src/TracerVarianceEquation.jl b/src/TracerVarianceEquation.jl index 5da13ee2..e7bb6baf 100644 --- a/src/TracerVarianceEquation.jl +++ b/src/TracerVarianceEquation.jl @@ -80,7 +80,7 @@ function TracerVarianceTendency(model::NonhydrostaticModel, tracer_name; locatio model.velocities, model.tracers, model.auxiliary_fields, - model.diffusivity_fields, + model.closure_fields, model.clock, model.forcing[tracer_name]) @@ -121,7 +121,7 @@ julia> DIFF = TracerVarianceEquation.TracerVarianceDiffusion(model, :b) KernelFunctionOperation at (Center, Center, Center) ├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── kernel_function: c∇_dot_qᶜ (generic function with 1 method) -└── arguments: ("Smagorinsky", "NamedTuple", "Val", "Field", "Clock", "NamedTuple", "Nothing") +└── arguments: ("Oceananigans.TurbulenceClosures.Smagorinskys.Smagorinsky", "NamedTuple", "Val", "Field", "Clock", "NamedTuple", "Nothing") ``` """ function TracerVarianceDiffusion(model, tracer_name; location = (Center, Center, Center)) @@ -129,7 +129,7 @@ function TracerVarianceDiffusion(model, tracer_name; location = (Center, Center, tracer_index = findfirst(n -> n === tracer_name, propertynames(model.tracers)) dependencies = (model.closure, - model.diffusivity_fields, + model.closure_fields, Val(tracer_index), model.tracers[tracer_name], model.clock, @@ -141,16 +141,16 @@ end #+++ Tracer variance dissipation rate # Variance dissipation rate at fcc -@inline Axᶠᶜᶜ_δcᶠᶜᶜ_q₁ᶠᶜᶜ(i, j, k, grid, closure, diffusivity_fields, id, c, args...) = - - Axᶠᶜᶜ(i, j, k, grid) * δxᶠᵃᵃ(i, j, k, grid, c) * diffusive_flux_x(i, j, k, grid, closure, diffusivity_fields, id, c, args...) +@inline Axᶠᶜᶜ_δcᶠᶜᶜ_q₁ᶠᶜᶜ(i, j, k, grid, closure, closure_fields, id, c, args...) = + - Axᶠᶜᶜ(i, j, k, grid) * δxᶠᵃᵃ(i, j, k, grid, c) * diffusive_flux_x(i, j, k, grid, closure, closure_fields, id, c, args...) # Variance dissipation rate at cfc -@inline Ayᶜᶠᶜ_δcᶜᶠᶜ_q₂ᶜᶠᶜ(i, j, k, grid, closure, diffusivity_fields, id, c, args...) = - - Ayᶜᶠᶜ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, grid, c) * diffusive_flux_y(i, j, k, grid, closure, diffusivity_fields, id, c, args...) +@inline Ayᶜᶠᶜ_δcᶜᶠᶜ_q₂ᶜᶠᶜ(i, j, k, grid, closure, closure_fields, id, c, args...) = + - Ayᶜᶠᶜ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, grid, c) * diffusive_flux_y(i, j, k, grid, closure, closure_fields, id, c, args...) # Variance dissipation rate at ccf -@inline Azᶜᶜᶠ_δcᶜᶜᶠ_q₃ᶜᶜᶠ(i, j, k, grid, closure, diffusivity_fields, id, c, args...) = - - Azᶜᶜᶠ(i, j, k, grid) * δzᵃᵃᶠ(i, j, k, grid, c) * diffusive_flux_z(i, j, k, grid, closure, diffusivity_fields, id, c, args...) +@inline Azᶜᶜᶠ_δcᶜᶜᶠ_q₃ᶜᶜᶠ(i, j, k, grid, closure, closure_fields, id, c, args...) = + - Azᶜᶜᶠ(i, j, k, grid) * δzᵃᵃᶠ(i, j, k, grid, c) * diffusive_flux_z(i, j, k, grid, closure, closure_fields, id, c, args...) @inline tracer_variance_dissipation_rate_ccc(i, j, k, grid, args...) = 2 * (ℑxᶜᵃᵃ(i, j, k, grid, Axᶠᶜᶜ_δcᶠᶜᶜ_q₁ᶠᶜᶜ, args...) + # F, C, C → C, C, C @@ -190,7 +190,7 @@ julia> χ = TracerVarianceEquation.TracerVarianceDissipationRate(model, :b) KernelFunctionOperation at (Center, Center, Center) ├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── kernel_function: tracer_variance_dissipation_rate_ccc (generic function with 1 method) -└── arguments: ("Smagorinsky", "NamedTuple", "Val", "Field", "Clock", "NamedTuple", "Nothing") +└── arguments: ("Oceananigans.TurbulenceClosures.Smagorinskys.Smagorinsky", "NamedTuple", "Val", "Field", "Clock", "NamedTuple", "Nothing") julia> b̄ = Field(Average(model.tracers.b, dims=(1,2))); @@ -200,7 +200,7 @@ julia> χb = TracerVarianceEquation.TracerVarianceDissipationRate(model, :b, tra KernelFunctionOperation at (Center, Center, Center) ├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── kernel_function: tracer_variance_dissipation_rate_ccc (generic function with 1 method) -└── arguments: ("Smagorinsky", "NamedTuple", "Val", "Oceananigans.AbstractOperations.BinaryOperation", "Clock", "NamedTuple", "Nothing") +└── arguments: ("Oceananigans.TurbulenceClosures.Smagorinskys.Smagorinsky", "NamedTuple", "Val", "Oceananigans.AbstractOperations.BinaryOperation", "Clock", "NamedTuple", "Nothing") ``` """ function TracerVarianceDissipationRate(model, tracer_name; tracer = nothing, location = (Center, Center, Center)) @@ -213,7 +213,7 @@ function TracerVarianceDissipationRate(model, tracer_name; tracer = nothing, loc tracer = tracer == nothing ? model.tracers[tracer_name] : tracer return KernelFunctionOperation{Center, Center, Center}(tracer_variance_dissipation_rate_ccc, model.grid, model.closure, - model.diffusivity_fields, + model.closure_fields, Val(tracer_index), tracer, model.clock, diff --git a/src/TurbulentKineticEnergyEquation.jl b/src/TurbulentKineticEnergyEquation.jl index 0f04c47a..5d062322 100644 --- a/src/TurbulentKineticEnergyEquation.jl +++ b/src/TurbulentKineticEnergyEquation.jl @@ -13,6 +13,7 @@ using Oceananigans.AbstractOperations using Oceananigans.AbstractOperations: KernelFunctionOperation using Oceananigans.Grids: Center using Oceananigans.Fields: ZeroField +using Oceananigans: fields using Oceanostics: validate_location, CustomKFO using Oceanostics.KineticEnergyEquation: KineticEnergyIsotropicDissipationRate @@ -88,7 +89,7 @@ KernelFunctionOperation at (Center, Center, Center) @inline function TurbulentKineticEnergyIsotropicDissipationRate(model; U=ZeroField(), V=ZeroField(), W=ZeroField(), kwargs...) u, v, w = model.velocities - return TurbulentKineticEnergyIsotropicDissipationRate((u - U), (v - V), (w - W), model.closure, model.diffusivity_fields, model.clock; kwargs...) + return TurbulentKineticEnergyIsotropicDissipationRate((u - U), (v - V), (w - W), model.closure, model.closure_fields, fields(model), model.clock; kwargs...) end const IsotropicDissipationRate = TurbulentKineticEnergyIsotropicDissipationRate diff --git a/test/test_budgets.jl b/test/test_budgets.jl index 481b7963..d18fb3a1 100644 --- a/test/test_budgets.jl +++ b/test/test_budgets.jl @@ -50,7 +50,7 @@ function test_tracer_variance_budget(; arch, N=16, rtol=0.01, stop_time=0.1, clo c.data.parent .-= mean(c) u.data.parent .-= mean(u) - κ = diffusivity(model.closure, model.diffusivity_fields, Val(:c)) + κ = diffusivity(model.closure, model.closure_fields, Val(:c)) κ = κ isa Tuple ? Field(sum(κ)) : κ Δt = min(minimum_zspacing(grid)^2/maximum(κ)/10, minimum_zspacing(grid)/maximum(u) / 10) simulation = Simulation(model; Δt, stop_time) diff --git a/test/test_canonical_flows.jl b/test/test_canonical_flows.jl index 3ad6cff8..137cf522 100644 --- a/test/test_canonical_flows.jl +++ b/test/test_canonical_flows.jl @@ -2,7 +2,7 @@ using Test using CUDA: has_cuda_gpu, @allowscalar using Oceananigans -using Oceananigans.TurbulenceClosures.Smagorinskys: LagrangianAveraging +using Oceananigans.TurbulenceClosures.Smagorinskys: LagrangianAveraging, DynamicCoefficient using Oceanostics @@ -25,8 +25,8 @@ grid_noise(x, y, z) = randn() #+++ Test options closures = (ScalarDiffusivity(ν=1e-6, κ=1e-7), SmagorinskyLilly(), - Smagorinsky(coefficient=DynamicCoefficient(averaging=(1, 2))), - Smagorinsky(coefficient=DynamicCoefficient(averaging=LagrangianAveraging())), + DynamicSmagorinsky(averaging=(1, 2)), + DynamicSmagorinsky(averaging=LagrangianAveraging()), (ScalarDiffusivity(ν=1e-6, κ=1e-7), AnisotropicMinimumDissipation()),) grids = Dict("regular grid" => regular_grid, @@ -54,9 +54,9 @@ function test_uniform_strain_flow(grid; model_type=NonhydrostaticModel, closure= idxs = (model.grid.Nx÷2, model.grid.Ny÷2, model.grid.Nz÷2) # Get a value far from boundaries if model.closure isa Tuple - ν_field = Field(sum(viscosity(model.closure, model.diffusivity_fields))) + ν_field = Field(sum(viscosity(model.closure, model.closure_fields))) else - ν_field = viscosity(model.closure, model.diffusivity_fields) + ν_field = viscosity(model.closure, model.closure_fields) end @allowscalar begin @@ -91,9 +91,9 @@ function test_solid_body_rotation_flow(grid; model_type=NonhydrostaticModel, clo idxs = (model.grid.Nx÷2, model.grid.Ny÷2, model.grid.Nz÷2) # Get a value far from boundaries if model.closure isa Tuple - ν_field = Field(sum(viscosity(model.closure, model.diffusivity_fields))) + ν_field = Field(sum(viscosity(model.closure, model.closure_fields))) else - ν_field = viscosity(model.closure, model.diffusivity_fields) + ν_field = viscosity(model.closure, model.closure_fields) end @allowscalar begin @@ -125,9 +125,9 @@ function test_uniform_shear_flow(grid; model_type=NonhydrostaticModel, closure=S idxs = (model.grid.Nx÷2, model.grid.Ny÷2, model.grid.Nz÷2) # Get a value far from boundaries if model.closure isa Tuple - ν_field = Field(sum(viscosity(model.closure, model.diffusivity_fields))) + ν_field = Field(sum(viscosity(model.closure, model.closure_fields))) else - ν_field = viscosity(model.closure, model.diffusivity_fields) + ν_field = viscosity(model.closure, model.closure_fields) end @allowscalar begin diff --git a/test/test_kinetic_energy_equation.jl b/test/test_kinetic_energy_equation.jl index f235a039..651c7ea4 100644 --- a/test/test_kinetic_energy_equation.jl +++ b/test/test_kinetic_energy_equation.jl @@ -72,9 +72,9 @@ function test_ke_dissipation_rate_terms(grid; model_type=NonhydrostaticModel, cl ε_iso1_field = Field(ε_iso1) @test ε_iso1 isa KineticEnergyEquation.KineticEnergyIsotropicDissipationRate - # Test the full signature: KineticEnergyIsotropicDissipationRate(u, v, w, closure, diffusivity_fields, clock; location) + # Test the full signature: KineticEnergyIsotropicDissipationRate(u, v, w, closure, closure_fields, fields(model), clock; location) u, v, w = model.velocities - ε_iso2 = KineticEnergyEquation.KineticEnergyIsotropicDissipationRate(u, v, w, model.closure, model.diffusivity_fields, model.clock) + ε_iso2 = KineticEnergyEquation.KineticEnergyIsotropicDissipationRate(u, v, w, model.closure, model.closure_fields, fields(model), model.clock) ε_iso2_field = Field(ε_iso2) @test ε_iso2 isa KineticEnergyEquation.KineticEnergyIsotropicDissipationRate @@ -83,7 +83,7 @@ function test_ke_dissipation_rate_terms(grid; model_type=NonhydrostaticModel, cl # Test with different location parameters ε_iso1_ccc = KineticEnergyEquation.IsotropicDissipationRate(model; location = (Center, Center, Center)) - ε_iso2_ccc = KineticEnergyEquation.IsotropicDissipationRate(u, v, w, model.closure, model.diffusivity_fields, model.clock; location = (Center, Center, Center)) + ε_iso2_ccc = KineticEnergyEquation.IsotropicDissipationRate(u, v, w, model.closure, model.closure_fields, fields(model), model.clock; location = (Center, Center, Center)) @test ε_iso1_ccc isa KineticEnergyEquation.IsotropicDissipationRate @test ε_iso2_ccc isa KineticEnergyEquation.IsotropicDissipationRate @test all(interior(Field(ε_iso1_ccc)) .≈ interior(Field(ε_iso2_ccc))) diff --git a/test/test_progress_messengers.jl b/test/test_progress_messengers.jl index beae8515..d687de03 100644 --- a/test/test_progress_messengers.jl +++ b/test/test_progress_messengers.jl @@ -3,6 +3,7 @@ using CUDA: has_cuda_gpu using Oceananigans using Oceananigans.TurbulenceClosures.Smagorinskys: LagrangianAveraging +using Oceananigans.TurbulenceClosures: DynamicCoefficient, LagrangianAveraging using Oceanostics using Oceanostics.ProgressMessengers @@ -16,8 +17,8 @@ regular_grid = RectilinearGrid(arch, size=(N, N, N), extent=(1, 1, 1)) #+++ Testing options closures = (ScalarDiffusivity(ν=1e-6, κ=1e-7), SmagorinskyLilly(), - Smagorinsky(coefficient=DynamicCoefficient(averaging=(1, 2))), - Smagorinsky(coefficient=DynamicCoefficient(averaging=LagrangianAveraging())), + DynamicSmagorinsky(averaging=(1, 2)), + DynamicSmagorinsky(averaging=LagrangianAveraging()), (ScalarDiffusivity(ν=1e-6, κ=1e-7), AnisotropicMinimumDissipation()),) #--- diff --git a/test/test_tracer_diagnostics.jl b/test/test_tracer_diagnostics.jl index 5ed1d3c5..1785773d 100644 --- a/test/test_tracer_diagnostics.jl +++ b/test/test_tracer_diagnostics.jl @@ -62,7 +62,7 @@ function test_tracer_terms(model) ADV_field = Field(ADV) @test ADV_field isa Field - DIFF = TracerEquation.TracerDiffusion(model, :a, model.tracers.a, model.closure, model.diffusivity_fields, model.clock, fields(model), model.buoyancy) + DIFF = TracerEquation.TracerDiffusion(model, :a, model.tracers.a, model.closure, model.closure_fields, model.clock, fields(model), model.buoyancy) DIFF_field = Field(DIFF) @test DIFF isa TracerEquation.Diffusion @test DIFF isa TracerDiffusion @@ -75,7 +75,7 @@ function test_tracer_terms(model) @test DIFF_field isa Field DIFF = TracerEquation.ImmersedDiffusion(model, model.tracers.a, model.tracers.a.boundary_conditions.immersed, - model.closure, model.diffusivity_fields, Val(:a), model.clock, fields(model)) + model.closure, model.closure_fields, Val(:a), model.clock, fields(model)) DIFF_field = Field(DIFF) @test DIFF isa TracerEquation.ImmersedDiffusion @test DIFF isa TracerImmersedDiffusion @@ -88,7 +88,7 @@ function test_tracer_terms(model) @test DIFF_field isa Field DIFF = TracerEquation.TotalDiffusion(model, model.tracers.a, model.tracers.a.boundary_conditions.immersed, - model.closure, model.diffusivity_fields, Val(:a), model.clock, fields(model), model.buoyancy) + model.closure, model.closure_fields, Val(:a), model.clock, fields(model), model.buoyancy) DIFF_field = Field(DIFF) @test DIFF isa TracerEquation.TotalDiffusion @test DIFF isa TracerTotalDiffusion diff --git a/test/test_turbulent_kinetic_energy_equation.jl b/test/test_turbulent_kinetic_energy_equation.jl index f077d406..c7add443 100644 --- a/test/test_turbulent_kinetic_energy_equation.jl +++ b/test/test_turbulent_kinetic_energy_equation.jl @@ -109,9 +109,9 @@ function test_ke_dissipation_rate_terms(grid; model_type=NonhydrostaticModel, cl idxs = (model.grid.Nx÷2, model.grid.Ny÷2, model.grid.Nz÷2) if closure isa Tuple - ν_field = Field(sum(viscosity(closure, model.diffusivity_fields))) + ν_field = Field(sum(viscosity(closure, model.closure_fields))) else - ν_field = viscosity(closure, model.diffusivity_fields) + ν_field = viscosity(closure, model.closure_fields) end rtol = zspacings(grid, Center()) isa Number ? 1e-12 : 0.06 # less accurate for stretched grid diff --git a/test/test_utils.jl b/test/test_utils.jl index 09f57ada..5bc7f635 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -1,6 +1,7 @@ using CUDA: has_cuda_gpu using Oceananigans using Oceananigans.TurbulenceClosures.Smagorinskys: LagrangianAveraging +using Oceananigans.TurbulenceClosures: DynamicCoefficient, LagrangianAveraging using SeawaterPolynomials: RoquetEquationOfState, TEOS10EquationOfState #+++ Common grid setup @@ -40,8 +41,8 @@ model_kwargs = (buoyancy = BuoyancyForce(BuoyancyTracer()), # Common closures closures = (ScalarDiffusivity(ν=1e-6, κ=1e-7), SmagorinskyLilly(), - Smagorinsky(coefficient=DynamicCoefficient(averaging=(1, 2))), - Smagorinsky(coefficient=DynamicCoefficient(averaging=LagrangianAveraging())), + DynamicSmagorinsky(averaging=(1, 2)), + DynamicSmagorinsky(averaging=LagrangianAveraging()), (ScalarDiffusivity(ν=1e-6, κ=1e-7), AnisotropicMinimumDissipation()),) # Common buoyancy formulations