initarray! julia 1.8

docs/src/model_setup/lagrangian_particles.md

@@ -28,9 +28,11 @@ z₀ = -0.5 * ones(n_particles);
 lagrangian_particles = LagrangianParticles(x=x₀, y=y₀, z=z₀)
 
 # output
-10 Lagrangian particles with
+10 LagrangianParticles with eltype Particle:
 ├── 3 properties: (:x, :y, :z)
-└── 0 tracked fields: ()
+├── particle-wall restitution coefficient: 1.0
+├── 0 tracked fields: ()
+└── dynamics: no_dynamics
 ```
 
 then pass it to a model constructor
@@ -46,7 +48,7 @@ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
 ├── closure: Nothing
 ├── buoyancy: Nothing
 ├── coriolis: Nothing
-└── particles: 10 Lagrangian particles with 3 properties: (:x, :y, :z)
+└── particles: 10 LagrangianParticles with eltype Particle and properties (:x, :y, :z)
 ```
 
 !!! warn "Lagrangian particles on GPUs"
@@ -87,9 +89,11 @@ particles = StructArray{LagrangianMicrobe}((x₀, y₀, z₀, species, dna));
 lagrangian_particles = LagrangianParticles(particles)
 
 # output
-3 Lagrangian particles with
+3 LagrangianParticles with eltype LagrangianMicrobe:
 ├── 5 properties: (:x, :y, :z, :species, :dna)
