Remove ODEProblem

src/nbody_simulation_result.jl

@@ -289,14 +289,11 @@ function run_simulation(s::NBodySimulation, args...; kwargs...)
     end
 end
 
-function calculate_simulation(s::NBodySimulation, alg_type=Tsit5(), args...; kwargs...)
-    solution = solve(ODEProblem(s), alg_type, args...; kwargs...)
-    return SimulationResult(solution, s)
-end
-
-# this should be a method for integrators designed for the SecondOrderODEProblem (It is worth somehow to sort them from other algorithms)
-function calculate_simulation(s::NBodySimulation, alg_type::Union{VelocityVerlet,DPRKN6,Yoshida6}, args...; kwargs...)
+function calculate_simulation(s::NBodySimulation, alg_type=VelocityVerlet(), args...; kwargs...)
     cb = obtain_callbacks_for_so_ode_problem(s)
+    if !haskey(kwargs, :dt)
+        kwargs = (kwargs..., dt=(s.tspan[2]-s.tspan[1])/10000)
+    end
     solution = solve(SecondOrderODEProblem(s), alg_type, args...; callback=cb, kwargs...)
     return SimulationResult(solution, s)
 end

src/nbody_to_ode.jl

@@ -291,26 +291,6 @@ function obtain_data_for_nosehoover_thermostating(simulation::NBodySimulation{<:
     (ms, kb, n, nc, γind, p)
 end
 
-function DiffEqBase.ODEProblem(simulation::NBodySimulation{<:PotentialNBodySystem})
-    (u0, v0, n) = gather_bodies_initial_coordinates(simulation)
-
-    acceleration_functions = gather_accelerations_for_potentials(simulation)
-
-    function ode_system!(du, u, p, t)
-        du[:, 1:n] = @view u[:, n + 1:2n];
-
-        @inbounds for i = 1:n
-            a = MVector(0.0, 0.0, 0.0)
-            for acceleration! in acceleration_functions
-                acceleration!(a, u[:, 1:n], u[:, n + 1:end], t, i);
-            end
-            du[:, n + i] .= a
-        end
-    end
-
-    return ODEProblem(ode_system!, hcat(u0, v0), simulation.tspan)
-end
-
 function DiffEqBase.SecondOrderODEProblem(simulation::NBodySimulation{<:PotentialNBodySystem})
 
     (u0, v0, n) = gather_bodies_initial_coordinates(simulation)

test/electrostatics_test.jl

@@ -17,11 +17,13 @@
         simulation = NBodySimulation(system, (0.0, t))
         sim_result = run_simulation(simulation)
 
-
+        # we can assume that the massive body does not move and that the second
+        # has a circular trajectory.
+        # test that the bodies return to their initial positions after one period
         solution = sim_result.solution;
         ε = 0.1 * r
-        for j = 1:2, i = 1:3
-            @test solution[1][i,j] ≈ solution[end][i,j] atol = ε
+        for i = 1:3*2
+            @test solution[1][2,i] ≈ solution[end][2,i] atol = ε
         end
 
         (qs_act, ms_act, indxs_act, exclude_act) = NBodySimulator.obtain_data_for_electrostatic_interaction(simulation.system)
@@ -80,22 +82,22 @@
         q = 1.0
         count = 1
         dL = L / (ceil(n^(1 / 3)) + 1)
-        for x = dL / 2:dL:L, y = dL / 2:dL:L, z = dL / 2:dL:L  
+        for x = dL / 2:dL:L, y = dL / 2:dL:L, z = dL / 2:dL:L
             if count > n
                 break
             end
             r = SVector(x, y, z)
             v = SVector(.0, .0, .0)
             body = ChargedParticle(r, v, m, q)
             push!(bodies, body)
-            count += 1           
+            count += 1
         end
 
         k = 9e9
         τ = 0.01 * dL / sqrt(2 * k * q * q / (dL * m))
         t1 = 0.0
         t2 = 1000 * τ
-        
+
         potential = ElectrostaticParameters(k, 0.45 * L)
         system = PotentialNBodySystem(bodies, Dict(:electrostatic => potential))
         pbc = CubicPeriodicBoundaryConditions(L)
@@ -106,6 +108,6 @@
         e_tot_2 = total_energy(result, t2)
 
         ε = 0.001
-        @test (e_tot_2 - e_tot_1) / e_tot_1 ≈ 0.0 atol = ε 
+        @test (e_tot_2 - e_tot_1) / e_tot_1 ≈ 0.0 atol = ε
     end
-end
+end

test/gravitational_test.jl

@@ -1,4 +1,4 @@
-using OrdinaryDiffEq 
+using OrdinaryDiffEq
 
 @testset "Gravitational Functional Test" begin
     G = 1
@@ -14,9 +14,10 @@ using OrdinaryDiffEq
         simulation = NBodySimulation(system, tspan)
         sim_result = run_simulation(simulation)
         solution_simo_3 = sim_result.solution;
+        # test that the bodies return to their initial positions after one period
         ε = 0.1
-        for j = 1:3, i = 1:3
-            @test solution_simo_3[1][i,j] ≈ solution_simo_3[end][i,j] atol = ε
+        for i = 1:3*3
+            @test solution_simo_3[1][2,i] ≈ solution_simo_3[end][2,i] atol = ε
         end
 
         @testset "Analyzing simulation result" begin
@@ -35,7 +36,7 @@ using OrdinaryDiffEq
             e_kin = 1.218
             @test e_kin ≈ kinetic_energy(sim_result, t1) atol = ε
         end
-    
+
 
         @testset "Using convertion into SecondOrderODEProblem" begin
             sim_result = run_simulation(simulation, DPRKN6())
@@ -75,20 +76,21 @@ using OrdinaryDiffEq
         sim_result = run_simulation(simulation, Tsit5(), abstol=1e-10, reltol=1e-10)
         solution_simo_5 = sim_result.solution;
 
+        # test that the bodies return to their initial positions after one period
         ε = 0.01
-        for j = 1:5, i = 1:3
-            @test solution_simo_5[1][i,j] ≈ solution_simo_5[end][i,j] atol = ε
+        for i = 1:3*5
+            @test solution_simo_5[1][2,i] ≈ solution_simo_5[end][2,i] atol = ε
         end
     end
 
     @testset "Constructing electorstatic potential parameters entity" begin
         default_potential = GravitationalParameters()
         @test 6.67408e-11 == default_potential.G
-    
+
         io = IOBuffer()
         potential1 = GravitationalParameters()
         potential2 = GravitationalParameters(35.67)
         @test sprint(io -> show(io, potential1)) == "Gravitational:\n\tG:6.67408e-11\n"
         @test sprint(io -> show(io, potential2)) == "Gravitational:\n\tG:35.67\n"
     end
-end
+end