From 2733f2aee3df657c7870f71d9b86a24c63573acc Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 14:56:18 +0200 Subject: [PATCH 1/8] implement basic gravity structs --- examples/n_body/n_body_newtonian_gravity.jl | 34 +++++++++ examples/n_body/n_body_system.jl | 43 +++++++---- src/TrixiParticles.jl | 1 + src/general/general.jl | 1 + src/general/gravity.jl | 85 +++++++++++++++++++++ test/examples/examples.jl | 37 +++++++++ test/general/general.jl | 1 + test/general/gravity.jl | 60 +++++++++++++++ 8 files changed, 248 insertions(+), 14 deletions(-) create mode 100644 examples/n_body/n_body_newtonian_gravity.jl create mode 100644 src/general/gravity.jl create mode 100644 test/general/gravity.jl diff --git a/examples/n_body/n_body_newtonian_gravity.jl b/examples/n_body/n_body_newtonian_gravity.jl new file mode 100644 index 0000000000..e043401d0b --- /dev/null +++ b/examples/n_body/n_body_newtonian_gravity.jl @@ -0,0 +1,34 @@ +# ========================================================================================== +# A minimal two-body simulation using the `NewtonianGravity` force model. +# ========================================================================================== + +using TrixiParticles +using OrdinaryDiffEqSymplecticRK + +include("n_body_system.jl") + +# ========================================================================================== +# ==== Systems + +coordinates = [-0.5 0.5; + 0.0 0.0] + +velocity = [0.0 0.0; + -sqrt(0.5) sqrt(0.5)] + +masses = [1.0, 1.0] + +initial_condition = InitialCondition(; coordinates, velocity, density=1.0, mass=masses) + +gravity = NewtonianGravity(; gravitational_constant=1.0) +particle_system = NBodySystem(initial_condition, gravity) + +# ========================================================================================== +# ==== Simulation + +semi = Semidiscretization(particle_system, neighborhood_search=nothing) + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +sol = solve(ode, SymplecticEuler(), dt=0.001, save_everystep=false); diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 8ff21c950c..b14cd63255 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -1,24 +1,37 @@ using TrixiParticles using LinearAlgebra -struct NBodySystem{NDIMS, ELTYPE <: Real, IC} <: TrixiParticles.AbstractSystem{NDIMS} +struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSystem{NDIMS} initial_condition :: IC mass :: Array{ELTYPE, 1} # [particle] + gravity :: GR G :: ELTYPE buffer :: Nothing - function NBodySystem(initial_condition, G) + function NBodySystem(initial_condition, gravity::NewtonianGravity) mass = copy(initial_condition.mass) + gravitational_constant = convert(eltype(mass), gravity.gravitational_constant) new{size(initial_condition.coordinates, 1), - eltype(mass), typeof(initial_condition)}(initial_condition, mass, G, nothing) + eltype(mass), typeof(initial_condition), typeof(gravity)}(initial_condition, + mass, gravity, + gravitational_constant, + nothing) end end +function NBodySystem(initial_condition, gravitational_constant) + gravity = NewtonianGravity(; gravitational_constant) + + return NBodySystem(initial_condition, gravity) +end + TrixiParticles.timer_name(::NBodySystem) = "nbody" @inline Base.eltype(system::NBodySystem{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE +@inline TrixiParticles.gravity_model(system::NBodySystem) = system.gravity + function TrixiParticles.write_u0!(u0, system::NBodySystem) u0 .= system.initial_condition.coordinates @@ -42,15 +55,15 @@ end function TrixiParticles.compact_support(system::NBodySystem, neighbor::NBodySystem) - # There is no cutoff. All particles interact with each other. - return Inf + return TrixiParticles.gravity_model(system, neighbor).cutoff_radius end function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::NBodySystem, neighbor_system::NBodySystem, semi) - (; mass, G) = neighbor_system + (; mass) = neighbor_system + gravity = TrixiParticles.gravity_model(particle_system, neighbor_system) system_coords = TrixiParticles.current_coordinates(u_particle_system, particle_system) neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system, neighbor_system) @@ -62,15 +75,11 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, # No interaction of a particle with itself particle_system === neighbor_system && particle === neighbor && return - # Original version - # dv = -G * mass[neighbor] * pos_diff / norm(pos_diff)^3 - - # Multiplying by pos_diff later makes this slightly faster - # Multiplying by (1 / norm) is also faster than dividing by norm - tmp = -G * mass[neighbor] * (1 / distance^3) + dv_gravity = TrixiParticles.gravity_acceleration(gravity, pos_diff, distance, + mass[neighbor]) @inbounds for i in 1:ndims(particle_system) - dv[i, particle] += tmp * pos_diff[i] + dv[i, particle] += dv_gravity[i] end end @@ -79,6 +88,8 @@ end function energy(v_ode, u_ode, system, semi) (; mass) = system + gravity = TrixiParticles.gravity_model(system) + (; gravitational_constant, softening_length, cutoff_radius) = gravity e = zero(eltype(system)) @@ -96,7 +107,11 @@ function energy(v_ode, u_ode, system, semi) pos_diff = particle_coords - neighbor_coords distance = norm(pos_diff) - e -= mass[particle] * mass[neighbor] / distance + if distance <= cutoff_radius + softened_distance = sqrt(distance^2 + softening_length^2) + e -= gravitational_constant * mass[particle] * mass[neighbor] / + softened_distance + end end end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bdf70324cd..e5fd5cd4cb 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -104,6 +104,7 @@ export WindingNumberHormann, WindingNumberJacobson export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, sample_boundary, planar_geometry_to_face export SourceTermDamping +export NewtonianGravity export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection export nparticles, eachparticle diff --git a/src/general/general.jl b/src/general/general.jl index 3e9c868507..5d52996837 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -4,6 +4,7 @@ include("density_calculators.jl") include("corrections.jl") include("smoothing_kernels.jl") +include("gravity.jl") include("initial_condition.jl") include("buffer.jl") include("interpolation.jl") diff --git a/src/general/gravity.jl b/src/general/gravity.jl new file mode 100644 index 0000000000..8299e6d971 --- /dev/null +++ b/src/general/gravity.jl @@ -0,0 +1,85 @@ +abstract type AbstractGravityModel end + +@doc raw""" + NewtonianGravity(; gravitational_constant, softening_length=0, cutoff_radius=Inf) + +Model for Newtonian pairwise self-gravity. + +# Keywords +- `gravitational_constant`: Strength of the pairwise gravity interaction. +- `softening_length=0`: Distance regularization used by pairwise interaction kernels. +- `cutoff_radius=Inf`: Maximum interaction distance used by pairwise interaction kernels. +""" +struct NewtonianGravity{ELTYPE <: Real} <: AbstractGravityModel + gravitational_constant :: ELTYPE + softening_length :: ELTYPE + cutoff_radius :: ELTYPE + + function NewtonianGravity(; gravitational_constant, + softening_length=zero(gravitational_constant), + cutoff_radius=oftype(float(gravitational_constant), Inf)) + gravitational_constant_, softening_length_, cutoff_radius_ = promote(gravitational_constant, + softening_length, + cutoff_radius) + + if gravitational_constant_ < zero(gravitational_constant_) + throw(ArgumentError("`gravitational_constant` must be non-negative")) + end + + if softening_length_ < zero(softening_length_) + throw(ArgumentError("`softening_length` must be non-negative")) + end + + if cutoff_radius_ <= zero(cutoff_radius_) + throw(ArgumentError("`cutoff_radius` must be positive")) + end + + return new{typeof(gravitational_constant_)}(gravitational_constant_, + softening_length_, + cutoff_radius_) + end +end + +@inline gravity_model(system) = nothing + +@inline function gravity_model(particle_system, neighbor_system) + return gravity_model(particle_system) +end + +@inline function gravity_acceleration(gravity::NewtonianGravity, pos_diff, distance, + neighbor_mass) + (; gravitational_constant, softening_length, cutoff_radius) = gravity + + if distance > cutoff_radius || iszero(distance) + return zero(pos_diff) + end + + distance_square = distance^2 + softening_length^2 + inverse_distance_cube = inv(distance_square * sqrt(distance_square)) + + return -gravitational_constant * neighbor_mass * inverse_distance_cube * pos_diff +end + +@inline function gravity_acceleration(::Nothing, pos_diff, distance, neighbor_mass) + return zero(pos_diff) +end + +@inline function gravity_interaction!(dv_particle, gravity, particle_system, + neighbor_system, particle, neighbor, + pos_diff, distance, m_a, m_b) + return dv_particle +end + +@inline function gravity_interaction!(dv_particle, gravity::NewtonianGravity, + particle_system, neighbor_system, particle, + neighbor, pos_diff, distance, m_a, m_b) + dv_particle[] += gravity_acceleration(gravity, pos_diff, distance, m_b) + + return dv_particle +end + +@inline function gravity_interaction!(dv_particle, ::Nothing, particle_system, + neighbor_system, particle, neighbor, + pos_diff, distance, m_a, m_b) + return dv_particle +end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 03d1e5e433..4c84f74df8 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -501,6 +501,14 @@ end @testset verbose=true "N-Body" begin + @trixi_testset "n_body/n_body_newtonian_gravity.jl" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "n_body", + "n_body_newtonian_gravity.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol) == 0 + end + @trixi_testset "n_body/n_body_solar_system.jl" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "n_body", @@ -515,6 +523,35 @@ "n_body_benchmark_trixi.jl")) [ r"WARNING: Method definition interact!.*\n" ] + + v_ode = vec(velocity) + u_ode = vec(coordinates) + dv_ode = similar(v_ode) + du_ode = similar(u_ode) + p = (; semi, split_integration_data=nothing) + + # Keep the optimized benchmark path allocation-free. This is a stable proxy + # for avoiding performance regressions without adding timing assertions to CI. + filename = tempname() + try + open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.disable_debug_timings() + end + end + + TrixiParticles.kick!(dv_ode, v_ode, u_ode, p, 0.0) + TrixiParticles.drift!(du_ode, v_ode, u_ode, p, 0.0) + + @test @allocated(TrixiParticles.kick!(dv_ode, v_ode, u_ode, p, 0.0)) == 0 + @test @allocated(TrixiParticles.drift!(du_ode, v_ode, u_ode, p, 0.0)) == 0 + finally + open(filename, "w") do f + redirect_stderr(f) do + TrixiParticles.enable_debug_timings() + end + end + end end @trixi_testset "n_body/n_body_benchmark_reference.jl" begin diff --git a/test/general/general.jl b/test/general/general.jl index 2eb24d50a0..a01a39f7d5 100644 --- a/test/general/general.jl +++ b/test/general/general.jl @@ -1,5 +1,6 @@ include("initial_condition.jl") include("smoothing_kernels.jl") +include("gravity.jl") include("density_calculator.jl") include("semidiscretization.jl") include("interpolation.jl") diff --git a/test/general/gravity.jl b/test/general/gravity.jl new file mode 100644 index 0000000000..bf66edd373 --- /dev/null +++ b/test/general/gravity.jl @@ -0,0 +1,60 @@ +@trixi_testset "Gravity" begin + gravity = NewtonianGravity(; gravitational_constant=1.0, + softening_length=0.1, + cutoff_radius=2.0) + + @test gravity isa TrixiParticles.AbstractGravityModel + @test gravity.gravitational_constant == 1.0 + @test gravity.softening_length == 0.1 + @test gravity.cutoff_radius == 2.0 + + gravity_float32 = NewtonianGravity(; gravitational_constant=1.0f0) + + @test gravity_float32.gravitational_constant === 1.0f0 + @test gravity_float32.softening_length === 0.0f0 + @test gravity_float32.cutoff_radius === Inf32 + + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=-1.0) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + softening_length=-0.1) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + cutoff_radius=0.0) + + struct GravitySystem <: TrixiParticles.AbstractSystem{2} end + system = GravitySystem() + neighbor_system = GravitySystem() + + @test TrixiParticles.gravity_model(system) === nothing + @test TrixiParticles.gravity_model(system, neighbor_system) === nothing + + dv_particle = Ref(SVector(1.0, -2.0)) + returned = TrixiParticles.gravity_interaction!(dv_particle, nothing, + system, neighbor_system, + 1, 2, SVector(1.0, 0.0), + 1.0, 1.0, 2.0) + + @test returned === dv_particle + @test dv_particle[] == SVector(1.0, -2.0) + + unsoftened_gravity = NewtonianGravity(; gravitational_constant=1.0) + dv_particle = Ref(SVector(1.0, -2.0)) + returned = TrixiParticles.gravity_interaction!(dv_particle, unsoftened_gravity, + system, neighbor_system, + 1, 2, SVector(2.0, 0.0), + 2.0, 1.0, 3.0) + + @test returned === dv_particle + @test dv_particle[] == SVector(0.25, -2.0) + + softened_gravity = NewtonianGravity(; gravitational_constant=1.0, + softening_length=1.0, + cutoff_radius=2.0) + acceleration = TrixiParticles.gravity_acceleration(softened_gravity, + SVector(2.0, 0.0), + 2.0, 3.0) + + @test acceleration ≈ SVector(-6 / (5 * sqrt(5)), 0.0) + @test TrixiParticles.gravity_acceleration(softened_gravity, + SVector(3.0, 0.0), + 3.0, 3.0) == SVector(0.0, 0.0) +end From f4cdc1e27c4f9f8529f395b03b5173e5fa71c640 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 15:24:14 +0200 Subject: [PATCH 2/8] Optimize n-body gravity fast path --- examples/n_body/n_body_system.jl | 53 ++++++++++++++++++++++++++------ 1 file changed, 43 insertions(+), 10 deletions(-) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index b14cd63255..08c41619e6 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -1,11 +1,12 @@ using TrixiParticles using LinearAlgebra -struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSystem{NDIMS} +struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR, FAST} <: + TrixiParticles.AbstractSystem{NDIMS} initial_condition :: IC mass :: Array{ELTYPE, 1} # [particle] - gravity :: GR G :: ELTYPE + gravity :: GR buffer :: Nothing function NBodySystem(initial_condition, gravity::NewtonianGravity) @@ -13,10 +14,12 @@ struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSyst gravitational_constant = convert(eltype(mass), gravity.gravitational_constant) new{size(initial_condition.coordinates, 1), - eltype(mass), typeof(initial_condition), typeof(gravity)}(initial_condition, - mass, gravity, - gravitational_constant, - nothing) + eltype(mass), typeof(initial_condition), typeof(gravity), + iszero(gravity.softening_length) && isinf(gravity.cutoff_radius)}(initial_condition, + mass, + gravitational_constant, + gravity, + nothing) end end @@ -58,10 +61,40 @@ function TrixiParticles.compact_support(system::NBodySystem, return TrixiParticles.gravity_model(system, neighbor).cutoff_radius end -function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, - v_neighbor_system, u_neighbor_system, - particle_system::NBodySystem, - neighbor_system::NBodySystem, semi) +@inline function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + particle_system::NBodySystem{NDIMS, ELTYPE, IC, + GR, true}, + neighbor_system::NBodySystem, + semi) where {NDIMS, ELTYPE, IC, GR} + (; mass, G) = neighbor_system + + system_coords = TrixiParticles.current_coordinates(u_particle_system, particle_system) + neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + TrixiParticles.foreach_point_neighbor(particle_system, neighbor_system, + system_coords, neighbor_coords, + semi) do particle, neighbor, pos_diff, distance + # No interaction of a particle with itself + particle_system === neighbor_system && particle === neighbor && return + + tmp = -G * mass[neighbor] * (1 / distance^3) + + @inbounds for i in 1:ndims(particle_system) + dv[i, particle] += tmp * pos_diff[i] + end + end + + return dv +end + +@inline function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + particle_system::NBodySystem{NDIMS, ELTYPE, IC, + GR, false}, + neighbor_system::NBodySystem, + semi) where {NDIMS, ELTYPE, IC, GR} (; mass) = neighbor_system gravity = TrixiParticles.gravity_model(particle_system, neighbor_system) From a5892ef1541374d700736a53278e0e18184073ac Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 15:24:14 +0200 Subject: [PATCH 3/8] Optimize n-body gravity dispatch --- examples/n_body/n_body_system.jl | 11 ++++--- src/general/gravity.jl | 55 +++++++++++++++++++++++++++----- 2 files changed, 53 insertions(+), 13 deletions(-) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index b14cd63255..d5359a11b9 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -4,8 +4,8 @@ using LinearAlgebra struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSystem{NDIMS} initial_condition :: IC mass :: Array{ELTYPE, 1} # [particle] - gravity :: GR G :: ELTYPE + gravity :: GR buffer :: Nothing function NBodySystem(initial_condition, gravity::NewtonianGravity) @@ -14,8 +14,9 @@ struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSyst new{size(initial_condition.coordinates, 1), eltype(mass), typeof(initial_condition), typeof(gravity)}(initial_condition, - mass, gravity, + mass, gravitational_constant, + gravity, nothing) end end @@ -75,11 +76,11 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, # No interaction of a particle with itself particle_system === neighbor_system && particle === neighbor && return - dv_gravity = TrixiParticles.gravity_acceleration(gravity, pos_diff, distance, - mass[neighbor]) + factor = TrixiParticles.gravity_acceleration_factor(gravity, distance, + mass[neighbor]) @inbounds for i in 1:ndims(particle_system) - dv[i, particle] += dv_gravity[i] + dv[i, particle] += factor * pos_diff[i] end end diff --git a/src/general/gravity.jl b/src/general/gravity.jl index 8299e6d971..629ce0c242 100644 --- a/src/general/gravity.jl +++ b/src/general/gravity.jl @@ -10,7 +10,7 @@ Model for Newtonian pairwise self-gravity. - `softening_length=0`: Distance regularization used by pairwise interaction kernels. - `cutoff_radius=Inf`: Maximum interaction distance used by pairwise interaction kernels. """ -struct NewtonianGravity{ELTYPE <: Real} <: AbstractGravityModel +struct NewtonianGravity{ELTYPE <: Real, SOFTENED, CUTOFF} <: AbstractGravityModel gravitational_constant :: ELTYPE softening_length :: ELTYPE cutoff_radius :: ELTYPE @@ -34,9 +34,11 @@ struct NewtonianGravity{ELTYPE <: Real} <: AbstractGravityModel throw(ArgumentError("`cutoff_radius` must be positive")) end - return new{typeof(gravitational_constant_)}(gravitational_constant_, - softening_length_, - cutoff_radius_) + return new{typeof(gravitational_constant_), + !iszero(softening_length_), + !isinf(cutoff_radius_)}(gravitational_constant_, + softening_length_, + cutoff_radius_) end end @@ -48,16 +50,53 @@ end @inline function gravity_acceleration(gravity::NewtonianGravity, pos_diff, distance, neighbor_mass) - (; gravitational_constant, softening_length, cutoff_radius) = gravity - - if distance > cutoff_radius || iszero(distance) + if iszero(distance) return zero(pos_diff) end + return gravity_acceleration_factor(gravity, distance, neighbor_mass) * pos_diff +end + +@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, false, + false}, + distance, neighbor_mass) where {ELTYPE} + (; gravitational_constant) = gravity + + return -gravitational_constant * neighbor_mass * (1 / distance^3) +end + +@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, true, + false}, + distance, neighbor_mass) where {ELTYPE} + (; gravitational_constant, softening_length) = gravity + + distance_square = distance^2 + softening_length^2 + inverse_distance_cube = inv(distance_square * sqrt(distance_square)) + + return -gravitational_constant * neighbor_mass * inverse_distance_cube +end + +@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, false, + true}, + distance, neighbor_mass) where {ELTYPE} + (; gravitational_constant, cutoff_radius) = gravity + + distance > cutoff_radius && return zero(distance) + + return -gravitational_constant * neighbor_mass * (1 / distance^3) +end + +@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, true, + true}, + distance, neighbor_mass) where {ELTYPE} + (; gravitational_constant, softening_length, cutoff_radius) = gravity + + distance > cutoff_radius && return zero(distance) + distance_square = distance^2 + softening_length^2 inverse_distance_cube = inv(distance_square * sqrt(distance_square)) - return -gravitational_constant * neighbor_mass * inverse_distance_cube * pos_diff + return -gravitational_constant * neighbor_mass * inverse_distance_cube end @inline function gravity_acceleration(::Nothing, pos_diff, distance, neighbor_mass) From 4213802737f5251ecae166e81c5b3eef3dfa3232 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 16:56:53 +0200 Subject: [PATCH 4/8] review --- examples/n_body/n_body_system.jl | 19 ++++++++++----- src/general/gravity.jl | 29 ++++++++++------------ test/examples/examples.jl | 37 ++++++++++++++++++++++++++++ test/general/gravity.jl | 41 +++++++++++++++++++++++++++++--- 4 files changed, 101 insertions(+), 25 deletions(-) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index d5359a11b9..b8633e04fe 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -10,14 +10,20 @@ struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSyst function NBodySystem(initial_condition, gravity::NewtonianGravity) mass = copy(initial_condition.mass) - gravitational_constant = convert(eltype(mass), gravity.gravitational_constant) + mass_eltype = eltype(mass) + gravitational_constant = convert(mass_eltype, gravity.gravitational_constant) + softening_length = convert(mass_eltype, gravity.softening_length) + cutoff_radius = convert(mass_eltype, gravity.cutoff_radius) + gravity_ = NewtonianGravity(; gravitational_constant, softening_length, + cutoff_radius) + gravitational_constant = gravity_.gravitational_constant new{size(initial_condition.coordinates, 1), - eltype(mass), typeof(initial_condition), typeof(gravity)}(initial_condition, - mass, - gravitational_constant, - gravity, - nothing) + eltype(mass), typeof(initial_condition), typeof(gravity_)}(initial_condition, + mass, + gravitational_constant, + gravity_, + nothing) end end @@ -75,6 +81,7 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, semi) do particle, neighbor, pos_diff, distance # No interaction of a particle with itself particle_system === neighbor_system && particle === neighbor && return + iszero(distance) && return factor = TrixiParticles.gravity_acceleration_factor(gravity, distance, mass[neighbor]) diff --git a/src/general/gravity.jl b/src/general/gravity.jl index 629ce0c242..75942fd8ea 100644 --- a/src/general/gravity.jl +++ b/src/general/gravity.jl @@ -18,20 +18,23 @@ struct NewtonianGravity{ELTYPE <: Real, SOFTENED, CUTOFF} <: AbstractGravityMode function NewtonianGravity(; gravitational_constant, softening_length=zero(gravitational_constant), cutoff_radius=oftype(float(gravitational_constant), Inf)) - gravitational_constant_, softening_length_, cutoff_radius_ = promote(gravitational_constant, - softening_length, - cutoff_radius) - - if gravitational_constant_ < zero(gravitational_constant_) - throw(ArgumentError("`gravitational_constant` must be non-negative")) + gravitational_constant_, softening_length_, + cutoff_radius_ = promote(gravitational_constant, + softening_length, + cutoff_radius) + + if !isfinite(gravitational_constant_) || + gravitational_constant_ < zero(gravitational_constant_) + throw(ArgumentError("`gravitational_constant` must be non-negative and finite")) end - if softening_length_ < zero(softening_length_) - throw(ArgumentError("`softening_length` must be non-negative")) + if !isfinite(softening_length_) || + softening_length_ < zero(softening_length_) + throw(ArgumentError("`softening_length` must be non-negative and finite")) end - if cutoff_radius_ <= zero(cutoff_radius_) - throw(ArgumentError("`cutoff_radius` must be positive")) + if isnan(cutoff_radius_) || cutoff_radius_ <= zero(cutoff_radius_) + throw(ArgumentError("`cutoff_radius` must be positive or `Inf`")) end return new{typeof(gravitational_constant_), @@ -103,12 +106,6 @@ end return zero(pos_diff) end -@inline function gravity_interaction!(dv_particle, gravity, particle_system, - neighbor_system, particle, neighbor, - pos_diff, distance, m_a, m_b) - return dv_particle -end - @inline function gravity_interaction!(dv_particle, gravity::NewtonianGravity, particle_system, neighbor_system, particle, neighbor, pos_diff, distance, m_a, m_b) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 4c84f74df8..108948a18f 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -507,6 +507,43 @@ "n_body_newtonian_gravity.jl")) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol) == 0 + + coordinates32 = Float32[0 1; + 0 0] + velocity32 = zeros(Float32, 2, 2) + masses32 = Float32[1, 2] + initial_condition32 = InitialCondition(; coordinates=coordinates32, + velocity=velocity32, + density=1.0f0, + mass=masses32, + particle_spacing=-1.0f0) + gravity32 = NewtonianGravity(; gravitational_constant=1.0, + softening_length=0.25, + cutoff_radius=2.0) + particle_system32 = NBodySystem(initial_condition32, gravity32) + + @test particle_system32.G === 1.0f0 + @test particle_system32.gravity.gravitational_constant === 1.0f0 + @test particle_system32.gravity.softening_length === 0.25f0 + @test particle_system32.gravity.cutoff_radius === 2.0f0 + + duplicate_coordinates = zeros(Float32, 2, 2) + duplicate_ic = InitialCondition(; coordinates=duplicate_coordinates, + velocity=velocity32, + density=1.0f0, + mass=masses32, + particle_spacing=-1.0f0) + duplicate_system = NBodySystem(duplicate_ic, 1.0f0) + duplicate_semi = Semidiscretization(duplicate_system, + neighborhood_search=nothing) + duplicate_ode = semidiscretize(duplicate_semi, (0.0f0, 1.0f0)) + v_duplicate, u_duplicate = duplicate_ode.u0.x + dv_duplicate = similar(v_duplicate) + TrixiParticles.kick!(dv_duplicate, v_duplicate, u_duplicate, + (; semi=duplicate_semi, + split_integration_data=nothing), 0.0f0) + + @test all(iszero, dv_duplicate) end @trixi_testset "n_body/n_body_solar_system.jl" begin diff --git a/test/general/gravity.jl b/test/general/gravity.jl index bf66edd373..1061a10841 100644 --- a/test/general/gravity.jl +++ b/test/general/gravity.jl @@ -15,10 +15,16 @@ @test gravity_float32.cutoff_radius === Inf32 @test_throws ArgumentError NewtonianGravity(; gravitational_constant=-1.0) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=Inf) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=NaN) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, softening_length=-0.1) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + softening_length=NaN) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, cutoff_radius=0.0) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + cutoff_radius=NaN) struct GravitySystem <: TrixiParticles.AbstractSystem{2} end system = GravitySystem() @@ -36,6 +42,15 @@ @test returned === dv_particle @test dv_particle[] == SVector(1.0, -2.0) + struct UnknownGravity <: TrixiParticles.AbstractGravityModel end + + @test_throws MethodError TrixiParticles.gravity_interaction!(dv_particle, + UnknownGravity(), + system, neighbor_system, + 1, 2, + SVector(1.0, 0.0), + 1.0, 1.0, 2.0) + unsoftened_gravity = NewtonianGravity(; gravitational_constant=1.0) dv_particle = Ref(SVector(1.0, -2.0)) returned = TrixiParticles.gravity_interaction!(dv_particle, unsoftened_gravity, @@ -45,16 +60,36 @@ @test returned === dv_particle @test dv_particle[] == SVector(0.25, -2.0) + @test TrixiParticles.gravity_acceleration(unsoftened_gravity, + SVector(0.0, 0.0), + 0.0, 3.0) == SVector(0.0, 0.0) softened_gravity = NewtonianGravity(; gravitational_constant=1.0, - softening_length=1.0, - cutoff_radius=2.0) + softening_length=1.0) acceleration = TrixiParticles.gravity_acceleration(softened_gravity, SVector(2.0, 0.0), 2.0, 3.0) @test acceleration ≈ SVector(-6 / (5 * sqrt(5)), 0.0) - @test TrixiParticles.gravity_acceleration(softened_gravity, + + cutoff_gravity = NewtonianGravity(; gravitational_constant=1.0, + cutoff_radius=2.0) + @test TrixiParticles.gravity_acceleration(cutoff_gravity, + SVector(2.0, 0.0), + 2.0, 3.0) == SVector(-0.75, 0.0) + @test TrixiParticles.gravity_acceleration(cutoff_gravity, + SVector(3.0, 0.0), + 3.0, 3.0) == SVector(0.0, 0.0) + + softened_cutoff_gravity = NewtonianGravity(; gravitational_constant=1.0, + softening_length=1.0, + cutoff_radius=2.0) + acceleration = TrixiParticles.gravity_acceleration(softened_cutoff_gravity, + SVector(2.0, 0.0), + 2.0, 3.0) + + @test acceleration ≈ SVector(-6 / (5 * sqrt(5)), 0.0) + @test TrixiParticles.gravity_acceleration(softened_cutoff_gravity, SVector(3.0, 0.0), 3.0, 3.0) == SVector(0.0, 0.0) end From cfd4d8856115d981792d10fa07b249b99b63c09e Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 19:18:45 +0200 Subject: [PATCH 5/8] review --- examples/n_body/n_body_system.jl | 7 +++--- test/examples/examples.jl | 43 +++----------------------------- test/examples/n_body_system.jl | 40 +++++++++++++++++++++++++++++ 3 files changed, 48 insertions(+), 42 deletions(-) create mode 100644 test/examples/n_body_system.jl diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index b8633e04fe..cc59280df6 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -4,9 +4,10 @@ using LinearAlgebra struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSystem{NDIMS} initial_condition :: IC mass :: Array{ELTYPE, 1} # [particle] - G :: ELTYPE - gravity :: GR - buffer :: Nothing + # Kept for compatibility with n-body benchmark code that reads `system.G`. + G :: ELTYPE + gravity :: GR + buffer :: Nothing function NBodySystem(initial_condition, gravity::NewtonianGravity) mass = copy(initial_condition.mass) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 108948a18f..55a9b3b14d 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -501,49 +501,14 @@ end @testset verbose=true "N-Body" begin + include("n_body_system.jl") + @trixi_testset "n_body/n_body_newtonian_gravity.jl" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "n_body", "n_body_newtonian_gravity.jl")) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol) == 0 - - coordinates32 = Float32[0 1; - 0 0] - velocity32 = zeros(Float32, 2, 2) - masses32 = Float32[1, 2] - initial_condition32 = InitialCondition(; coordinates=coordinates32, - velocity=velocity32, - density=1.0f0, - mass=masses32, - particle_spacing=-1.0f0) - gravity32 = NewtonianGravity(; gravitational_constant=1.0, - softening_length=0.25, - cutoff_radius=2.0) - particle_system32 = NBodySystem(initial_condition32, gravity32) - - @test particle_system32.G === 1.0f0 - @test particle_system32.gravity.gravitational_constant === 1.0f0 - @test particle_system32.gravity.softening_length === 0.25f0 - @test particle_system32.gravity.cutoff_radius === 2.0f0 - - duplicate_coordinates = zeros(Float32, 2, 2) - duplicate_ic = InitialCondition(; coordinates=duplicate_coordinates, - velocity=velocity32, - density=1.0f0, - mass=masses32, - particle_spacing=-1.0f0) - duplicate_system = NBodySystem(duplicate_ic, 1.0f0) - duplicate_semi = Semidiscretization(duplicate_system, - neighborhood_search=nothing) - duplicate_ode = semidiscretize(duplicate_semi, (0.0f0, 1.0f0)) - v_duplicate, u_duplicate = duplicate_ode.u0.x - dv_duplicate = similar(v_duplicate) - TrixiParticles.kick!(dv_duplicate, v_duplicate, u_duplicate, - (; semi=duplicate_semi, - split_integration_data=nothing), 0.0f0) - - @test all(iszero, dv_duplicate) end @trixi_testset "n_body/n_body_solar_system.jl" begin @@ -567,8 +532,8 @@ du_ode = similar(u_ode) p = (; semi, split_integration_data=nothing) - # Keep the optimized benchmark path allocation-free. This is a stable proxy - # for avoiding performance regressions without adding timing assertions to CI. + # Keep the benchmark system's generic `kick!`/`drift!` path allocation-free. + # This avoids noisy timing assertions while still catching practical regressions. filename = tempname() try open(filename, "w") do f diff --git a/test/examples/n_body_system.jl b/test/examples/n_body_system.jl new file mode 100644 index 0000000000..42eb241d2f --- /dev/null +++ b/test/examples/n_body_system.jl @@ -0,0 +1,40 @@ +@trixi_testset "n_body/n_body_system.jl" begin + include(joinpath(examples_dir(), "n_body", "n_body_system.jl")) + + coordinates32 = Float32[0 1; + 0 0] + velocity32 = zeros(Float32, 2, 2) + masses32 = Float32[1, 2] + initial_condition32 = InitialCondition(; coordinates=coordinates32, + velocity=velocity32, + density=1.0f0, + mass=masses32, + particle_spacing=-1.0f0) + gravity32 = NewtonianGravity(; gravitational_constant=1.0, + softening_length=0.25, + cutoff_radius=2.0) + particle_system32 = NBodySystem(initial_condition32, gravity32) + + @test particle_system32.G === 1.0f0 + @test particle_system32.gravity.gravitational_constant === 1.0f0 + @test particle_system32.gravity.softening_length === 0.25f0 + @test particle_system32.gravity.cutoff_radius === 2.0f0 + + duplicate_coordinates = zeros(Float32, 2, 2) + duplicate_ic = InitialCondition(; coordinates=duplicate_coordinates, + velocity=velocity32, + density=1.0f0, + mass=masses32, + particle_spacing=-1.0f0) + duplicate_system = NBodySystem(duplicate_ic, 1.0f0) + duplicate_semi = Semidiscretization(duplicate_system, + neighborhood_search=nothing) + duplicate_ode = semidiscretize(duplicate_semi, (0.0f0, 1.0f0)) + v_duplicate, u_duplicate = duplicate_ode.u0.x + dv_duplicate = similar(v_duplicate) + TrixiParticles.kick!(dv_duplicate, v_duplicate, u_duplicate, + (; semi=duplicate_semi, + split_integration_data=nothing), 0.0f0) + + @test all(iszero, dv_duplicate) +end From 30630cdbb025fd2c5b9e89620764d747bde4cbe9 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 22:14:29 +0200 Subject: [PATCH 6/8] Move base gravity docs and tests --- docs/make.jl | 1 + docs/src/general/gravity.md | 6 +++++ test/examples/examples.jl | 29 --------------------- test/examples/n_body_system.jl | 46 ++++++++++++++++++++++++++++++++++ 4 files changed, 53 insertions(+), 29 deletions(-) create mode 100644 docs/src/general/gravity.md diff --git a/docs/make.jl b/docs/make.jl index 3c04970b77..c6119e09c2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -113,6 +113,7 @@ makedocs(sitename="TrixiParticles.jl", "Interpolation" => joinpath("general", "interpolation.md"), "Density Calculators" => joinpath("general", "density_calculators.md"), "Smoothing Kernels" => joinpath("general", "smoothing_kernels.md"), + "Gravity" => joinpath("general", "gravity.md"), "Neighborhood Search" => joinpath("general", "neighborhood_search.md"), "Util" => joinpath("general", "util.md") ], diff --git a/docs/src/general/gravity.md b/docs/src/general/gravity.md new file mode 100644 index 0000000000..b175cfd099 --- /dev/null +++ b/docs/src/general/gravity.md @@ -0,0 +1,6 @@ +# [Gravity](@id gravity) + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("general", "gravity.jl")] +``` diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 55a9b3b14d..2f4b21a3c2 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -525,35 +525,6 @@ "n_body_benchmark_trixi.jl")) [ r"WARNING: Method definition interact!.*\n" ] - - v_ode = vec(velocity) - u_ode = vec(coordinates) - dv_ode = similar(v_ode) - du_ode = similar(u_ode) - p = (; semi, split_integration_data=nothing) - - # Keep the benchmark system's generic `kick!`/`drift!` path allocation-free. - # This avoids noisy timing assertions while still catching practical regressions. - filename = tempname() - try - open(filename, "w") do f - redirect_stderr(f) do - TrixiParticles.disable_debug_timings() - end - end - - TrixiParticles.kick!(dv_ode, v_ode, u_ode, p, 0.0) - TrixiParticles.drift!(du_ode, v_ode, u_ode, p, 0.0) - - @test @allocated(TrixiParticles.kick!(dv_ode, v_ode, u_ode, p, 0.0)) == 0 - @test @allocated(TrixiParticles.drift!(du_ode, v_ode, u_ode, p, 0.0)) == 0 - finally - open(filename, "w") do f - redirect_stderr(f) do - TrixiParticles.enable_debug_timings() - end - end - end end @trixi_testset "n_body/n_body_benchmark_reference.jl" begin diff --git a/test/examples/n_body_system.jl b/test/examples/n_body_system.jl index 42eb241d2f..9ad76de32d 100644 --- a/test/examples/n_body_system.jl +++ b/test/examples/n_body_system.jl @@ -37,4 +37,50 @@ split_integration_data=nothing), 0.0f0) @test all(iszero, dv_duplicate) + + energy_coordinates = Float64[0 3; + 0 0] + energy_velocity = zeros(2, 2) + energy_masses = [2.0, 4.0] + energy_ic = InitialCondition(; coordinates=energy_coordinates, + velocity=energy_velocity, density=1.0, + mass=energy_masses, particle_spacing=-1.0) + softened_energy_system = NBodySystem(energy_ic, + NewtonianGravity(; gravitational_constant=5.0, + softening_length=4.0, + cutoff_radius=10.0)) + softened_energy_semi = Semidiscretization(softened_energy_system, + neighborhood_search=nothing) + softened_energy_ode = semidiscretize(softened_energy_semi, (0.0, 1.0)) + + @test energy(softened_energy_ode.u0.x..., + softened_energy_system, softened_energy_semi) ≈ -8.0 + + cutoff_energy_system = NBodySystem(energy_ic, + NewtonianGravity(; gravitational_constant=5.0, + cutoff_radius=2.0)) + cutoff_energy_semi = Semidiscretization(cutoff_energy_system, + neighborhood_search=nothing) + cutoff_energy_ode = semidiscretize(cutoff_energy_semi, (0.0, 1.0)) + + @test energy(cutoff_energy_ode.u0.x..., cutoff_energy_system, + cutoff_energy_semi) == 0.0 + + function test_nbody_rhs_allocations() + coordinates = Float64[0 1; + 0 0] + velocity = zeros(2, 2) + masses = [1.0, 2.0] + initial_condition = InitialCondition(; coordinates, velocity, density=1.0, + mass=masses, particle_spacing=-1.0) + particle_system = NBodySystem(initial_condition, 1.0) + semi = Semidiscretization(particle_system, neighborhood_search=nothing) + ode = semidiscretize(semi, (0.0, 1.0)) + + sol = solve(ode, SymplecticEuler(), dt=0.1, save_everystep=false) + + @test count_rhs_allocations(sol) == 0 + end + + test_nbody_rhs_allocations() end From f3efab281d6dc4751933420abf8d4dc5ccb84c69 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 20:15:58 +0200 Subject: [PATCH 7/8] initial implementation --- examples/n_body/n_body_system.jl | 14 +-- src/TrixiParticles.jl | 2 +- src/general/gravity.jl | 171 +++++++++++++++++++++++-------- test/examples/n_body_system.jl | 5 +- test/general/gravity.jl | 30 +++++- 5 files changed, 159 insertions(+), 63 deletions(-) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index cc59280df6..031454e421 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -11,12 +11,7 @@ struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSyst function NBodySystem(initial_condition, gravity::NewtonianGravity) mass = copy(initial_condition.mass) - mass_eltype = eltype(mass) - gravitational_constant = convert(mass_eltype, gravity.gravitational_constant) - softening_length = convert(mass_eltype, gravity.softening_length) - cutoff_radius = convert(mass_eltype, gravity.cutoff_radius) - gravity_ = NewtonianGravity(; gravitational_constant, softening_length, - cutoff_radius) + gravity_ = TrixiParticles.copy_gravity_model(gravity, eltype(mass)) gravitational_constant = gravity_.gravitational_constant new{size(initial_condition.coordinates, 1), @@ -98,7 +93,7 @@ end function energy(v_ode, u_ode, system, semi) (; mass) = system gravity = TrixiParticles.gravity_model(system) - (; gravitational_constant, softening_length, cutoff_radius) = gravity + (; gravitational_constant, softening, cutoff_radius) = gravity e = zero(eltype(system)) @@ -117,9 +112,8 @@ function energy(v_ode, u_ode, system, semi) distance = norm(pos_diff) if distance <= cutoff_radius - softened_distance = sqrt(distance^2 + softening_length^2) - e -= gravitational_constant * mass[particle] * mass[neighbor] / - softened_distance + e -= gravitational_constant * mass[particle] * mass[neighbor] * + TrixiParticles.gravity_potential_factor(softening, distance) end end end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index e5fd5cd4cb..b5c07207ae 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -104,7 +104,7 @@ export WindingNumberHormann, WindingNumberJacobson export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, sample_boundary, planar_geometry_to_face export SourceTermDamping -export NewtonianGravity +export NewtonianGravity, NoSoftening, PlummerSoftening export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection export nparticles, eachparticle diff --git a/src/general/gravity.jl b/src/general/gravity.jl index 75942fd8ea..5ff21835b2 100644 --- a/src/general/gravity.jl +++ b/src/general/gravity.jl @@ -1,26 +1,74 @@ abstract type AbstractGravityModel end +abstract type AbstractGravitySoftening end @doc raw""" - NewtonianGravity(; gravitational_constant, softening_length=0, cutoff_radius=Inf) + NoSoftening() + +Newtonian gravity without softening. + +Calling `softening(distance)` returns the scalar force factor ``1 / r^3``. +Calling `softening(pos_diff)` or `softening(pos_diff, distance)` returns the +corresponding vector factor. +""" +struct NoSoftening <: AbstractGravitySoftening end + +@doc raw""" + PlummerSoftening(softening_length) + PlummerSoftening(; softening_length) + +Plummer gravitational softening with length ``\epsilon``. + +Calling `softening(distance)` returns the scalar force factor +``1 / (r^2 + \epsilon^2)^{3/2}``. Calling `softening(pos_diff)` or +`softening(pos_diff, distance)` returns the corresponding vector factor. +""" +struct PlummerSoftening{ELTYPE <: Real} <: AbstractGravitySoftening + softening_length::ELTYPE + + function PlummerSoftening(softening_length) + if !isfinite(softening_length) || softening_length < zero(softening_length) + throw(ArgumentError("`softening_length` must be non-negative and finite")) + end + + return new{typeof(softening_length)}(softening_length) + end +end + +function PlummerSoftening(; softening_length) + return PlummerSoftening(softening_length) +end + +@doc raw""" + NewtonianGravity(; gravitational_constant, softening=NoSoftening(), cutoff_radius=Inf) Model for Newtonian pairwise self-gravity. # Keywords - `gravitational_constant`: Strength of the pairwise gravity interaction. -- `softening_length=0`: Distance regularization used by pairwise interaction kernels. +- `softening=NoSoftening()`: Gravitational softening model. - `cutoff_radius=Inf`: Maximum interaction distance used by pairwise interaction kernels. """ -struct NewtonianGravity{ELTYPE <: Real, SOFTENED, CUTOFF} <: AbstractGravityModel +struct NewtonianGravity{ELTYPE <: Real, SOFTENING <: AbstractGravitySoftening, + CUTOFF} <: AbstractGravityModel gravitational_constant :: ELTYPE - softening_length :: ELTYPE + softening :: SOFTENING cutoff_radius :: ELTYPE function NewtonianGravity(; gravitational_constant, - softening_length=zero(gravitational_constant), + softening=NoSoftening(), + softening_length=nothing, cutoff_radius=oftype(float(gravitational_constant), Inf)) - gravitational_constant_, softening_length_, + if softening_length !== nothing + softening isa NoSoftening || + throw(ArgumentError("`softening` and `softening_length` cannot both be set")) + softening = iszero(softening_length) ? NoSoftening() : + PlummerSoftening(softening_length) + end + + gravitational_constant_, _, cutoff_radius_ = promote(gravitational_constant, - softening_length, + softening_length_for_promotion(softening, + gravitational_constant), cutoff_radius) if !isfinite(gravitational_constant_) || @@ -28,78 +76,111 @@ struct NewtonianGravity{ELTYPE <: Real, SOFTENED, CUTOFF} <: AbstractGravityMode throw(ArgumentError("`gravitational_constant` must be non-negative and finite")) end - if !isfinite(softening_length_) || - softening_length_ < zero(softening_length_) - throw(ArgumentError("`softening_length` must be non-negative and finite")) - end - if isnan(cutoff_radius_) || cutoff_radius_ <= zero(cutoff_radius_) throw(ArgumentError("`cutoff_radius` must be positive or `Inf`")) end - return new{typeof(gravitational_constant_), - !iszero(softening_length_), + softening_ = copy_softening_model(softening, typeof(gravitational_constant_)) + + return new{typeof(gravitational_constant_), typeof(softening_), !isinf(cutoff_radius_)}(gravitational_constant_, - softening_length_, + softening_, cutoff_radius_) end end -@inline gravity_model(system) = nothing +@inline function (softening::AbstractGravitySoftening)(distance) + return gravity_force_factor(softening, distance) +end -@inline function gravity_model(particle_system, neighbor_system) - return gravity_model(particle_system) +@inline function (softening::AbstractGravitySoftening)(pos_diff::AbstractVector) + return softening(pos_diff, norm(pos_diff)) end -@inline function gravity_acceleration(gravity::NewtonianGravity, pos_diff, distance, - neighbor_mass) +@inline function (softening::AbstractGravitySoftening)(pos_diff, distance) if iszero(distance) return zero(pos_diff) end - return gravity_acceleration_factor(gravity, distance, neighbor_mass) * pos_diff + return gravity_force_factor(softening, distance) * pos_diff end -@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, false, - false}, - distance, neighbor_mass) where {ELTYPE} - (; gravitational_constant) = gravity +@inline softening_length_for_promotion(::NoSoftening, + gravitational_constant) = zero(gravitational_constant) +@inline softening_length_for_promotion(softening, + gravitational_constant) = softening.softening_length + +function copy_gravity_model(gravity::NewtonianGravity, ::Type{ELTYPE}) where {ELTYPE} + return NewtonianGravity(; + gravitational_constant=convert(ELTYPE, + gravity.gravitational_constant), + softening=copy_softening_model(gravity.softening, ELTYPE), + cutoff_radius=convert(ELTYPE, gravity.cutoff_radius)) +end - return -gravitational_constant * neighbor_mass * (1 / distance^3) +@inline copy_softening_model(::NoSoftening, ::Type{ELTYPE}) where {ELTYPE} = NoSoftening() +@inline function copy_softening_model(softening::PlummerSoftening, + ::Type{ELTYPE}) where {ELTYPE} + return PlummerSoftening(convert(ELTYPE, softening.softening_length)) end -@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, true, - false}, - distance, neighbor_mass) where {ELTYPE} - (; gravitational_constant, softening_length) = gravity +@inline function gravity_force_factor(::NoSoftening, distance) + return inv(distance^3) +end + +@inline function gravity_force_factor(softening::PlummerSoftening, distance) + (; softening_length) = softening distance_square = distance^2 + softening_length^2 - inverse_distance_cube = inv(distance_square * sqrt(distance_square)) - return -gravitational_constant * neighbor_mass * inverse_distance_cube + return inv(distance_square * sqrt(distance_square)) end -@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, false, - true}, - distance, neighbor_mass) where {ELTYPE} - (; gravitational_constant, cutoff_radius) = gravity +@inline function gravity_potential_factor(::NoSoftening, distance) + return inv(distance) +end - distance > cutoff_radius && return zero(distance) +@inline function gravity_potential_factor(softening::PlummerSoftening, distance) + (; softening_length) = softening + + return inv(sqrt(distance^2 + softening_length^2)) +end + +@inline gravity_model(system) = nothing + +@inline function gravity_model(particle_system, neighbor_system) + return gravity_model(particle_system) +end + +@inline function gravity_acceleration(gravity::NewtonianGravity, pos_diff, distance, + neighbor_mass) + if iszero(distance) + return zero(pos_diff) + end - return -gravitational_constant * neighbor_mass * (1 / distance^3) + return gravity_acceleration_factor(gravity, distance, neighbor_mass) * pos_diff end -@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, true, +@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, SOFTENING, + false}, + distance, + neighbor_mass) where {ELTYPE, SOFTENING} + (; gravitational_constant, softening) = gravity + + return -gravitational_constant * neighbor_mass * + gravity_force_factor(softening, distance) +end + +@inline function gravity_acceleration_factor(gravity::NewtonianGravity{ELTYPE, SOFTENING, true}, - distance, neighbor_mass) where {ELTYPE} - (; gravitational_constant, softening_length, cutoff_radius) = gravity + distance, + neighbor_mass) where {ELTYPE, SOFTENING} + (; gravitational_constant, softening, cutoff_radius) = gravity distance > cutoff_radius && return zero(distance) - distance_square = distance^2 + softening_length^2 - inverse_distance_cube = inv(distance_square * sqrt(distance_square)) - - return -gravitational_constant * neighbor_mass * inverse_distance_cube + return -gravitational_constant * neighbor_mass * + gravity_force_factor(softening, distance) end @inline function gravity_acceleration(::Nothing, pos_diff, distance, neighbor_mass) diff --git a/test/examples/n_body_system.jl b/test/examples/n_body_system.jl index 9ad76de32d..d41e01e8c6 100644 --- a/test/examples/n_body_system.jl +++ b/test/examples/n_body_system.jl @@ -11,13 +11,14 @@ mass=masses32, particle_spacing=-1.0f0) gravity32 = NewtonianGravity(; gravitational_constant=1.0, - softening_length=0.25, + softening=PlummerSoftening(0.25), cutoff_radius=2.0) particle_system32 = NBodySystem(initial_condition32, gravity32) @test particle_system32.G === 1.0f0 @test particle_system32.gravity.gravitational_constant === 1.0f0 - @test particle_system32.gravity.softening_length === 0.25f0 + @test particle_system32.gravity.softening isa PlummerSoftening + @test particle_system32.gravity.softening.softening_length === 0.25f0 @test particle_system32.gravity.cutoff_radius === 2.0f0 duplicate_coordinates = zeros(Float32, 2, 2) diff --git a/test/general/gravity.jl b/test/general/gravity.jl index 1061a10841..c1e5bc26dd 100644 --- a/test/general/gravity.jl +++ b/test/general/gravity.jl @@ -1,24 +1,44 @@ @trixi_testset "Gravity" begin gravity = NewtonianGravity(; gravitational_constant=1.0, - softening_length=0.1, + softening=PlummerSoftening(0.1), cutoff_radius=2.0) @test gravity isa TrixiParticles.AbstractGravityModel @test gravity.gravitational_constant == 1.0 - @test gravity.softening_length == 0.1 + @test gravity.softening isa PlummerSoftening + @test gravity.softening.softening_length == 0.1 @test gravity.cutoff_radius == 2.0 gravity_float32 = NewtonianGravity(; gravitational_constant=1.0f0) @test gravity_float32.gravitational_constant === 1.0f0 - @test gravity_float32.softening_length === 0.0f0 + @test gravity_float32.softening isa NoSoftening @test gravity_float32.cutoff_radius === Inf32 + gravity_softening_length = NewtonianGravity(; gravitational_constant=1.0, + softening_length=0.1) + @test gravity_softening_length.softening isa PlummerSoftening + @test gravity_softening_length.softening.softening_length == 0.1 + + @test NoSoftening()(2.0) == 0.125 + @test NoSoftening()(SVector(2.0, 0.0)) == SVector(0.25, 0.0) + @test NoSoftening()(SVector(2.0, 0.0), 2.0) == SVector(0.25, 0.0) + @test PlummerSoftening(1.0)(2.0) ≈ inv(5 * sqrt(5)) + @test PlummerSoftening(1.0)(SVector(2.0, 0.0)) ≈ + SVector(2 / (5 * sqrt(5)), 0.0) + @test PlummerSoftening(; softening_length=1.0)(SVector(2.0, 0.0), 2.0) ≈ + SVector(2 / (5 * sqrt(5)), 0.0) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=-1.0) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=Inf) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=NaN) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, softening_length=-0.1) + @test_throws ArgumentError PlummerSoftening(-0.1) + @test_throws ArgumentError PlummerSoftening(Inf) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + softening=PlummerSoftening(0.1), + softening_length=0.1) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, softening_length=NaN) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, @@ -65,7 +85,7 @@ 0.0, 3.0) == SVector(0.0, 0.0) softened_gravity = NewtonianGravity(; gravitational_constant=1.0, - softening_length=1.0) + softening=PlummerSoftening(1.0)) acceleration = TrixiParticles.gravity_acceleration(softened_gravity, SVector(2.0, 0.0), 2.0, 3.0) @@ -82,7 +102,7 @@ 3.0, 3.0) == SVector(0.0, 0.0) softened_cutoff_gravity = NewtonianGravity(; gravitational_constant=1.0, - softening_length=1.0, + softening=PlummerSoftening(1.0), cutoff_radius=2.0) acceleration = TrixiParticles.gravity_acceleration(softened_cutoff_gravity, SVector(2.0, 0.0), From 15e2e8d6c26c8c0c3ac9019e2b18417abc9ac3c7 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 22 May 2026 21:41:36 +0200 Subject: [PATCH 8/8] implement --- examples/n_body/n_body_system.jl | 17 ++++++++++++++++- src/general/gravity.jl | 17 ++--------------- test/examples/n_body_system.jl | 2 +- test/general/gravity.jl | 14 +++----------- 4 files changed, 22 insertions(+), 28 deletions(-) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 031454e421..f70b8308cd 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -11,7 +11,7 @@ struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSyst function NBodySystem(initial_condition, gravity::NewtonianGravity) mass = copy(initial_condition.mass) - gravity_ = TrixiParticles.copy_gravity_model(gravity, eltype(mass)) + gravity_ = nbody_gravity_model(gravity, eltype(mass)) gravitational_constant = gravity_.gravitational_constant new{size(initial_condition.coordinates, 1), @@ -23,6 +23,21 @@ struct NBodySystem{NDIMS, ELTYPE <: Real, IC, GR} <: TrixiParticles.AbstractSyst end end +function nbody_gravity_model(gravity::NewtonianGravity, ::Type{ELTYPE}) where {ELTYPE} + return NewtonianGravity(; + gravitational_constant=convert(ELTYPE, + gravity.gravitational_constant), + softening=nbody_gravity_softening(gravity.softening, ELTYPE), + cutoff_radius=convert(ELTYPE, gravity.cutoff_radius)) +end + +@inline nbody_gravity_softening(::NoSoftening, + ::Type{ELTYPE}) where {ELTYPE} = NoSoftening() +@inline function nbody_gravity_softening(softening::PlummerSoftening, + ::Type{ELTYPE}) where {ELTYPE} + return PlummerSoftening(convert(ELTYPE, softening.softening_length)) +end + function NBodySystem(initial_condition, gravitational_constant) gravity = NewtonianGravity(; gravitational_constant) diff --git a/src/general/gravity.jl b/src/general/gravity.jl index 5ff21835b2..2d58ced644 100644 --- a/src/general/gravity.jl +++ b/src/general/gravity.jl @@ -56,14 +56,9 @@ struct NewtonianGravity{ELTYPE <: Real, SOFTENING <: AbstractGravitySoftening, function NewtonianGravity(; gravitational_constant, softening=NoSoftening(), - softening_length=nothing, cutoff_radius=oftype(float(gravitational_constant), Inf)) - if softening_length !== nothing - softening isa NoSoftening || - throw(ArgumentError("`softening` and `softening_length` cannot both be set")) - softening = iszero(softening_length) ? NoSoftening() : - PlummerSoftening(softening_length) - end + softening isa AbstractGravitySoftening || + throw(ArgumentError("`softening` must be a gravitational softening model")) gravitational_constant_, _, cutoff_radius_ = promote(gravitational_constant, @@ -110,14 +105,6 @@ end @inline softening_length_for_promotion(softening, gravitational_constant) = softening.softening_length -function copy_gravity_model(gravity::NewtonianGravity, ::Type{ELTYPE}) where {ELTYPE} - return NewtonianGravity(; - gravitational_constant=convert(ELTYPE, - gravity.gravitational_constant), - softening=copy_softening_model(gravity.softening, ELTYPE), - cutoff_radius=convert(ELTYPE, gravity.cutoff_radius)) -end - @inline copy_softening_model(::NoSoftening, ::Type{ELTYPE}) where {ELTYPE} = NoSoftening() @inline function copy_softening_model(softening::PlummerSoftening, ::Type{ELTYPE}) where {ELTYPE} diff --git a/test/examples/n_body_system.jl b/test/examples/n_body_system.jl index d41e01e8c6..0afca8d2c3 100644 --- a/test/examples/n_body_system.jl +++ b/test/examples/n_body_system.jl @@ -48,7 +48,7 @@ mass=energy_masses, particle_spacing=-1.0) softened_energy_system = NBodySystem(energy_ic, NewtonianGravity(; gravitational_constant=5.0, - softening_length=4.0, + softening=PlummerSoftening(4.0), cutoff_radius=10.0)) softened_energy_semi = Semidiscretization(softened_energy_system, neighborhood_search=nothing) diff --git a/test/general/gravity.jl b/test/general/gravity.jl index c1e5bc26dd..a03e5bf048 100644 --- a/test/general/gravity.jl +++ b/test/general/gravity.jl @@ -15,11 +15,6 @@ @test gravity_float32.softening isa NoSoftening @test gravity_float32.cutoff_radius === Inf32 - gravity_softening_length = NewtonianGravity(; gravitational_constant=1.0, - softening_length=0.1) - @test gravity_softening_length.softening isa PlummerSoftening - @test gravity_softening_length.softening.softening_length == 0.1 - @test NoSoftening()(2.0) == 0.125 @test NoSoftening()(SVector(2.0, 0.0)) == SVector(0.25, 0.0) @test NoSoftening()(SVector(2.0, 0.0), 2.0) == SVector(0.25, 0.0) @@ -33,14 +28,11 @@ @test_throws ArgumentError NewtonianGravity(; gravitational_constant=Inf) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=NaN) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, - softening_length=-0.1) + softening=nothing) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + softening=0.1) @test_throws ArgumentError PlummerSoftening(-0.1) @test_throws ArgumentError PlummerSoftening(Inf) - @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, - softening=PlummerSoftening(0.1), - softening_length=0.1) - @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, - softening_length=NaN) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, cutoff_radius=0.0) @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0,