-└── 0 tracked fields: ()
+├── particle-wall restitution coefficient: 1.0
+├── 0 tracked fields: ()
+└── dynamics: no_dynamics
 ```
 
 !!! warn "Custom properties on GPUs"

src/BoundaryConditions/fill_halo_regions.jl

@@ -9,6 +9,7 @@ using KernelAbstractions.Extras.LoopInfo: @unroll
 #####
 
 fill_halo_regions!(::Nothing, args...) = nothing
+fill_halo_regions!(::NamedTuple{(), Tuple{}}, args...) = nothing
 
 """
     fill_halo_regions!(fields::Union{Tuple, NamedTuple}, arch, args...)

src/BuoyancyModels/linear_equation_of_state.jl

@@ -9,8 +9,8 @@ struct LinearEquationOfState{FT} <: AbstractEquationOfState
 end
 
 Base.summary(eos::LinearEquationOfState) =
-    string("LinearEquationOfState(thermal_expansion=", scalar_summary(eos.thermal_expansion),
-                               ", haline_contraction=", scalar_summary(eos.haline_contraction), ")")
+    string("LinearEquationOfState(thermal_expansion=", prettysummary(eos.thermal_expansion),
+                               ", haline_contraction=", prettysummary(eos.haline_contraction), ")")
 
 Base.show(io, eos::LinearEquationOfState) = print(io, summary(eos))
 

src/BuoyancyModels/seawater_buoyancy.jl

@@ -1,5 +1,5 @@
 using Oceananigans.BoundaryConditions: NoFluxBoundaryCondition
-using Oceananigans.Grids: scalar_summary
+using Oceananigans.Utils: prettysummary
 
 """
     SeawaterBuoyancy{FT, EOS, T, S} <: AbstractBuoyancyModel{EOS}
@@ -20,7 +20,7 @@ required_tracers(::SeawaterBuoyancy{FT, EOS, <:Nothing, <:Number}) where {FT, EO
 required_tracers(::SeawaterBuoyancy{FT, EOS, <:Number, <:Nothing}) where {FT, EOS} = (:S,) # active salinity only
 
 Base.nameof(::Type{SeawaterBuoyancy}) = "SeawaterBuoyancy"
-Base.summary(b::SeawaterBuoyancy) = string(nameof(typeof(b)), " with g=", scalar_summary(b.gravitational_acceleration),
+Base.summary(b::SeawaterBuoyancy) = string(nameof(typeof(b)), " with g=", prettysummary(b.gravitational_acceleration),
                                            " and ", summary(b.equation_of_state))
 
 function Base.show(io::IO, b::SeawaterBuoyancy{FT}) where FT

src/Coriolis/beta_plane.jl

@@ -52,8 +52,8 @@ end
 @inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::BetaPlane, U) where FT = zero(FT)
 
 function Base.summary(βplane::BetaPlane{FT}) where FT 
-    fstr = scalar_summary(βplane.f₀)
-    βstr = scalar_summary(βplane.β)
+    fstr = prettysummary(βplane.f₀)
+    βstr = prettysummary(βplane.β)
     return "BetaPlane{$FT}(f₀=$fstr, β=$βstr)"
 end
 

src/Coriolis/f_plane.jl

@@ -1,4 +1,4 @@
-using Oceananigans.Grids: scalar_summary
+using Oceananigans.Utils: prettysummary
 
 """
     FPlane{FT} <: AbstractRotation
@@ -44,7 +44,7 @@ end
 @inline z_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = zero(eltype(grid))
 
 function Base.summary(fplane::FPlane{FT}) where FT 
-    fstr = scalar_summary(fplane.f)
+    fstr = prettysummary(fplane.f)
     return "FPlane{$FT}(f=$fstr)"
 end
 

src/Fields/field.jl

@@ -508,10 +508,12 @@ const MinimumReduction = typeof(Base.minimum!)
 const AllReduction     = typeof(Base.all!)
 const AnyReduction     = typeof(Base.any!)
 
-initialize_reduced_field!(::SumReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), Base.add_sum, true, interior(c))
-initialize_reduced_field!(::ProdReduction, f, r::ReducedField, c) = Base.initarray!(interior(r), Base.mul_prod, true, interior(c))
-initialize_reduced_field!(::AllReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), &, true, interior(c))
-initialize_reduced_field!(::AnyReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), |, true, interior(c))
+check_version(VERSION) = VERSION.minor > 7
+
+initialize_reduced_field!(::SumReduction,  f, r::ReducedField, c) = check_version(VERSION) ? Base.initarray!(interior(r), f, Base.add_sum, true, interior(c))  : Base.initarray!(interior(r), Base.add_sum, true, interior(c))
+initialize_reduced_field!(::ProdReduction, f, r::ReducedField, c) = check_version(VERSION) ? Base.initarray!(interior(r), f, Base.mul_prod, true, interior(c)) : Base.initarray!(interior(r), Base.mul_prod, true, interior(c))
+initialize_reduced_field!(::AllReduction,  f, r::ReducedField, c) = check_version(VERSION) ? Base.initarray!(interior(r), f, &, true, interior(c))             : Base.initarray!(interior(r), &, true, interior(c))
+initialize_reduced_field!(::AnyReduction,  f, r::ReducedField, c) = check_version(VERSION) ? Base.initarray!(interior(r), f, |, true, interior(c))             : Base.initarray!(interior(r), |, true, interior(c))
 
 initialize_reduced_field!(::MaximumReduction, f, r::ReducedField, c) = Base.mapfirst!(f, interior(r), interior(c))
 initialize_reduced_field!(::MinimumReduction, f, r::ReducedField, c) = Base.mapfirst!(f, interior(r), interior(c))

src/Fields/interpolate.jl

@@ -1,19 +1,39 @@
 using Oceananigans.Grids: XRegRectilinearGrid, YRegRectilinearGrid, ZRegRectilinearGrid
 using CUDA: allowscalar 
 
-@inline function linear_interpolate_sorted_vector(vec, val)
+function binary_search(list, query; rev=false, lt=<, by=identity)
+    if issorted(list) || issorted(list; rev=true)
+        low = !rev ? 0 : length(list) - 1
+        high = !rev ? length(list) - 1 : 0
+        middle(l, h) = round(Int, (l + h)//2)
+        query = by(query)
+
+        while !rev ? low + 1 < high : high + 1 < low
+            mid = middle(low, high)
+            by(list[mid+1]) === query && return mid:mid
+            if lt(by(list[mid+1]), query)
+                low = mid
+            else
+                high = mid
+            end
+        end
+        return !rev ? (low:high).+1 : (high:low).+1
 
-    allowscalar(true)
-    y2 = searchsortedfirst(vec, val)
-    y1 = searchsortedlast(vec, val)
-    x2 = vec[y2]
-    x1 = vec[y1]
-    allowscalar(false)
+    else
+        throw(error("List not sorted, unable to search value"))
+    end
+end
+
+@inline function fractional_index(vec, val)
+
+    y₁, y₂ = binary_search(vec, val)
+    x₂ = vec[y₂]
+    x₁ = vec[y₁]
 
-    if y1 == y2
-        return y1
+    if y₁ == y₂
+        return y₁
     else
-        return (y2 - y1) / (x2 - x1) * (val - x1) + y1
+        return (y₂ - y₁) / (x₂ - x₁) * (val - x₁) + y₁
     end
 end
 
@@ -22,20 +42,20 @@ end
 #### Use other methods if a more accurate interpolation is required
 ####
 
-@inline fractional_x_index(x, ::Face,   grid::RectilinearGrid)     = linear_interpolate_sorted_vector(grid.xᶠᵃᵃ, x)
-@inline fractional_x_index(x, ::Center, grid::RectilinearGrid)     = linear_interpolate_sorted_vector(grid.xᶜᵃᵃ, x)
-@inline fractional_y_index(y, ::Face,   grid::RectilinearGrid)     = linear_interpolate_sorted_vector(grid.yᵃᶠᵃ, y)
-@inline fractional_y_index(y, ::Center, grid::RectilinearGrid)     = linear_interpolate_sorted_vector(grid.yᵃᶜᵃ, y)
+@inline fractional_x_index(x, ::Face,   grid::RectilinearGrid) = fractional_index(grid.xᶠᵃᵃ, x)
+@inline fractional_x_index(x, ::Center, grid::RectilinearGrid) = fractional_index(grid.xᶜᵃᵃ, x)
+@inline fractional_y_index(y, ::Face,   grid::RectilinearGrid) = fractional_index(grid.yᵃᶠᵃ, y)
+@inline fractional_y_index(y, ::Center, grid::RectilinearGrid) = fractional_index(grid.yᵃᶜᵃ, y)
 
 @inline fractional_x_index(x, ::Face,   grid::XRegRectilinearGrid) = @inbounds (x - grid.xᶠᵃᵃ[1]) / grid.Δxᶠᵃᵃ
 @inline fractional_x_index(x, ::Center, grid::XRegRectilinearGrid) = @inbounds (x - grid.xᶜᵃᵃ[1]) / grid.Δxᶜᵃᵃ
 @inline fractional_y_index(y, ::Face,   grid::YRegRectilinearGrid) = @inbounds (y - grid.yᵃᶠᵃ[1]) / grid.Δyᵃᶠᵃ
 @inline fractional_y_index(y, ::Center, grid::YRegRectilinearGrid) = @inbounds (y - grid.yᵃᶜᵃ[1]) / grid.Δyᵃᶜᵃ
 
-@inline fractional_x_index(λ, ::Face,   grid::LatitudeLongitudeGrid) = linear_interpolate_sorted_vector(grid.λᶠᵃᵃ, λ)
-@inline fractional_x_index(λ, ::Center, grid::LatitudeLongitudeGrid) = linear_interpolate_sorted_vector(grid.λᶜᵃᵃ, λ)
-@inline fractional_y_index(φ, ::Face,   grid::LatitudeLongitudeGrid) = linear_interpolate_sorted_vector(grid.φᵃᶠᵃ, φ)
-@inline fractional_y_index(φ, ::Center, grid::LatitudeLongitudeGrid) = linear_interpolate_sorted_vector(grid.φᵃᶜᵃ, φ)
+@inline fractional_x_index(λ, ::Face,   grid::LatitudeLongitudeGrid) = fractional_index(grid.λᶠᵃᵃ, λ)
+@inline fractional_x_index(λ, ::Center, grid::LatitudeLongitudeGrid) = fractional_index(grid.λᶜᵃᵃ, λ)
+@inline fractional_y_index(φ, ::Face,   grid::LatitudeLongitudeGrid) = fractional_index(grid.φᵃᶠᵃ, φ)
+@inline fractional_y_index(φ, ::Center, grid::LatitudeLongitudeGrid) = fractional_index(grid.φᵃᶜᵃ, φ)
 
 @inline fractional_x_index(λ, ::Face,   grid::XRegLatLonGrid) = @inbounds (λ - grid.λᶠᵃᵃ[1]) / grid.Δλᶠᵃᵃ
 @inline fractional_x_index(λ, ::Center, grid::XRegLatLonGrid) = @inbounds (λ - grid.λᶜᵃᵃ[1]) / grid.Δλᶜᵃᵃ
@@ -44,8 +64,8 @@ end
 
 const ZReg = Union{ZRegRectilinearGrid, ZRegLatLonGrid}
 
-@inline fractional_z_index(z, ::Face,   grid) = linear_interpolate_sorted_vector(grid.zᵃᵃᶠ, z)
-@inline fractional_z_index(z, ::Center, grid) = linear_interpolate_sorted_vector(grid.zᵃᵃᶜ, z)
+@inline fractional_z_index(z, ::Face,   grid) = fractional_index(grid.zᵃᵃᶠ, z)
+@inline fractional_z_index(z, ::Center, grid) = fractional_index(grid.zᵃᵃᶜ, z)
 
 @inline fractional_z_index(z, ::Face,   grid::ZReg) = @inbounds (z - grid.zᵃᵃᶠ[1]) / grid.Δzᵃᵃᶠ
 @inline fractional_z_index(z, ::Center, grid::ZReg) = @inbounds (z - grid.zᵃᵃᶜ[1]) / grid.Δzᵃᵃᶜ

src/Fields/show_fields.jl

@@ -1,5 +1,5 @@
 using Printf
-using Oceananigans.Grids: size_summary, scalar_summary
+using Oceananigans.Grids: size_summary
 using Oceananigans.Utils: prettysummary
 using Oceananigans.BoundaryConditions: bc_str
 

src/Grids/grid_utils.jl

@@ -403,7 +403,7 @@ Base.summary(::ZDirection) = "ZDirection"
 #####
 
 size_summary(sz) = string(sz[1], "×", sz[2], "×", sz[3])
-scalar_summary(σ) = writeshortest(σ, false, false, true, -1, UInt8('e'), false, UInt8('.'), false, true)
+prettysummary(σ::AbstractFloat, plus=false) = writeshortest(σ, plus, false, true, -1, UInt8('e'), false, UInt8('.'), false, true)
 dimension_summary(topo::Flat, name, args...) = "Flat $name"
 
 function domain_summary(topo, name, left, right)
@@ -414,10 +414,10 @@ function domain_summary(topo, name, left, right)
                   topo isa FullyConnected ? "FullyConnected " :
                   topo isa LeftConnected ? "LeftConnected  " :
                   "RightConnected "
-    
-    prefix = string(topo_string, name, " ∈ [",
-                    scalar_summary(left), ", ",
-                    scalar_summary(right), interval)
+
+    return string(topo_string, name, " ∈ [",
+                  prettysummary(left), ", ",
+                  prettysummary(right), interval)
 end
 
 function dimension_summary(topo, name, left, right, spacing, pad_domain=0)
@@ -426,7 +426,7 @@ function dimension_summary(topo, name, left, right, spacing, pad_domain=0)
     return string(prefix, padding, coordinate_summary(spacing, name))
 end
 
-coordinate_summary(Δ::Number, name) = @sprintf("regularly spaced with Δ%s=%s", name, scalar_summary(Δ))
+coordinate_summary(Δ::Number, name) = @sprintf("regularly spaced with Δ%s=%s", name, prettysummary(Δ))
 coordinate_summary(Δ::AbstractVector, name) = @sprintf("variably spaced with min(Δ%s)=%s, max(Δ%s)=%s",
-                                                       name, scalar_summary(minimum(parent(Δ))),
-                                                       name, scalar_summary(maximum(parent(Δ))))
+                                                       name, prettysummary(minimum(parent(Δ))),
+                                                       name, prettysummary(maximum(parent(Δ))))

src/LagrangianParticleTracking/LagrangianParticleTracking.jl

@@ -2,13 +2,19 @@ module LagrangianParticleTracking
 
 export LagrangianParticles, update_particle_properties!
 
+using Printf
 using Adapt
 using KernelAbstractions
 using StructArrays
 
 using Oceananigans.Grids
-using Oceananigans.Architectures: device
-using Oceananigans.Fields: interpolate, datatuple, compute!, location
+using Oceananigans.Grids: xnode, ynode, znode
+using Oceananigans.Grids: AbstractUnderlyingGrid, AbstractGrid, hack_cosd
+using Oceananigans.ImmersedBoundaries
+using Oceananigans.ImmersedBoundaries: immersed_cell
+using Oceananigans.Architectures: device, architecture
+using Oceananigans.Fields: interpolate, datatuple, compute!, location, fractional_indices, fractional_y_index
+using Oceananigans.Utils: prettysummary, launch!
 
 import Base: size, length, show
 
@@ -20,6 +26,11 @@ struct Particle{T} <: AbstractParticle
     z :: T
 end
 
+Base.show(io::IO, p::Particle) = print(io, "Particle at (",
+                                       @sprintf("%-8s", prettysummary(p.x, true) * ", "),
+                                       @sprintf("%-8s", prettysummary(p.y, true) * ", "),
+                                       @sprintf("%-8s", prettysummary(p.z, true) * ")"))
+
 struct LagrangianParticles{P, R, T, D, Π}
         properties :: P
        restitution :: R
@@ -79,13 +90,22 @@ end
 size(lagrangian_particles::LagrangianParticles) = size(lagrangian_particles.properties)
 length(lagrangian_particles::LagrangianParticles) = length(lagrangian_particles.properties)
 
+Base.summary(particles::LagrangianParticles) =
+    string(length(particles), " LagrangianParticles with eltype ", nameof(eltype(particles.properties)),
+           " and properties ", propertynames(particles.properties))
+
 function Base.show(io::IO, lagrangian_particles::LagrangianParticles)
     particles = lagrangian_particles.properties
+    Tparticle = nameof(eltype(particles))
     properties = propertynames(particles)
     fields = lagrangian_particles.tracked_fields
-    print(io, "$(length(particles)) Lagrangian particles with\n",
-        "├── $(length(properties)) properties: $properties\n",
-        "└── $(length(fields)) tracked fields: $(propertynames(fields))")
+    Nparticles = length(particles)
+
+    print(io, Nparticles, " LagrangianParticles with eltype ", Tparticle, ":", '\n',
+        "├── ", length(properties), " properties: ", properties, '\n',
+        "├── particle-wall restitution coefficient: ", lagrangian_particles.restitution, '\n',
+        "├── ", length(fields), " tracked fields: ", propertynames(fields), '\n',
+        "└── dynamics: ", prettysummary(lagrangian_particles.dynamics, false))
 end
 
 include("update_particle_properties.jl")