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/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..f70b8308cd 100644 --- a/examples/n_body/n_body_system.jl +++ b/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 diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bdf70324cd..b5c07207ae 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, NoSoftening, PlummerSoftening 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..2d58ced644 --- /dev/null +++ b/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 diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 03d1e5e433..2f4b21a3c2 100644 --- a/test/examples/examples.jl +++ b/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", diff --git a/test/examples/n_body_system.jl b/test/examples/n_body_system.jl new file mode 100644 index 0000000000..0afca8d2c3 --- /dev/null +++ b/test/examples/n_body_system.jl @@ -0,0 +1,87 @@ +@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=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 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) + 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) + + 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=PlummerSoftening(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 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..a03e5bf048 --- /dev/null +++ b/test/general/gravity.jl @@ -0,0 +1,107 @@ +@trixi_testset "Gravity" begin + gravity = NewtonianGravity(; gravitational_constant=1.0, + softening=PlummerSoftening(0.1), + cutoff_radius=2.0) + + @test gravity isa TrixiParticles.AbstractGravityModel + @test gravity.gravitational_constant == 1.0 + @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 isa NoSoftening + @test gravity_float32.cutoff_radius === Inf32 + + @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=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, + cutoff_radius=0.0) + @test_throws ArgumentError NewtonianGravity(; gravitational_constant=1.0, + cutoff_radius=NaN) + + 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) + + 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, + 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) + @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=PlummerSoftening(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) + + 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=PlummerSoftening(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