Softening model

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")
                  ],

docs/src/general/gravity.md

@@ -0,0 +1,6 @@
+# [Gravity](@id gravity)
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("general", "gravity.jl")]
+```

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);

examples/n_body/n_body_system.jl

@@ -1,24 +1,55 @@
 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]
-    G                 :: ELTYPE
-    buffer            :: Nothing
+    # Kept for compatibility with n-body benchmark code that reads `system.G`.
+    G       :: ELTYPE
+    gravity :: GR
+    buffer  :: Nothing
 
-    function NBodySystem(initial_condition, G)
+    function NBodySystem(initial_condition, gravity::NewtonianGravity)
         mass = copy(initial_condition.mass)
+        gravity_ = nbody_gravity_model(gravity, eltype(mass))
+        gravitational_constant = 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,
+                                                                       gravitational_constant,
+                                                                       gravity_,
+                                                                       nothing)
     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)
+
+    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 +73,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)
@@ -61,16 +92,13 @@ 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
 
-        # 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)
+        factor = TrixiParticles.gravity_acceleration_factor(gravity, distance,
+                                                            mass[neighbor])
 
         @inbounds for i in 1:ndims(particle_system)
-            dv[i, particle] += tmp * pos_diff[i]
+            dv[i, particle] += factor * pos_diff[i]
         end
     end
 
@@ -79,6 +107,8 @@ end
 
 function energy(v_ode, u_ode, system, semi)
     (; mass) = system
+    gravity = TrixiParticles.gravity_model(system)
+    (; gravitational_constant, softening, cutoff_radius) = gravity
 
     e = zero(eltype(system))
 
@@ -96,7 +126,10 @@ 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
+                e -= gravitational_constant * mass[particle] * mass[neighbor] *
+                     TrixiParticles.gravity_potential_factor(softening, distance)
+            end
         end
     end
 

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, NoSoftening, PlummerSoftening
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection
 export nparticles, eachparticle

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")

src/general/gravity.jl

@@ -0,0 +1,189 @@
+abstract type AbstractGravityModel end
+abstract type AbstractGravitySoftening end
+
+@doc raw"""
+    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=NoSoftening()`: Gravitational softening model.
+- `cutoff_radius=Inf`: Maximum interaction distance used by pairwise interaction kernels.
+"""
+struct NewtonianGravity{ELTYPE <: Real, SOFTENING <: AbstractGravitySoftening,
+                        CUTOFF} <: AbstractGravityModel
+    gravitational_constant :: ELTYPE
+    softening              :: SOFTENING
+    cutoff_radius          :: ELTYPE
+
+    function NewtonianGravity(; gravitational_constant,
+                              softening=NoSoftening(),
+                              cutoff_radius=oftype(float(gravitational_constant), Inf))
+        softening isa AbstractGravitySoftening ||
+            throw(ArgumentError("`softening` must be a gravitational softening model"))
+
+        gravitational_constant_, _,
+        cutoff_radius_ = promote(gravitational_constant,
+                                 softening_length_for_promotion(softening,
+                                                                gravitational_constant),
+                                 cutoff_radius)
+
+        if !isfinite(gravitational_constant_) ||
+           gravitational_constant_ < zero(gravitational_constant_)
+            throw(ArgumentError("`gravitational_constant` 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
+
+        softening_ = copy_softening_model(softening, typeof(gravitational_constant_))
+
+        return new{typeof(gravitational_constant_), typeof(softening_),
+                   !isinf(cutoff_radius_)}(gravitational_constant_,
+                                           softening_,
+                                           cutoff_radius_)
+    end
+end
+
+@inline function (softening::AbstractGravitySoftening)(distance)
+    return gravity_force_factor(softening, distance)
+end
+
+@inline function (softening::AbstractGravitySoftening)(pos_diff::AbstractVector)
+    return softening(pos_diff, norm(pos_diff))
+end
+
+@inline function (softening::AbstractGravitySoftening)(pos_diff, distance)
+    if iszero(distance)
+        return zero(pos_diff)
+    end
+
+    return gravity_force_factor(softening, distance) * pos_diff
+end
+
+@inline softening_length_for_promotion(::NoSoftening,
+                                       gravitational_constant) = zero(gravitational_constant)
+@inline softening_length_for_promotion(softening,
+                                       gravitational_constant) = softening.softening_length
+
+@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_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
+
+    return inv(distance_square * sqrt(distance_square))
+end
+
+@inline function gravity_potential_factor(::NoSoftening, distance)
+    return inv(distance)
+end
+
+@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 gravity_acceleration_factor(gravity, distance, neighbor_mass) * pos_diff
+end
+
+@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, SOFTENING}
+    (; gravitational_constant, softening, cutoff_radius) = gravity
+
+    distance > cutoff_radius && return zero(distance)
+
+    return -gravitational_constant * neighbor_mass *
+           gravity_force_factor(softening, distance)
+end
+
+@inline function gravity_acceleration(::Nothing, pos_diff, distance, neighbor_mass)
+    return zero(pos_diff)
+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

test/examples/examples.jl

@@ -501,6 +501,16 @@
     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
+        end
+
         @trixi_testset "n_body/n_body_solar_system.jl" begin
             @trixi_test_nowarn trixi_include(@__MODULE__,
                                              joinpath(examples_dir(), "n_body",