implement

Author: svchb

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)
 

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}

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)

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,