From 89ba3bd231e55a6c8d96cf23095536ccb7f190cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 12 Jun 2020 00:59:31 +0300 Subject: [PATCH 1/3] Remove ODEProblem --- src/nbody_to_ode.jl | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/nbody_to_ode.jl b/src/nbody_to_ode.jl index e17fee9..955858a 100644 --- a/src/nbody_to_ode.jl +++ b/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) From bff56e25e06f8a44ae4dc8c05252e6727a696154 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 12 Jun 2020 01:00:06 +0300 Subject: [PATCH 2/3] Change default to VelocityVerlet --- src/nbody_simulation_result.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/nbody_simulation_result.jl b/src/nbody_simulation_result.jl index ca255fe..4402101 100644 --- a/src/nbody_simulation_result.jl +++ b/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 From 83ba3dec6cccf25976e14b61295f56a08fcf52ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 12 Jun 2020 01:00:51 +0300 Subject: [PATCH 3/3] Fix tests Add some explanations for tests and remove extra whitespace --- test/electrostatics_test.jl | 18 ++++++++++-------- test/gravitational_test.jl | 18 ++++++++++-------- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/test/electrostatics_test.jl b/test/electrostatics_test.jl index 65d497b..c0107ba 100644 --- a/test/electrostatics_test.jl +++ b/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,7 +82,7 @@ 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 @@ -88,14 +90,14 @@ 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 \ No newline at end of file +end diff --git a/test/gravitational_test.jl b/test/gravitational_test.jl index 6d24f1e..acf1567 100644 --- a/test/gravitational_test.jl +++ b/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 \ No newline at end of file +end