From e871835fa0a5d8cf6e1516a5b40f271184d15b46 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 30 Mar 2022 18:54:19 -0400 Subject: [PATCH 01/30] fixed interpolation --- src/Fields/interpolate.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Fields/interpolate.jl b/src/Fields/interpolate.jl index 723beb2ce8..0365696480 100644 --- a/src/Fields/interpolate.jl +++ b/src/Fields/interpolate.jl @@ -1,16 +1,16 @@ using Oceananigans.Grids: XRegRectilinearGrid, YRegRectilinearGrid, ZRegRectilinearGrid +using CUDA: allowscalar @inline function linear_interpolate_sorted_vector(vec, val) - @CUDA.allowscalar begin - y₂ = searchsortedfirst(vec, val) - y₁ = searchsortedlast(vec, val) + allowscalar(true) + y2 = searchsortedfirst(vec, val) + y1 = searchsortedlast(vec, val) + x2 = vec[y2] + x1 = vec[y1] + allowscalar(false) - x₂ = vec[y₂] - x₁ = vec[y₁] - end - - return (y₂ - y₁) / (x₂ - x₁) * (x - x₁) + y₁ + return (y2 - y1) / (x2 - x1) * (val - x1) + y1 end #### From 8d0924e7624f221b6d6895b3bada495cc2f3b8b0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Wed, 30 Mar 2022 19:15:31 -0400 Subject: [PATCH 02/30] Update src/Fields/interpolate.jl Co-authored-by: Gregory L. Wagner --- src/Fields/interpolate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Fields/interpolate.jl b/src/Fields/interpolate.jl index 0365696480..c425e13872 100644 --- a/src/Fields/interpolate.jl +++ b/src/Fields/interpolate.jl @@ -5,7 +5,7 @@ using CUDA: allowscalar allowscalar(true) y2 = searchsortedfirst(vec, val) - y1 = searchsortedlast(vec, val) + y1 = searchsortedlast(vec, val) x2 = vec[y2] x1 = vec[y1] allowscalar(false) From bb27a781bd7288a0ad624cf8200e3a61cc7c6e0b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:17:06 -0600 Subject: [PATCH 03/30] Test Lagrangian particles on stretched grid --- test/test_lagrangian_particle_tracking.jl | 103 +++++++++++++++------- 1 file changed, 69 insertions(+), 34 deletions(-) diff --git a/test/test_lagrangian_particle_tracking.jl b/test/test_lagrangian_particle_tracking.jl index 196c9750a1..08c3a995e3 100644 --- a/test/test_lagrangian_particle_tracking.jl +++ b/test/test_lagrangian_particle_tracking.jl @@ -2,6 +2,7 @@ include("dependencies_for_runtests.jl") using NCDatasets using StructArrays +using Oceananigans.Architectures: arch_array struct TestParticle{T} x :: T @@ -13,56 +14,87 @@ struct TestParticle{T} s :: T end -function run_simple_particle_tracking_tests(arch, timestepper) +function particle_tracking_simulation(; grid, particles, timestepper=:RungeKutta3, velocities=nothing) + model = NonhydrostaticModel(; grid, timestepper, velocities, particles) + set!(model, u=1, v=1) + sim = Simulation(model, Δt=1e-2, stop_iteration=1) + + jld2_filepath = "test_particles.jld2" + sim.output_writers[:particles_jld2] = + JLD2OutputWriter(model, (; particles=model.particles), + prefix="test_particles", schedule=IterationInterval(1), force=true) + + nc_filepath = "test_particles.nc" + sim.output_writers[:particles_nc] = + NetCDFOutputWriter(model, model.particles, filepath=nc_filepath, schedule=IterationInterval(1), mode="c") + + sim.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(1), + dir = ".", prefix = "particles_checkpoint") + + return sim +end + +function run_simple_particle_tracking_tests(arch, timestepper; vertically_stretched=false) topo = (Periodic, Periodic, Bounded) - domain = (x=(-1, 1), y=(-1, 1), z=(-1, 1)) - grid = RectilinearGrid(arch, topology=topo, size=(5, 5, 5); domain...) + + Nx = Ny = Nz = 5 + if vertically_stretched + # Slightly stretched at the top + z = [-1, -0.5, 0.0, 0.4, 0.7, 1] + else + z = (-1, 1) + end + + grid = RectilinearGrid(arch; topology=topo, size=(Nx, Ny, Nz), + x=(-1, 1), y=(-1, 1), z) P = 10 + + ##### + ##### Test default particle + ##### + xs = arch_array(arch, zeros(P)) ys = arch_array(arch, zeros(P)) zs = arch_array(arch, 0.5*ones(P)) - # Test first constructor - lagrangian_particles = LagrangianParticles(x=xs, y=ys, z=zs) - @test lagrangian_particles isa LagrangianParticles + particles = LagrangianParticles(x=xs, y=ys, z=zs) + @test particles isa LagrangianParticles + + sim = particle_tracking_simulation(; grid, particles, timestepper) + model = sim.model + run!(sim) + + # Just test we run without errors + @test length(model.particles) == P + @test propertynames(model.particles.properties) == (:x, :y, :z) - us = convert(array_type(arch), zeros(P)) - vs = convert(array_type(arch), zeros(P)) - ws = convert(array_type(arch), zeros(P)) - ss = convert(array_type(arch), zeros(P)) + ##### + ##### Test custom particle "SpeedTrackingParticle" + ##### + + xs = arch_array(arch, zeros(P)) + ys = arch_array(arch, zeros(P)) + zs = arch_array(arch, 0.5 * ones(P)) + us = arch_array(arch, zeros(P)) + vs = arch_array(arch, zeros(P)) + ws = arch_array(arch, zeros(P)) + ss = arch_array(arch, zeros(P)) + # Test custom constructor particles = StructArray{TestParticle}((xs, ys, zs, us, vs, ws, ss)) velocities = VelocityFields(grid) u, v, w = velocities speed = Field(√(u*u + v*v + w*w)) - tracked_fields = merge(velocities, (; s=speed)) # Test second constructor - lagrangian_particles = LagrangianParticles(particles; tracked_fields) - @test lagrangian_particles isa LagrangianParticles - - model = NonhydrostaticModel(grid=grid, timestepper=timestepper, - velocities=velocities, particles=lagrangian_particles) - - set!(model, u=1, v=1) - - sim = Simulation(model, Δt=1e-2, stop_iteration=1) - - jld2_filepath = "test_particles.jld2" - sim.output_writers[:particles_jld2] = - JLD2OutputWriter(model, (; particles=model.particles), - prefix="test_particles", schedule=IterationInterval(1)) - - nc_filepath = "test_particles.nc" - sim.output_writers[:particles_nc] = - NetCDFOutputWriter(model, model.particles, filepath=nc_filepath, schedule=IterationInterval(1)) - - sim.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(1), - dir = ".", prefix = "particles_checkpoint") + particles = LagrangianParticles(particles; tracked_fields) + @test particles isa LagrangianParticles + sim = particle_tracking_simulation(; grid, particles, timestepper, velocities) + model = sim.model run!(sim) @test length(model.particles) == P @@ -184,7 +216,10 @@ end @testset "Lagrangian particle tracking" begin for arch in archs, timestepper in (:QuasiAdamsBashforth2, :RungeKutta3) - @info " Testing Lagrangian particle tacking [$(typeof(arch)), $timestepper]..." - run_simple_particle_tracking_tests(arch, timestepper) + @info " Testing uniform grid Lagrangian particle tacking [$(typeof(arch)), $timestepper]..." + run_simple_particle_tracking_tests(arch, timestepper; vertically_stretched=false) + + @info " Testing stretched grid Lagrangian particle tacking [$(typeof(arch)), $timestepper]..." + run_simple_particle_tracking_tests(arch, timestepper; vertically_stretched=true) end end From c909cb008be03c8e9c861cfa4ad787ccd586d72f Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:17:19 -0600 Subject: [PATCH 04/30] Adds a validation experiment for Lagrangian particles --- .../particles_in_convection.jl | 56 +++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 validation/lagrangian_particles/particles_in_convection.jl diff --git a/validation/lagrangian_particles/particles_in_convection.jl b/validation/lagrangian_particles/particles_in_convection.jl new file mode 100644 index 0000000000..6cc1c2d892 --- /dev/null +++ b/validation/lagrangian_particles/particles_in_convection.jl @@ -0,0 +1,56 @@ +using Random +using Printf +using Oceananigans +using Oceananigans.Units: minute, minutes, hour + +# Stretched grid +Lz = 32 # (m) domain depth +Lx = Ly = 64 # domain width +Nx = Ny = 32 +Nz = 24 # number of points in the vertical direction + +refinement = 1.2 # controls spacing near surface (higher means finer spaced) +stretching = 12 # controls rate of stretching at bottom +h(k) = (k - 1) / Nz +ζ₀(k) = 1 + (h(k) - 1) / refinement +Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching)) + +# Vertically-stretched +z(k) = Lz * (ζ₀(k) * Σ(k) - 1) + +# Vertically-uniform +#z = (-24, 0) + +grid = RectilinearGrid(size = (Nx, Ny, Nz), halo=(3, 3, 3), + x = (-Lx/2, Lx/2), + y = (-Ly/2, Lx/2), + z = (-24, 0)) + +# 10 Lagrangian particles +Nparticles = 10 +x₀ = Lx / 10 * rand(Nparticles) +y₀ = Ly / 10 * rand(Nparticles) +z₀ = - Lz / 10 * rand(Nparticles) .+ 3Lz / 5 +particles = LagrangianParticles(x=x₀, y=y₀, z=z₀, restitution=0) + +# Convection +b_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(1e-8)) + +model = NonhydrostaticModel(; grid, particles, + advection = UpwindBiasedFifthOrder(), + timestepper = :RungeKutta3, + tracers = :b, + buoyancy = BuoyancyTracer(), + closure = AnisotropicMinimumDissipation(), + boundary_conditions = (; b=b_bcs)) + +bᵢ(x, y, z) = 1e-5 * z + 1e-9 * rand() +set!(model, b=bᵢ) + +simulation = Simulation(model, Δt=10.0, stop_time=20minutes) +wizard = TimeStepWizard(cfl=0.5, max_change=1.1, max_Δt=1minute) + +simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + +run!(simulation) + From 5c35f5dcf1ca109938a66913b05cdbc2f1616509 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:22:42 -0600 Subject: [PATCH 05/30] Couple bugfixes in new test --- test/test_lagrangian_particle_tracking.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_lagrangian_particle_tracking.jl b/test/test_lagrangian_particle_tracking.jl index 08c3a995e3..75ced8b2e1 100644 --- a/test/test_lagrangian_particle_tracking.jl +++ b/test/test_lagrangian_particle_tracking.jl @@ -31,7 +31,7 @@ function particle_tracking_simulation(; grid, particles, timestepper=:RungeKutta sim.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(1), dir = ".", prefix = "particles_checkpoint") - return sim + return sim, jld2_filepath, nc_filepath end function run_simple_particle_tracking_tests(arch, timestepper; vertically_stretched=false) @@ -61,7 +61,7 @@ function run_simple_particle_tracking_tests(arch, timestepper; vertically_stretc particles = LagrangianParticles(x=xs, y=ys, z=zs) @test particles isa LagrangianParticles - sim = particle_tracking_simulation(; grid, particles, timestepper) + sim, _, _ = particle_tracking_simulation(; grid, particles, timestepper) model = sim.model run!(sim) @@ -93,7 +93,7 @@ function run_simple_particle_tracking_tests(arch, timestepper; vertically_stretc particles = LagrangianParticles(particles; tracked_fields) @test particles isa LagrangianParticles - sim = particle_tracking_simulation(; grid, particles, timestepper, velocities) + sim, jld2_filepath, nc_filepath = particle_tracking_simulation(; grid, particles, timestepper, velocities) model = sim.model run!(sim) From 4917bf2e8c04a4b6e78011bfa2d0de63ad191ef6 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:31:52 -0600 Subject: [PATCH 06/30] Updates validation --- .../particles_in_convection.jl | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/validation/lagrangian_particles/particles_in_convection.jl b/validation/lagrangian_particles/particles_in_convection.jl index 6cc1c2d892..bba7659842 100644 --- a/validation/lagrangian_particles/particles_in_convection.jl +++ b/validation/lagrangian_particles/particles_in_convection.jl @@ -26,13 +26,19 @@ grid = RectilinearGrid(size = (Nx, Ny, Nz), halo=(3, 3, 3), y = (-Ly/2, Lx/2), z = (-24, 0)) +@info "Build a grid:" +@show grid + # 10 Lagrangian particles Nparticles = 10 -x₀ = Lx / 10 * rand(Nparticles) -y₀ = Ly / 10 * rand(Nparticles) -z₀ = - Lz / 10 * rand(Nparticles) .+ 3Lz / 5 +x₀ = Lx / 10 * (2rand(Nparticles) .- 1) +y₀ = Ly / 10 * (2rand(Nparticles) .- 1) +z₀ = - Lz / 10 * rand(Nparticles) particles = LagrangianParticles(x=x₀, y=y₀, z=z₀, restitution=0) +@info "Initialized Lagrangian particles" +@show particles + # Convection b_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(1e-8)) @@ -44,6 +50,9 @@ model = NonhydrostaticModel(; grid, particles, closure = AnisotropicMinimumDissipation(), boundary_conditions = (; b=b_bcs)) +@info "Constructed a model" +@show model + bᵢ(x, y, z) = 1e-5 * z + 1e-9 * rand() set!(model, b=bᵢ) From 1f4cf69d1970225a9f5a0e04f58fe96eabd159e8 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:44:12 -0600 Subject: [PATCH 07/30] Improves show method --- .../LagrangianParticleTracking.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index ae392edc43..adfd85fce5 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -9,6 +9,7 @@ using StructArrays using Oceananigans.Grids using Oceananigans.Architectures: device using Oceananigans.Fields: interpolate, datatuple, compute!, location +using Oceananigans.Utils: prettysummary import Base: size, length, show @@ -81,11 +82,16 @@ length(lagrangian_particles::LagrangianParticles) = length(lagrangian_particles. 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") From 8452ade4db46d2ca65f49ca188437b4ddcacc018 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:50:21 -0600 Subject: [PATCH 08/30] Updates validation experiment --- .../lagrangian_particles/particles_in_convection.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/validation/lagrangian_particles/particles_in_convection.jl b/validation/lagrangian_particles/particles_in_convection.jl index bba7659842..be5404d40a 100644 --- a/validation/lagrangian_particles/particles_in_convection.jl +++ b/validation/lagrangian_particles/particles_in_convection.jl @@ -15,16 +15,14 @@ h(k) = (k - 1) / Nz ζ₀(k) = 1 + (h(k) - 1) / refinement Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching)) -# Vertically-stretched -z(k) = Lz * (ζ₀(k) * Σ(k) - 1) +# Vertically-stretched and uniform options +z_stretched(k) = Lz * (ζ₀(k) * Σ(k) - 1) +z_uniform = (-Lz, 0) -# Vertically-uniform -#z = (-24, 0) - -grid = RectilinearGrid(size = (Nx, Ny, Nz), halo=(3, 3, 3), +grid = RectilinearGrid(; size = (Nx, Ny, Nz), halo=(3, 3, 3), x = (-Lx/2, Lx/2), y = (-Ly/2, Lx/2), - z = (-24, 0)) + z = z_stretched) @info "Build a grid:" @show grid From f35fa9bba6894eff010ead7183988c6d8fab7be2 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 06:50:30 -0600 Subject: [PATCH 09/30] Better show again --- src/LagrangianParticleTracking/LagrangianParticleTracking.jl | 4 ++++ src/Models/NonhydrostaticModels/show_nonhydrostatic_model.jl | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index adfd85fce5..8c67634e3f 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -80,6 +80,10 @@ 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)) diff --git a/src/Models/NonhydrostaticModels/show_nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/show_nonhydrostatic_model.jl index 850f11328f..046c0810cc 100644 --- a/src/Models/NonhydrostaticModels/show_nonhydrostatic_model.jl +++ b/src/Models/NonhydrostaticModels/show_nonhydrostatic_model.jl @@ -25,7 +25,7 @@ function Base.show(io::IO, model::NonhydrostaticModel) particles = model.particles.properties properties = propertynames(particles) print(io, "├── coriolis: ", summary(model.coriolis), '\n') - print(io, "└── particles: $(length(particles)) Lagrangian particles with $(length(properties)) properties: $properties") + print(io, "└── particles: ", summary(model.particles)) end end From fcc6fffc98bdf141935bac6ddef69102d3d3a88e Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 07:02:17 -0600 Subject: [PATCH 10/30] Updates show for Particle so we can see their location --- src/Fields/interpolate.jl | 22 +++++++++---------- src/Grids/grid_utils.jl | 12 +++++----- .../LagrangianParticleTracking.jl | 6 +++++ src/Utils/prettysummary.jl | 3 +-- .../particles_in_convection.jl | 4 ++-- 5 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/Fields/interpolate.jl b/src/Fields/interpolate.jl index c425e13872..d2e8196d89 100644 --- a/src/Fields/interpolate.jl +++ b/src/Fields/interpolate.jl @@ -1,7 +1,7 @@ using Oceananigans.Grids: XRegRectilinearGrid, YRegRectilinearGrid, ZRegRectilinearGrid using CUDA: allowscalar -@inline function linear_interpolate_sorted_vector(vec, val) +@inline function fractional_index(vec, val) allowscalar(true) y2 = searchsortedfirst(vec, val) @@ -18,20 +18,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.Δλᶜᵃᵃ @@ -40,8 +40,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ᵃᵃᶜ diff --git a/src/Grids/grid_utils.jl b/src/Grids/grid_utils.jl index f62f282e11..2e86655eca 100644 --- a/src/Grids/grid_utils.jl +++ b/src/Grids/grid_utils.jl @@ -395,7 +395,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) @@ -404,8 +404,8 @@ function domain_summary(topo, name, left, right) "Bounded " return string(topo_string, name, " ∈ [", - scalar_summary(left), ", ", - scalar_summary(right), interval) + prettysummary(left), ", ", + prettysummary(right), interval) end function dimension_summary(topo, name, left, right, spacing, pad_domain=0) @@ -414,7 +414,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(Δ)))) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index 8c67634e3f..b41e7275fe 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -2,6 +2,7 @@ module LagrangianParticleTracking export LagrangianParticles, update_particle_properties! +using Printf using Adapt using KernelAbstractions using StructArrays @@ -21,6 +22,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 diff --git a/src/Utils/prettysummary.jl b/src/Utils/prettysummary.jl index 428d71a4d2..382ce8ed71 100644 --- a/src/Utils/prettysummary.jl +++ b/src/Utils/prettysummary.jl @@ -1,4 +1,4 @@ -using Oceananigans.Grids: scalar_summary +import Oceananigans.Grids: prettysummary prettysummary(x, args...) = summary(x) @@ -18,7 +18,6 @@ function prettysummary(f::Function, showmethods=true) end end -prettysummary(x::AbstractFloat, args...) = scalar_summary(x) prettysummary(x::Int, args...) = string(x) # This is very important diff --git a/validation/lagrangian_particles/particles_in_convection.jl b/validation/lagrangian_particles/particles_in_convection.jl index be5404d40a..c0fdb794be 100644 --- a/validation/lagrangian_particles/particles_in_convection.jl +++ b/validation/lagrangian_particles/particles_in_convection.jl @@ -54,10 +54,10 @@ model = NonhydrostaticModel(; grid, particles, bᵢ(x, y, z) = 1e-5 * z + 1e-9 * rand() set!(model, b=bᵢ) -simulation = Simulation(model, Δt=10.0, stop_time=20minutes) +simulation = Simulation(model, Δt=10.0, stop_iteration=1) wizard = TimeStepWizard(cfl=0.5, max_change=1.1, max_Δt=1minute) simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) -run!(simulation) +# run!(simulation) From 18407783d75333737ba04d589ca34a1c78b8de22 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 31 Mar 2022 07:08:28 -0600 Subject: [PATCH 11/30] scalar_summary -> prettysummary --- src/BuoyancyModels/linear_equation_of_state.jl | 4 ++-- src/BuoyancyModels/seawater_buoyancy.jl | 4 ++-- src/Coriolis/beta_plane.jl | 4 ++-- src/Coriolis/f_plane.jl | 4 ++-- src/Fields/show_fields.jl | 2 +- src/Simulations/time_step_wizard.jl | 8 ++++---- 6 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/BuoyancyModels/linear_equation_of_state.jl b/src/BuoyancyModels/linear_equation_of_state.jl index 1c059cb33b..1b570586fe 100644 --- a/src/BuoyancyModels/linear_equation_of_state.jl +++ b/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)) diff --git a/src/BuoyancyModels/seawater_buoyancy.jl b/src/BuoyancyModels/seawater_buoyancy.jl index 2fd48c1308..fc763a9896 100644 --- a/src/BuoyancyModels/seawater_buoyancy.jl +++ b/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 diff --git a/src/Coriolis/beta_plane.jl b/src/Coriolis/beta_plane.jl index f1b70db91a..dbf5949b17 100644 --- a/src/Coriolis/beta_plane.jl +++ b/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 diff --git a/src/Coriolis/f_plane.jl b/src/Coriolis/f_plane.jl index e0126d757a..1069a9eda5 100644 --- a/src/Coriolis/f_plane.jl +++ b/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 diff --git a/src/Fields/show_fields.jl b/src/Fields/show_fields.jl index 34fb2a9789..319a17fc11 100644 --- a/src/Fields/show_fields.jl +++ b/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 diff --git a/src/Simulations/time_step_wizard.jl b/src/Simulations/time_step_wizard.jl index 01271fe64e..970695da24 100644 --- a/src/Simulations/time_step_wizard.jl +++ b/src/Simulations/time_step_wizard.jl @@ -1,5 +1,5 @@ using Oceananigans: TurbulenceClosures -using Oceananigans.Grids: scalar_summary +using Oceananigans.Grids: prettysummary mutable struct TimeStepWizard{FT, C, D} cfl :: FT @@ -15,9 +15,9 @@ end infinite_diffusion_timescale(args...) = Inf # its not very limiting Base.summary(wizard::TimeStepWizard) = string("TimeStepWizard(", - "cfl=", scalar_summary(wizard.cfl), - ", max_Δt=", scalar_summary(wizard.max_Δt), - ", min_Δt=", scalar_summary(wizard.min_Δt), ")") + "cfl=", prettysummary(wizard.cfl), + ", max_Δt=", prettysummary(wizard.max_Δt), + ", min_Δt=", prettysummary(wizard.min_Δt), ")") """ TimeStepWizard(cfl=0.2, diffusive_cfl=Inf, max_change=1.1, min_change=0.5, max_Δt=Inf, min_Δt=0.0) From 5800c1e69f5e92ab84327b00d3bdb0484cb48925 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Wed, 13 Jul 2022 10:19:09 -0400 Subject: [PATCH 12/30] Update prescribed_hydrostatic_velocity_fields.jl --- .../prescribed_hydrostatic_velocity_fields.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl index c89cdf8990..4a575aebec 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl @@ -7,6 +7,7 @@ using Oceananigans.Fields: AbstractField, FunctionField import Oceananigans.BoundaryConditions: fill_halo_regions! import Oceananigans.Models.NonhydrostaticModels: extract_boundary_conditions +import Oceananigans.Utils: datatuple using Adapt @@ -69,6 +70,8 @@ end @inline fill_halo_regions!(::PrescribedVelocityFields, args...) = nothing @inline fill_halo_regions!(::FunctionField, args...) = nothing +@inline datatuple(obj::PrescribedVelocityFields) = (; u = datatuple(obj.u), v = datatuple(obj.v), w = datatuple(obj.w)) + ab2_step_velocities!(::PrescribedVelocityFields, args...) = [NoneEvent()] ab2_step_free_surface!(::Nothing, args...) = NoneEvent() compute_w_from_continuity!(::PrescribedVelocityFields, args...) = nothing From 93f0d28c562319919b9a7f4c58ed4980cdf1e24b Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 10:29:04 -0400 Subject: [PATCH 13/30] typo --- src/TimeSteppers/quasi_adams_bashforth_2.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TimeSteppers/quasi_adams_bashforth_2.jl b/src/TimeSteppers/quasi_adams_bashforth_2.jl index cb8b00d9ed..c7e5542620 100644 --- a/src/TimeSteppers/quasi_adams_bashforth_2.jl +++ b/src/TimeSteppers/quasi_adams_bashforth_2.jl @@ -94,7 +94,7 @@ function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities. calculate_pressure_correction!(model, Δt) - @apply_regionally correct_velocties_and_store_tendecies!(model, Δt) + @apply_regionally correct_velocities_and_store_tendecies!(model, Δt) tick!(model.clock, Δt) update_state!(model) @@ -103,7 +103,7 @@ function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt return nothing end -function correct_velocties_and_store_tendecies!(model, Δt) +function correct_velocities_and_store_tendecies!(model, Δt) pressure_correct_velocities!(model, Δt) store_tendencies!(model) end From 23be5196b2a81864a6b16a2425aae8a9c5457561 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 10:51:01 -0400 Subject: [PATCH 14/30] add preliminary IBG --- .../LagrangianParticleTracking.jl | 5 ++-- .../update_particle_properties.jl | 28 ++++++++++++------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index b41e7275fe..7e73aba7b9 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -8,9 +8,10 @@ using KernelAbstractions using StructArrays using Oceananigans.Grids -using Oceananigans.Architectures: device +using Oceananigans.ImmersedBoundaries +using Oceananigans.Architectures: device, architecture using Oceananigans.Fields: interpolate, datatuple, compute!, location -using Oceananigans.Utils: prettysummary +using Oceananigans.Utils: prettysummary, launch! import Base: size, length, show diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 683a59d0df..dc749b53c3 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -1,5 +1,3 @@ -using Oceananigans.Grids: architecture - """ enforce_boundary_conditions(x, xₗ, xᵣ, ::Bounded) @@ -38,6 +36,20 @@ end @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz], restitution) end +@kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundayGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} + p = @index(Global) + + # Advect particles using forward Euler. + @inbounds particles.x[p] += interpolate(velocities.u, Face(), Center(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + + # Enforce boundary conditions for particles. + @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], grid.xᶠᵃᵃ[1], grid.xᶠᵃᵃ[grid.Nx], restitution) + @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], grid.yᵃᶠᵃ[1], grid.yᵃᶠᵃ[grid.Ny], restitution) + @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz], restitution) +end + @kernel function update_field_property!(particle_property, particles, grid, field, LX, LY, LZ) p = @index(Global) @@ -77,15 +89,11 @@ function update_particle_properties!(lagrangian_particles, model, Δt) - # advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) - - advect_particles_event = launch!(arch, model.grid, worksize, _advect_particles!, - lagrangian_particles.properties, lagrangian_particles.restitution, model.grid, Δt, - datatuple(model.velocities); dependencies=Event(device(arch))) + advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) - # advect_particles_event = advect_particles_kernel!(lagrangian_particles.properties, lagrangian_particles.restitution, model.grid, Δt, - # datatuple(model.velocities), - # dependencies=Event(device(arch))) + advect_particles_event = advect_particles_kernel!(lagrangian_particles.properties, lagrangian_particles.restitution, model.grid, Δt, + datatuple(model.velocities), + dependencies=Event(device(arch))) wait(device(arch), advect_particles_event) From c2900fd7bb584c8f3f957213b46abb36faf2103a Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 10:53:22 -0400 Subject: [PATCH 15/30] fixed bugs --- test/test_lagrangian_particle_tracking.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_lagrangian_particle_tracking.jl b/test/test_lagrangian_particle_tracking.jl index e1c023b26d..d52847108f 100644 --- a/test/test_lagrangian_particle_tracking.jl +++ b/test/test_lagrangian_particle_tracking.jl @@ -22,7 +22,7 @@ function particle_tracking_simulation(; grid, particles, timestepper=:RungeKutta jld2_filepath = "test_particles.jld2" sim.output_writers[:particles_jld2] = JLD2OutputWriter(model, (; particles=model.particles), - prefix="test_particles", schedule=IterationInterval(1), force=true) + filename="test_particles", schedule=IterationInterval(1), overwrite_existing=true) nc_filepath = "test_particles.nc" sim.output_writers[:particles_nc] = From 74d62c92f65f012ac6d81d4c50b9dc1ef576d91e Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 10:56:51 -0400 Subject: [PATCH 16/30] IBG before particles --- .../update_particle_properties.jl | 8 ++++---- src/Oceananigans.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index dc749b53c3..2c7efbd667 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -36,13 +36,13 @@ end @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz], restitution) end -@kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundayGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} +@kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} p = @index(Global) # Advect particles using forward Euler. - @inbounds particles.x[p] += interpolate(velocities.u, Face(), Center(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + @inbounds particles.x[p] += interpolate(velocities.u, Face(), Center(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt # Enforce boundary conditions for particles. @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], grid.xᶠᵃᵃ[1], grid.xᶠᵃᵃ[grid.Nx], restitution) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index f8a88f65c3..ed290164c4 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -198,10 +198,10 @@ include("Coriolis/Coriolis.jl") include("BuoyancyModels/BuoyancyModels.jl") include("StokesDrift.jl") include("TurbulenceClosures/TurbulenceClosures.jl") -include("LagrangianParticleTracking/LagrangianParticleTracking.jl") include("Forcings/Forcings.jl") include("ImmersedBoundaries/ImmersedBoundaries.jl") +include("LagrangianParticleTracking/LagrangianParticleTracking.jl") include("TimeSteppers/TimeSteppers.jl") include("Models/Models.jl") From 66650c5184674d26ae864f89464427e9f6216a4f Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 11:35:03 -0400 Subject: [PATCH 17/30] added immersed restitution --- .../LagrangianParticleTracking.jl | 3 +- .../update_particle_properties.jl | 50 ++++++++++++++++--- 2 files changed, 46 insertions(+), 7 deletions(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index 7e73aba7b9..af9f10e303 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -8,9 +8,10 @@ using KernelAbstractions using StructArrays using Oceananigans.Grids +using Oceananigans.Grids: return_metrics using Oceananigans.ImmersedBoundaries using Oceananigans.Architectures: device, architecture -using Oceananigans.Fields: interpolate, datatuple, compute!, location +using Oceananigans.Fields: interpolate, datatuple, compute!, location, fractional_indices using Oceananigans.Utils: prettysummary, launch! import Base: size, length, show diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 2c7efbd667..c584b8d890 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -31,9 +31,9 @@ end @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt # Enforce boundary conditions for particles. - @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], grid.xᶠᵃᵃ[1], grid.xᶠᵃᵃ[grid.Nx], restitution) - @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], grid.yᵃᶠᵃ[1], grid.yᵃᶠᵃ[grid.Ny], restitution) - @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz], restitution) + @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], grid.xᶠᵃᵃ[1], grid.xᶠᵃᵃ[grid.Nx+1], restitution) + @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], grid.yᵃᶠᵃ[1], grid.yᵃᶠᵃ[grid.Ny+1], restitution) + @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz+1], restitution) end @kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} @@ -44,10 +44,48 @@ end @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt + x, y, z = return_face_metrics(grid) + # Enforce boundary conditions for particles. - @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], grid.xᶠᵃᵃ[1], grid.xᶠᵃᵃ[grid.Nx], restitution) - @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], grid.yᵃᶠᵃ[1], grid.yᵃᶠᵃ[grid.Ny], restitution) - @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz], restitution) + @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], x[1], x[grid.Nx], restitution) + @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], y[1], y[grid.Ny], restitution) + @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], z[1], z[grid.Nz], restitution) + + # Enforce Immersed boundary conditions for particles. + @inbounds enforce_immersed_boundary_condition!(particles, p, grid, restitution) +end + +@inline return_face_metrics(g::LatitudeLongitudeGrid) = (g.λᶠᵃᵃ, g.φᵃᶠᵃ, g.zᵃᵃᶠ) +@inline return_face_metrics(g::RectilinearGrid) = (g.xᶠᵃᵃ, g.yᵃᶠᵃ, g.zᵃᵃᶠ) +@inline return_face_metrics(g::ImmersedBoundaryGrid) = return_face_metrics(g.underlying_grid) + +@inline positive_x_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i-1, j, k, grid) +@inline positive_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j-1, k, grid) +@inline positive_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k-1, grid) + +@inline negative_x_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i+1, j, k, grid) +@inline negative_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j+1, k, grid) +@inline negative_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k+1, grid) + +@inline function enforce_immersed_boundary_condition_x(particles, p, grid, restitution) + xₚ, yₚ, zₚ = (particles.x[p], particles.y[p], particles.z[p]) + i, j, k = fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid) + + xᶠ⁻ = xnode(Face(), i, grid) + yᶠ⁻ = ynode(Face(), j, grid) + zᶠ⁻ = znode(Face(), k, grid) + + xᶠ⁺ = xnode(Face(), i + 1, grid) + yᶠ⁺ = ynode(Face(), j + 1, grid) + zᶠ⁺ = znode(Face(), k + 1, grid) + + positive_x_immersed_boundary(i, j, k, grid) && xₚ < xᶠ⁻ && xₚ = xᶠ⁻ + (xᶠ⁻ - xₚ) * restitution + positive_y_immersed_boundary(i, j, k, grid) && yₚ < yᶠ⁻ && yₚ = yᶠ⁻ + (yᶠ⁻ - yₚ) * restitution + positive_z_immersed_boundary(i, j, k, grid) && zₚ < zᶠ⁻ && zₚ = zᶠ⁻ + (zᶠ⁻ - zₚ) * restitution + + negative_x_immersed_boundary(i, j, k, grid) && xₚ > xᶠ⁺ && xₚ = xᶠ⁺ - (xₚ - xᶠ⁺) * restitution + negative_y_immersed_boundary(i, j, k, grid) && yₚ > yᶠ⁺ && yₚ = yᶠ⁺ - (yₚ - yᶠ⁺) * restitution + negative_z_immersed_boundary(i, j, k, grid) && zₚ > zᶠ⁺ && zₚ = zᶠ⁺ - (zₚ - zᶠ⁺) * restitution end @kernel function update_field_property!(particle_property, particles, grid, field, LX, LY, LZ) From 6b1ced758dda245613d77d85930f692e6878fed4 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 11:36:32 -0400 Subject: [PATCH 18/30] fixed test bug --- test/test_lagrangian_particle_tracking.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_lagrangian_particle_tracking.jl b/test/test_lagrangian_particle_tracking.jl index d52847108f..6de8bdba6e 100644 --- a/test/test_lagrangian_particle_tracking.jl +++ b/test/test_lagrangian_particle_tracking.jl @@ -26,7 +26,7 @@ function particle_tracking_simulation(; grid, particles, timestepper=:RungeKutta nc_filepath = "test_particles.nc" sim.output_writers[:particles_nc] = - NetCDFOutputWriter(model, model.particles, filepath=nc_filepath, schedule=IterationInterval(1), mode="c") + NetCDFOutputWriter(model, model.particles, filename=nc_filepath, schedule=IterationInterval(1), mode="c") sim.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(1), dir = ".", prefix = "particles_checkpoint") From 96b033c7097e5923ce7aaf3599e957b539a17042 Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 11:43:16 -0400 Subject: [PATCH 19/30] immersed restitution --- .../update_particle_properties.jl | 45 ++++++++++++------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index c584b8d890..3bea280ddc 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -52,7 +52,7 @@ end @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], z[1], z[grid.Nz], restitution) # Enforce Immersed boundary conditions for particles. - @inbounds enforce_immersed_boundary_condition!(particles, p, grid, restitution) + @inbounds particles.x[p], particles.y[p], particles.z[p] = enforce_immersed_boundary_condition(particles, p, grid, restitution) end @inline return_face_metrics(g::LatitudeLongitudeGrid) = (g.λᶠᵃᵃ, g.φᵃᶠᵃ, g.zᵃᵃᶠ) @@ -67,25 +67,36 @@ end @inline negative_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j+1, k, grid) @inline negative_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k+1, grid) -@inline function enforce_immersed_boundary_condition_x(particles, p, grid, restitution) +@inline function enforce_immersed_boundary_condition(particles, p, grid, restitution) xₚ, yₚ, zₚ = (particles.x[p], particles.y[p], particles.z[p]) - i, j, k = fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid) + i, j, k = Int.(floor(fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid))) - xᶠ⁻ = xnode(Face(), i, grid) - yᶠ⁻ = ynode(Face(), j, grid) - zᶠ⁻ = znode(Face(), k, grid) - - xᶠ⁺ = xnode(Face(), i + 1, grid) - yᶠ⁺ = ynode(Face(), j + 1, grid) - zᶠ⁺ = znode(Face(), k + 1, grid) - - positive_x_immersed_boundary(i, j, k, grid) && xₚ < xᶠ⁻ && xₚ = xᶠ⁻ + (xᶠ⁻ - xₚ) * restitution - positive_y_immersed_boundary(i, j, k, grid) && yₚ < yᶠ⁻ && yₚ = yᶠ⁻ + (yᶠ⁻ - yₚ) * restitution - positive_z_immersed_boundary(i, j, k, grid) && zₚ < zᶠ⁻ && zₚ = zᶠ⁻ + (zᶠ⁻ - zₚ) * restitution + if positive_x_immersed_boundary(i, j, k, grid) + xᶠ⁻ = xnode(Face(), i, grid) + xₚ < xᶠ⁻ && xₚ = xᶠ⁻ + (xᶠ⁻ - xₚ) * restitution + end + if positive_y_immersed_boundary(i, j, k, grid) + yᶠ⁻ = ynode(Face(), j, grid) + yₚ < yᶠ⁻ && yₚ = yᶠ⁻ + (yᶠ⁻ - yₚ) * restitution + end + if positive_z_immersed_boundary(i, j, k, grid) + zᶠ⁻ = znode(Face(), k, grid) + zₚ < zᶠ⁻ && zₚ = zᶠ⁻ + (zᶠ⁻ - zₚ) * restitution + end + if negative_x_immersed_boundary(i, j, k, grid) + xᶠ⁺ = xnode(Face(), i + 1, grid) + xₚ > xᶠ⁺ && xₚ = xᶠ⁺ - (xₚ - xᶠ⁺) * restitution + end + if negative_y_immersed_boundary(i, j, k, grid) + yᶠ⁺ = ynode(Face(), j + 1, grid) + yₚ > yᶠ⁺ && yₚ = yᶠ⁺ - (yₚ - yᶠ⁺) * restitution + end + if negative_z_immersed_boundary(i, j, k, grid) + zᶠ⁺ = znode(Face(), k + 1, grid) + zₚ > zᶠ⁺ && zₚ = zᶠ⁺ - (zₚ - zᶠ⁺) * restitution + end - negative_x_immersed_boundary(i, j, k, grid) && xₚ > xᶠ⁺ && xₚ = xᶠ⁺ - (xₚ - xᶠ⁺) * restitution - negative_y_immersed_boundary(i, j, k, grid) && yₚ > yᶠ⁺ && yₚ = yᶠ⁺ - (yₚ - yᶠ⁺) * restitution - negative_z_immersed_boundary(i, j, k, grid) && zₚ > zᶠ⁺ && zₚ = zᶠ⁺ - (zₚ - zᶠ⁺) * restitution + return xₚ, yₚ, zₚ end @kernel function update_field_property!(particle_property, particles, grid, field, LX, LY, LZ) From 5c18a0c8fdfdda133f25d19e4ac8cac8d226870a Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 12:06:18 -0400 Subject: [PATCH 20/30] angular velocity for particles --- .../LagrangianParticleTracking.jl | 2 +- .../update_particle_properties.jl | 59 +++++++++++-------- 2 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index af9f10e303..ccb421bc81 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -8,7 +8,7 @@ using KernelAbstractions using StructArrays using Oceananigans.Grids -using Oceananigans.Grids: return_metrics +using Oceananigans.Grids: xnode, ynode, znode using Oceananigans.ImmersedBoundaries using Oceananigans.Architectures: device, architecture using Oceananigans.Fields: interpolate, datatuple, compute!, location, fractional_indices diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 3bea280ddc..763416eba1 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -39,6 +39,8 @@ end @kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} p = @index(Global) + adapted_velocities = calc_correct_velocities(velocities, grid.underlying_grid) + # Advect particles using forward Euler. @inbounds particles.x[p] += interpolate(velocities.u, Face(), Center(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt @@ -55,6 +57,9 @@ end @inbounds particles.x[p], particles.y[p], particles.z[p] = enforce_immersed_boundary_condition(particles, p, grid, restitution) end +@inline calc_correct_velocities(vel, g::RectilinearGrid) = vel +@inline calc_correct_velocities(vel, g::LatitudeLongitudeGrid) = (; u = vel.u / g.radius, v = vel.v / g.radius, w = vel.w) + @inline return_face_metrics(g::LatitudeLongitudeGrid) = (g.λᶠᵃᵃ, g.φᵃᶠᵃ, g.zᵃᵃᶠ) @inline return_face_metrics(g::RectilinearGrid) = (g.xᶠᵃᵃ, g.yᵃᶠᵃ, g.zᵃᵃᶠ) @inline return_face_metrics(g::ImmersedBoundaryGrid) = return_face_metrics(g.underlying_grid) @@ -67,34 +72,40 @@ end @inline negative_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j+1, k, grid) @inline negative_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k+1, grid) +@inline function positive_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) + if func(i, j, k, grid) + if xₚ < xᶠ⁻ + return xᶠ⁻ + (xᶠ⁻ - xₚ) * restitution + end + end +end + +@inline function negative_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) + if func(i, j, k, grid) + if xₚ < xᶠ⁻ + return xᶠ⁺ - (xₚ - xᶠ⁺) * restitution + end + end +end + @inline function enforce_immersed_boundary_condition(particles, p, grid, restitution) xₚ, yₚ, zₚ = (particles.x[p], particles.y[p], particles.z[p]) i, j, k = Int.(floor(fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid))) - if positive_x_immersed_boundary(i, j, k, grid) - xᶠ⁻ = xnode(Face(), i, grid) - xₚ < xᶠ⁻ && xₚ = xᶠ⁻ + (xᶠ⁻ - xₚ) * restitution - end - if positive_y_immersed_boundary(i, j, k, grid) - yᶠ⁻ = ynode(Face(), j, grid) - yₚ < yᶠ⁻ && yₚ = yᶠ⁻ + (yᶠ⁻ - yₚ) * restitution - end - if positive_z_immersed_boundary(i, j, k, grid) - zᶠ⁻ = znode(Face(), k, grid) - zₚ < zᶠ⁻ && zₚ = zᶠ⁻ + (zᶠ⁻ - zₚ) * restitution - end - if negative_x_immersed_boundary(i, j, k, grid) - xᶠ⁺ = xnode(Face(), i + 1, grid) - xₚ > xᶠ⁺ && xₚ = xᶠ⁺ - (xₚ - xᶠ⁺) * restitution - end - if negative_y_immersed_boundary(i, j, k, grid) - yᶠ⁺ = ynode(Face(), j + 1, grid) - yₚ > yᶠ⁺ && yₚ = yᶠ⁺ - (yₚ - yᶠ⁺) * restitution - end - if negative_z_immersed_boundary(i, j, k, grid) - zᶠ⁺ = znode(Face(), k + 1, grid) - zₚ > zᶠ⁺ && zₚ = zᶠ⁺ - (zₚ - zᶠ⁺) * restitution - end + xᶠ⁻ = xnode(Face(), i, grid) + yᶠ⁻ = ynode(Face(), j, grid) + zᶠ⁻ = znode(Face(), k, grid) + xᶠ⁺ = xnode(Face(), i + 1, grid) + yᶠ⁺ = ynode(Face(), j + 1, grid) + zᶠ⁺ = znode(Face(), k + 1, grid) + + xₚ = positive_correction(xₚ, xᶠ⁻, restitution, positive_x_immersed_boundary, i, j, k, grid) + yₚ = positive_correction(yₚ, yᶠ⁻, restitution, positive_y_immersed_boundary, i, j, k, grid) + zₚ = positive_correction(zₚ, zᶠ⁻, restitution, positive_z_immersed_boundary, i, j, k, grid) + + xₚ = negative_correction(xₚ, xᶠ⁻, restitution, negative_x_immersed_boundary, i, j, k, grid) + yₚ = negative_correction(yₚ, yᶠ⁻, restitution, negative_y_immersed_boundary, i, j, k, grid) + yₚ = negative_correction(zₚ, zᶠ⁻, restitution, negative_z_immersed_boundary, i, j, k, grid) return xₚ, yₚ, zₚ end From bc9a7f7b13ac587152a66b38c9761ec5e95fd4dd Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 13:25:57 -0400 Subject: [PATCH 21/30] correct kernel for immersed boundary --- .../LagrangianParticleTracking.jl | 2 +- .../update_particle_properties.jl | 114 ++++++++---------- 2 files changed, 54 insertions(+), 62 deletions(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index ccb421bc81..dbf0b4a408 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -8,7 +8,7 @@ using KernelAbstractions using StructArrays using Oceananigans.Grids -using Oceananigans.Grids: xnode, ynode, znode +using Oceananigans.Grids: xnode, ynode, znode, AbstractUnderlyingGrid, AbstractGrid using Oceananigans.ImmersedBoundaries using Oceananigans.Architectures: device, architecture using Oceananigans.Fields: interpolate, datatuple, compute!, location, fractional_indices diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 763416eba1..36fe786b62 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -22,48 +22,6 @@ along a `Periodic` dimension, put them on the other side. return x end -@kernel function _advect_particles!(particles, restitution, grid::RectilinearGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} - p = @index(Global) - - # Advect particles using forward Euler. - @inbounds particles.x[p] += interpolate(velocities.u, Face(), Center(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - - # Enforce boundary conditions for particles. - @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], grid.xᶠᵃᵃ[1], grid.xᶠᵃᵃ[grid.Nx+1], restitution) - @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], grid.yᵃᶠᵃ[1], grid.yᵃᶠᵃ[grid.Ny+1], restitution) - @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz+1], restitution) -end - -@kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} - p = @index(Global) - - adapted_velocities = calc_correct_velocities(velocities, grid.underlying_grid) - - # Advect particles using forward Euler. - @inbounds particles.x[p] += interpolate(velocities.u, Face(), Center(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - @inbounds particles.y[p] += interpolate(velocities.v, Center(), Face(), Center(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - @inbounds particles.z[p] += interpolate(velocities.w, Center(), Center(), Face(), grid.underlying_grid, particles.x[p], particles.y[p], particles.z[p]) * Δt - - x, y, z = return_face_metrics(grid) - - # Enforce boundary conditions for particles. - @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], x[1], x[grid.Nx], restitution) - @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], y[1], y[grid.Ny], restitution) - @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], z[1], z[grid.Nz], restitution) - - # Enforce Immersed boundary conditions for particles. - @inbounds particles.x[p], particles.y[p], particles.z[p] = enforce_immersed_boundary_condition(particles, p, grid, restitution) -end - -@inline calc_correct_velocities(vel, g::RectilinearGrid) = vel -@inline calc_correct_velocities(vel, g::LatitudeLongitudeGrid) = (; u = vel.u / g.radius, v = vel.v / g.radius, w = vel.w) - -@inline return_face_metrics(g::LatitudeLongitudeGrid) = (g.λᶠᵃᵃ, g.φᵃᶠᵃ, g.zᵃᵃᶠ) -@inline return_face_metrics(g::RectilinearGrid) = (g.xᶠᵃᵃ, g.yᵃᶠᵃ, g.zᵃᵃᶠ) -@inline return_face_metrics(g::ImmersedBoundaryGrid) = return_face_metrics(g.underlying_grid) - @inline positive_x_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i-1, j, k, grid) @inline positive_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j-1, k, grid) @inline positive_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k-1, grid) @@ -72,26 +30,22 @@ end @inline negative_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j+1, k, grid) @inline negative_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k+1, grid) -@inline function positive_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) - if func(i, j, k, grid) - if xₚ < xᶠ⁻ - return xᶠ⁻ + (xᶠ⁻ - xₚ) * restitution - end - end -end +@inline positive_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) = func(i, j, k, grid) && xₚ < xᶠ ? xᶠ + (xᶠ - xₚ) * restitution : xₚ +@inline negative_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) = func(i, j, k, grid) && xₚ > xᶠ ? xᶠ - (xₚ - xᶠ) * restitution : xₚ -@inline function negative_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) - if func(i, j, k, grid) - if xₚ < xᶠ⁻ - return xᶠ⁺ - (xₚ - xᶠ⁺) * restitution - end - end -end +""" + enforce_immersed_boundary_condition(particles, p, grid, restitution) +If a particle with position `x, y, z` is at the edge of an immersed boundary, correct the +position as to avoid +""" @inline function enforce_immersed_boundary_condition(particles, p, grid, restitution) xₚ, yₚ, zₚ = (particles.x[p], particles.y[p], particles.z[p]) - i, j, k = Int.(floor(fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid))) - + i, j, k = fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid.underlying_grid) + i = Base.unsafe_trunc(Int, i) + j = Base.unsafe_trunc(Int, j) + k = Base.unsafe_trunc(Int, k) + xᶠ⁻ = xnode(Face(), i, grid) yᶠ⁻ = ynode(Face(), j, grid) zᶠ⁻ = znode(Face(), k, grid) @@ -103,13 +57,51 @@ end yₚ = positive_correction(yₚ, yᶠ⁻, restitution, positive_y_immersed_boundary, i, j, k, grid) zₚ = positive_correction(zₚ, zᶠ⁻, restitution, positive_z_immersed_boundary, i, j, k, grid) - xₚ = negative_correction(xₚ, xᶠ⁻, restitution, negative_x_immersed_boundary, i, j, k, grid) - yₚ = negative_correction(yₚ, yᶠ⁻, restitution, negative_y_immersed_boundary, i, j, k, grid) - yₚ = negative_correction(zₚ, zᶠ⁻, restitution, negative_z_immersed_boundary, i, j, k, grid) + xₚ = negative_correction(xₚ, xᶠ⁺, restitution, negative_x_immersed_boundary, i, j, k, grid) + yₚ = negative_correction(yₚ, yᶠ⁺, restitution, negative_y_immersed_boundary, i, j, k, grid) + yₚ = negative_correction(zₚ, zᶠ⁺, restitution, negative_z_immersed_boundary, i, j, k, grid) return xₚ, yₚ, zₚ end +@inline function update_particle_position!(particles, p, restitution, grid::AbstractGrid{FT, TX, TY, TZ}, Δt, velocities) where {FT, TX, TY, TZ} + + # Advect particles using forward Euler. + @inbounds u = interpolate(velocities.u, Face(), Center(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) + @inbounds v = interpolate(velocities.v, Center(), Face(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) + @inbounds w = interpolate(velocities.w, Center(), Center(), Face(), grid, particles.x[p], particles.y[p], particles.z[p]) + + @inbounds particles.x[p] += calc_correct_velocity(u, grid) * Δt + @inbounds particles.y[p] += calc_correct_velocity(v, grid) * Δt + @inbounds particles.z[p] += calc_correct_velocity(w, grid) * Δt + + x, y, z = return_face_metrics(grid) + + # Enforce boundary conditions for particles. + @inbounds particles.x[p] = enforce_boundary_conditions(TX(), particles.x[p], x[1], x[grid.Nx], restitution) + @inbounds particles.y[p] = enforce_boundary_conditions(TY(), particles.y[p], y[1], y[grid.Ny], restitution) + @inbounds particles.z[p] = enforce_boundary_conditions(TZ(), particles.z[p], z[1], z[grid.Nz], restitution) +end + +@kernel function _advect_particles!(particles, restitution, grid::AbstractUnderlyingGrid, Δt, velocities) + p = @index(Global) + update_particle_position!(particles, p, restitution, grid, Δt, velocities) +end + +@kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid, Δt, velocities) + p = @index(Global) + update_particle_position!(particles, p, restitution, grid.underlying_grid, Δt, velocities) + @inbounds particles.x[p], particles.y[p], particles.z[p] = enforce_immersed_boundary_condition(particles, p, grid, restitution) +end + +# Linear velocity for RectilinearGrid, Angular velocity for LatitudeLongitudeGrid +@inline calc_correct_velocity(U, g::RectilinearGrid) = U +@inline calc_correct_velocity(U, g::LatitudeLongitudeGrid) = U / g.radius + +@inline return_face_metrics(g::LatitudeLongitudeGrid) = (g.λᶠᵃᵃ, g.φᵃᶠᵃ, g.zᵃᵃᶠ) +@inline return_face_metrics(g::RectilinearGrid) = (g.xᶠᵃᵃ, g.yᵃᶠᵃ, g.zᵃᵃᶠ) +@inline return_face_metrics(g::ImmersedBoundaryGrid) = return_face_metrics(g.underlying_grid) + @kernel function update_field_property!(particle_property, particles, grid, field, LX, LY, LZ) p = @index(Global) From 6087bf7b8395ba99f4c2d2c06aa9e7da1312b79c Mon Sep 17 00:00:00 2001 From: simone-silvestri Date: Wed, 13 Jul 2022 15:03:22 -0400 Subject: [PATCH 22/30] fixed bunch of stuff --- .../LagrangianParticleTracking.jl | 6 ++-- .../update_particle_properties.jl | 32 ++++++++++--------- test/test_lagrangian_particle_tracking.jl | 2 +- 3 files changed, 22 insertions(+), 18 deletions(-) diff --git a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl index dbf0b4a408..84e5923f2d 100644 --- a/src/LagrangianParticleTracking/LagrangianParticleTracking.jl +++ b/src/LagrangianParticleTracking/LagrangianParticleTracking.jl @@ -8,10 +8,12 @@ using KernelAbstractions using StructArrays using Oceananigans.Grids -using Oceananigans.Grids: xnode, ynode, znode, AbstractUnderlyingGrid, AbstractGrid +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 +using Oceananigans.Fields: interpolate, datatuple, compute!, location, fractional_indices, fractional_y_index using Oceananigans.Utils: prettysummary, launch! import Base: size, length, show diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 36fe786b62..8ad1f9b0c8 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -22,13 +22,13 @@ along a `Periodic` dimension, put them on the other side. return x end -@inline positive_x_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i-1, j, k, grid) -@inline positive_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j-1, k, grid) -@inline positive_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k-1, grid) +@inline positive_x_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i-1, j, k, grid) +@inline positive_y_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j-1, k, grid) +@inline positive_z_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j, k-1, grid) -@inline negative_x_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i+1, j, k, grid) -@inline negative_y_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j+1, k, grid) -@inline negative_z_immersed_boundary(i, j, k, grid) = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k+1, grid) +@inline negative_x_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i+1, j, k, grid) +@inline negative_y_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j+1, k, grid) +@inline negative_z_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j, k+1, grid) @inline positive_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) = func(i, j, k, grid) && xₚ < xᶠ ? xᶠ + (xᶠ - xₚ) * restitution : xₚ @inline negative_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) = func(i, j, k, grid) && xₚ > xᶠ ? xᶠ - (xₚ - xᶠ) * restitution : xₚ @@ -71,9 +71,12 @@ end @inbounds v = interpolate(velocities.v, Center(), Face(), Center(), grid, particles.x[p], particles.y[p], particles.z[p]) @inbounds w = interpolate(velocities.w, Center(), Center(), Face(), grid, particles.x[p], particles.y[p], particles.z[p]) - @inbounds particles.x[p] += calc_correct_velocity(u, grid) * Δt - @inbounds particles.y[p] += calc_correct_velocity(v, grid) * Δt - @inbounds particles.z[p] += calc_correct_velocity(w, grid) * Δt + j = fractional_y_index(particles.y[p], Center(), grid) + j = Base.unsafe_trunc(Int, j) + + @inbounds particles.x[p] += calc_correct_velocity_u(u, grid, j) * Δt + @inbounds particles.y[p] += calc_correct_velocity_v(v, grid) * Δt + @inbounds particles.z[p] += w * Δt x, y, z = return_face_metrics(grid) @@ -95,8 +98,11 @@ end end # Linear velocity for RectilinearGrid, Angular velocity for LatitudeLongitudeGrid -@inline calc_correct_velocity(U, g::RectilinearGrid) = U -@inline calc_correct_velocity(U, g::LatitudeLongitudeGrid) = U / g.radius +@inline calc_correct_velocity_u(U, g::RectilinearGrid, j) = U +@inline calc_correct_velocity_u(U, g::LatitudeLongitudeGrid, j) = U / (g.radius * hack_cosd(g.φᵃᶜᵃ[j])) * 360 / (2π) + +@inline calc_correct_velocity_v(V, g::RectilinearGrid) = V +@inline calc_correct_velocity_v(V, g::LatitudeLongitudeGrid) = V / g.radius * 360 / (2π) @inline return_face_metrics(g::LatitudeLongitudeGrid) = (g.λᶠᵃᵃ, g.φᵃᶠᵃ, g.zᵃᵃᶠ) @inline return_face_metrics(g::RectilinearGrid) = (g.xᶠᵃᵃ, g.yᵃᶠᵃ, g.zᵃᵃᶠ) @@ -104,7 +110,6 @@ end @kernel function update_field_property!(particle_property, particles, grid, field, LX, LY, LZ) p = @index(Global) - @inbounds particle_property[p] = interpolate(field, LX, LY, LZ, grid, particles.x[p], particles.y[p], particles.z[p]) end @@ -139,8 +144,6 @@ function update_particle_properties!(lagrangian_particles, model, Δt) # Advect particles - - advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) advect_particles_event = advect_particles_kernel!(lagrangian_particles.properties, lagrangian_particles.restitution, model.grid, Δt, @@ -148,7 +151,6 @@ function update_particle_properties!(lagrangian_particles, model, Δt) dependencies=Event(device(arch))) wait(device(arch), advect_particles_event) - return nothing end diff --git a/test/test_lagrangian_particle_tracking.jl b/test/test_lagrangian_particle_tracking.jl index 6de8bdba6e..e88fdad4de 100644 --- a/test/test_lagrangian_particle_tracking.jl +++ b/test/test_lagrangian_particle_tracking.jl @@ -26,7 +26,7 @@ function particle_tracking_simulation(; grid, particles, timestepper=:RungeKutta nc_filepath = "test_particles.nc" sim.output_writers[:particles_nc] = - NetCDFOutputWriter(model, model.particles, filename=nc_filepath, schedule=IterationInterval(1), mode="c") + NetCDFOutputWriter(model, model.particles, filename=nc_filepath, schedule=IterationInterval(1), overwrite_existing=true) sim.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(1), dir = ".", prefix = "particles_checkpoint") From ade05a2a13b920b6364c3f75d376d020c533d1f0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 13 Jul 2022 16:55:04 -0400 Subject: [PATCH 23/30] better immersed restitution --- .../update_particle_properties.jl | 54 ++++++++++--------- 1 file changed, 28 insertions(+), 26 deletions(-) diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 8ad1f9b0c8..029a2c1e88 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -22,16 +22,12 @@ along a `Periodic` dimension, put them on the other side. return x end -@inline positive_x_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i-1, j, k, grid) -@inline positive_y_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j-1, k, grid) -@inline positive_z_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j, k-1, grid) +# Fallback if we skip ine cell +@inline adjust_coord(x, args...) where N = x -@inline negative_x_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i+1, j, k, grid) -@inline negative_y_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j+1, k, grid) -@inline negative_z_immersed_boundary(i, j, k, grid) = !immersed_cell(i, j, k, grid) & immersed_cell(i, j, k+1, grid) - -@inline positive_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) = func(i, j, k, grid) && xₚ < xᶠ ? xᶠ + (xᶠ - xₚ) * restitution : xₚ -@inline negative_correction(xₚ, xᶠ, restitution, func, i, j, k, grid) = func(i, j, k, grid) && xₚ > xᶠ ? xᶠ - (xₚ - xᶠ) * restitution : xₚ +@inline adjust_coord(x, nodefunc, i, ::Val{0} , grid, rest) = x +@inline adjust_coord(x, nodefunc, i, ::Val{-1}, grid, rest) = nodefunc(Face(), i+1, grid) - (x - nodefunc(Face(), i+1, grid)) * restitution +@inline adjust_coord(x, nodefunc, i, ::Val{1} , grid, rest) = nodefunc(Face(), i, grid) + (nodefunc(Face(), i, grid) - x) * restitution """ enforce_immersed_boundary_condition(particles, p, grid, restitution) @@ -39,27 +35,26 @@ end If a particle with position `x, y, z` is at the edge of an immersed boundary, correct the position as to avoid """ -@inline function enforce_immersed_boundary_condition(particles, p, grid, restitution) +@inline function pop_immersed_particles(particles, p, grid, restitution, old_pos) xₚ, yₚ, zₚ = (particles.x[p], particles.y[p], particles.z[p]) i, j, k = fractional_indices(xₚ, yₚ, zₚ, (Center(), Center(), Center()), grid.underlying_grid) i = Base.unsafe_trunc(Int, i) j = Base.unsafe_trunc(Int, j) k = Base.unsafe_trunc(Int, k) - - xᶠ⁻ = xnode(Face(), i, grid) - yᶠ⁻ = ynode(Face(), j, grid) - zᶠ⁻ = znode(Face(), k, grid) - xᶠ⁺ = xnode(Face(), i + 1, grid) - yᶠ⁺ = ynode(Face(), j + 1, grid) - zᶠ⁺ = znode(Face(), k + 1, grid) - - xₚ = positive_correction(xₚ, xᶠ⁻, restitution, positive_x_immersed_boundary, i, j, k, grid) - yₚ = positive_correction(yₚ, yᶠ⁻, restitution, positive_y_immersed_boundary, i, j, k, grid) - zₚ = positive_correction(zₚ, zᶠ⁻, restitution, positive_z_immersed_boundary, i, j, k, grid) - - xₚ = negative_correction(xₚ, xᶠ⁺, restitution, negative_x_immersed_boundary, i, j, k, grid) - yₚ = negative_correction(yₚ, yᶠ⁺, restitution, negative_y_immersed_boundary, i, j, k, grid) - yₚ = negative_correction(zₚ, zᶠ⁺, restitution, negative_z_immersed_boundary, i, j, k, grid) + + if immersed_cell(i, j, k, grid) + iₒ, jₒ, kₒ = fractional_indices(old_pos..., (Center(), Center(), Center()), grid.underlying_grid) + iₒ = Base.unsafe_trunc(Int, iₒ) + jₒ = Base.unsafe_trunc(Int, jₒ) + kₒ = Base.unsafe_trunc(Int, kₒ) + + if !immersed_cell(iₒ, jₒ, kₒ, grid) + iᵈ, jᵈ, kᵈ = (i, j, k) .- (iₒ, jₒ, kₒ) + xₚ = adjust_coord(xₚ, xnode, iₒ, Val(iᵈ), grid) + yₚ = adjust_coord(yₚ, ynode, jₒ, Val(jᵈ), grid) + zₚ = adjust_coord(zₚ, znode, kₒ, Val(kᵈ), grid) + end + end return xₚ, yₚ, zₚ end @@ -93,8 +88,15 @@ end @kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid, Δt, velocities) p = @index(Global) + + old_pos = (particle.x[p], particle.y[p], particle.y[p]) + update_particle_position!(particles, p, restitution, grid.underlying_grid, Δt, velocities) - @inbounds particles.x[p], particles.y[p], particles.z[p] = enforce_immersed_boundary_condition(particles, p, grid, restitution) + x, y, z = pop_immersed_particles(particles, p, grid, restitution, old_pos) + + particles.x[p] = x + particles.y[p] = y + particles.z[p] = z end # Linear velocity for RectilinearGrid, Angular velocity for LatitudeLongitudeGrid From 46867267d35b7efc1d53f5ea4cd543a6a01d4c0b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 13 Jul 2022 16:58:21 -0400 Subject: [PATCH 24/30] little bit of comment --- .../update_particle_properties.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 029a2c1e88..9a8297a170 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -22,7 +22,7 @@ along a `Periodic` dimension, put them on the other side. return x end -# Fallback if we skip ine cell +# Fallback if we skip one cell @inline adjust_coord(x, args...) where N = x @inline adjust_coord(x, nodefunc, i, ::Val{0} , grid, rest) = x @@ -30,10 +30,10 @@ end @inline adjust_coord(x, nodefunc, i, ::Val{1} , grid, rest) = nodefunc(Face(), i, grid) + (nodefunc(Face(), i, grid) - x) * restitution """ - enforce_immersed_boundary_condition(particles, p, grid, restitution) + pop_immersed_boundary_condition(particles, p, grid, restitution) -If a particle with position `x, y, z` is at the edge of an immersed boundary, correct the -position as to avoid +If a particle with position `x, y, z` is inside and immersed boundary, correct the +position based on the previous position (we bounce back a certain restitution from the old cell) """ @inline function pop_immersed_particles(particles, p, grid, restitution, old_pos) xₚ, yₚ, zₚ = (particles.x[p], particles.y[p], particles.z[p]) From 4cc17cb49e1d0d74b47bcc56d144e79bf5046dd5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 13 Jul 2022 17:32:23 -0400 Subject: [PATCH 25/30] better adjust coord --- .../update_particle_properties.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 9a8297a170..7f34df5a7b 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -26,8 +26,8 @@ end @inline adjust_coord(x, args...) where N = x @inline adjust_coord(x, nodefunc, i, ::Val{0} , grid, rest) = x -@inline adjust_coord(x, nodefunc, i, ::Val{-1}, grid, rest) = nodefunc(Face(), i+1, grid) - (x - nodefunc(Face(), i+1, grid)) * restitution -@inline adjust_coord(x, nodefunc, i, ::Val{1} , grid, rest) = nodefunc(Face(), i, grid) + (nodefunc(Face(), i, grid) - x) * restitution +@inline adjust_coord(x, nodefunc, i, ::Val{-1}, grid, rest) = nodefunc(Face(), i+1, grid) - (x - nodefunc(Face(), i+1, grid)) * rest +@inline adjust_coord(x, nodefunc, i, ::Val{1} , grid, rest) = nodefunc(Face(), i, grid) + (nodefunc(Face(), i, grid) - x) * rest """ pop_immersed_boundary_condition(particles, p, grid, restitution) @@ -89,14 +89,14 @@ end @kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid, Δt, velocities) p = @index(Global) - old_pos = (particle.x[p], particle.y[p], particle.y[p]) + # old_pos = (particles.x[p], particles.y[p], particles.y[p]) update_particle_position!(particles, p, restitution, grid.underlying_grid, Δt, velocities) - x, y, z = pop_immersed_particles(particles, p, grid, restitution, old_pos) + # x, y, z = pop_immersed_particles(particles, p, grid, restitution, old_pos) - particles.x[p] = x - particles.y[p] = y - particles.z[p] = z + # particles.x[p] = x + # particles.y[p] = y + # particles.z[p] = z end # Linear velocity for RectilinearGrid, Angular velocity for LatitudeLongitudeGrid From a4ee00e9706d793befa2dc020375252b412e1fb9 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 13 Jul 2022 18:31:49 -0400 Subject: [PATCH 26/30] fixed all tests --- .../update_particle_properties.jl | 10 +++++----- test/test_lagrangian_particle_tracking.jl | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/LagrangianParticleTracking/update_particle_properties.jl b/src/LagrangianParticleTracking/update_particle_properties.jl index 7f34df5a7b..abd631e47f 100644 --- a/src/LagrangianParticleTracking/update_particle_properties.jl +++ b/src/LagrangianParticleTracking/update_particle_properties.jl @@ -89,14 +89,14 @@ end @kernel function _advect_particles!(particles, restitution, grid::ImmersedBoundaryGrid, Δt, velocities) p = @index(Global) - # old_pos = (particles.x[p], particles.y[p], particles.y[p]) + old_pos = (particles.x[p], particles.y[p], particles.y[p]) update_particle_position!(particles, p, restitution, grid.underlying_grid, Δt, velocities) - # x, y, z = pop_immersed_particles(particles, p, grid, restitution, old_pos) + x, y, z = pop_immersed_particles(particles, p, grid, restitution, old_pos) - # particles.x[p] = x - # particles.y[p] = y - # particles.z[p] = z + particles.x[p] = x + particles.y[p] = y + particles.z[p] = z end # Linear velocity for RectilinearGrid, Angular velocity for LatitudeLongitudeGrid diff --git a/test/test_lagrangian_particle_tracking.jl b/test/test_lagrangian_particle_tracking.jl index e88fdad4de..054b65c5fe 100644 --- a/test/test_lagrangian_particle_tracking.jl +++ b/test/test_lagrangian_particle_tracking.jl @@ -112,7 +112,7 @@ function run_simple_particle_tracking_tests(arch, timestepper; vertically_stretc sim.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(1), dir = ".", prefix = "particles_checkpoint") - sim, jld2_filepath, nc_filepath = particle_tracking_simulation(; grid, particles, timestepper, velocities) + sim, jld2_filepath, nc_filepath = particle_tracking_simulation(; grid, particles=lagrangian_particles, timestepper, velocities) model = sim.model run!(sim) From 3df934f07d05d3ba0ab760470ba34bd91cbfdbd6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 14 Jul 2022 21:03:26 -0400 Subject: [PATCH 27/30] Particles are fast!!! --- src/BoundaryConditions/fill_halo_regions.jl | 1 + .../HydrostaticFreeSurfaceModels.jl | 2 +- .../hydrostatic_free_surface_ab2_step.jl | 10 ---------- .../prescribed_hydrostatic_velocity_fields.jl | 13 +++++++++++++ 4 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/BoundaryConditions/fill_halo_regions.jl b/src/BoundaryConditions/fill_halo_regions.jl index df9ece67db..12ec51358b 100644 --- a/src/BoundaryConditions/fill_halo_regions.jl +++ b/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...) diff --git a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl index 9e194982aa..9e69d470e9 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl @@ -5,7 +5,7 @@ export ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface, PrescribedVelocityFields -using KernelAbstractions: @index, @kernel, Event, MultiEvent +using KernelAbstractions: @index, @kernel, Event, MultiEvent, NoneEvent using KernelAbstractions.Extras.LoopInfo: @unroll using Oceananigans.Utils diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 34ef020c92..bd1db74e32 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -3,8 +3,6 @@ using Oceananigans.Fields: location using Oceananigans.TimeSteppers: ab2_step_field! using Oceananigans.TurbulenceClosures: implicit_step! -using KernelAbstractions: NoneEvent - import Oceananigans.TimeSteppers: ab2_step! ##### @@ -22,10 +20,6 @@ function ab2_step!(model::HydrostaticFreeSurfaceModel, Δt, χ) # waiting all the ab2 steps (velocities, free_surface and tracers to complete) @apply_regionally wait(device(model.architecture), prognostic_field_events) - - # @apply_regionally local_ab2_step!(model, Δt, χ) - # ab2_step_free_surface!(model.free_surface, model, Δt, χ, nothing) - return nothing end @@ -44,10 +38,6 @@ function local_ab2_step!(model, Δt, χ) tuple(explicit_tracer_step_events...)) return prognostic_field_events - - # wait(device(model.architecture), MultiEvent(tuple(explicit_velocity_step_events..., explicit_tracer_step_events...))) - - # return nothing end ##### diff --git a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl index 4a575aebec..98168b625b 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl @@ -4,10 +4,13 @@ using Oceananigans.Grids: Center, Face using Oceananigans.Fields: AbstractField, FunctionField +using Oceananigans.TimeSteppers: tick! +using Oceananigans.LagrangianParticleTracking: update_particle_properties! import Oceananigans.BoundaryConditions: fill_halo_regions! import Oceananigans.Models.NonhydrostaticModels: extract_boundary_conditions import Oceananigans.Utils: datatuple +import Oceananigans.TimeSteppers: time_step! using Adapt @@ -93,3 +96,13 @@ Adapt.adapt_structure(to, velocities::PrescribedVelocityFields) = Adapt.adapt(to, velocities.v), Adapt.adapt(to, velocities.w), nothing) + +# If the model only tracks particles... do nothing but that!!! +const OnlyParticleTrackingModel = HydrostaticFreeSurfaceModel{TS, E, A, S, G, T, V, B, R, F, P, U, C} where + {TS, E, A, S, G, T, V, B, R, F, P<:LagrangianParticles, U<:PrescribedVelocityFields, C<:NamedTuple{(), Tuple{}}} + +function time_step!(model::OnlyParticleTrackingModel, Δt; euler=false) + model.timestepper.previous_Δt = Δt + tick!(model.clock, Δt) + update_particle_properties!(model, Δt) +end \ No newline at end of file From a59c8abb4c69da94a37c208baae60a355c5817f0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 15 Jul 2022 09:47:25 -0400 Subject: [PATCH 28/30] changes interpolation and bugfixes --- docs/src/model_setup/lagrangian_particles.md | 14 +++++--- src/Fields/interpolate.jl | 38 +++++++++++++++----- 2 files changed, 38 insertions(+), 14 deletions(-) diff --git a/docs/src/model_setup/lagrangian_particles.md b/docs/src/model_setup/lagrangian_particles.md index a0fb6c1336..dfdfd1c646 100644 --- a/docs/src/model_setup/lagrangian_particles.md +++ b/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" diff --git a/src/Fields/interpolate.jl b/src/Fields/interpolate.jl index 3ccdb51973..f8ee1ce634 100644 --- a/src/Fields/interpolate.jl +++ b/src/Fields/interpolate.jl @@ -1,19 +1,39 @@ using Oceananigans.Grids: XRegRectilinearGrid, YRegRectilinearGrid, ZRegRectilinearGrid using CUDA: allowscalar +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 + + else + throw(error("List not sorted, unable to search value")) + end +end + @inline function fractional_index(vec, val) - allowscalar(true) - y2 = searchsortedfirst(vec, val) - y1 = searchsortedlast(vec, val) - x2 = vec[y2] - x1 = vec[y1] - allowscalar(false) + 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 From 40d07450a978f9b44c062cbc1c78149bcc3f1d1c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 15 Jul 2022 10:16:08 -0400 Subject: [PATCH 29/30] first commit --- src/Fields/field.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/Fields/field.jl b/src/Fields/field.jl index 7d975ce18c..63ccc20117 100644 --- a/src/Fields/field.jl +++ b/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), identity, 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), identity, 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), identity, &, 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), identity, |, 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)) From 6fc81508a8f2e71e5a814fb4ae6c7f9c3c2517cf Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 15 Jul 2022 10:20:21 -0400 Subject: [PATCH 30/30] fixed problem --- src/Fields/field.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Fields/field.jl b/src/Fields/field.jl index 63ccc20117..75f5e85223 100644 --- a/src/Fields/field.jl +++ b/src/Fields/field.jl @@ -510,10 +510,10 @@ const AnyReduction = typeof(Base.any!) check_version(VERSION) = VERSION.minor > 7 -initialize_reduced_field!(::SumReduction, f, r::ReducedField, c) = check_version(VERSION) ? Base.initarray!(interior(r), identity, 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), identity, 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), identity, &, 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), identity, |, true, interior(c)) : Base.initarray!(interior(r), |, true, interior(c)) +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))