type stable function calling

From the timings I'm not sure if it's worth it, but it does fix the dynamic dispatching.

```julia
using LinearAlgebra
using NBodySimulator
using Plots
using StaticArrays

J_per_eV = 1.602176634e-19
kg_per_amu = 1.66054e-27

kb = 1.380649e-23 # J/K

N = 864
m = 6.6335209e-26 # kg

box_size = 3.47786e-9 # m

σ = 0.34e-9 # m
cutoff = 0.765e-9 # m
ϵ = 1.657e-21 # J

reference_temp = 94.4 # K
mean_v = √(kb * reference_temp / m) # m/s
thermostat_prob = 0.1

eq_steps = 2000
prod_steps = 5000
Δt = 1e-14 # s

stride = 10

function get_final_bodies(sr::NBodySimulator.SimulationResult)
  positions = get_position(sr, sr.solution.t[end])
  velocities = get_velocity(sr, sr.solution.t[end])
  masses = get_masses(sr.simulation.system)
  N = length(masses)
  bodies = MassBody[]
  for i ∈ 1:N
      push!(bodies, MassBody(SVector{3}(positions[:, i]), SVector{3}(velocities[:, i]), masses[i]))
  end
  return bodies
end

initial_bodies = generate_bodies_in_cell_nodes(N, m, mean_v, box_size)
potentials = Dict(:lennard_jones => LennardJonesParameters(ϵ, σ, cutoff))
eq_system = PotentialNBodySystem(initial_bodies, potentials)

boundary_conditions = CubicPeriodicBoundaryConditions(box_size)
eq_thermostat = AndersenThermostat(reference_temp, thermostat_prob / Δt)
eq_simulation = NBodySimulation(
  eq_system,
  (0.0, eq_steps * Δt),
  boundary_conditions,
  eq_thermostat,
  kb
)

simulator = VelocityVerlet()
eq_result = run_simulation(eq_simulation, simulator, dt=Δt)
eq_bodies = get_final_bodies(eq_result)
prod_system = PotentialNBodySystem(eq_bodies, potentials)

prod_simulation = NBodySimulation(
  prod_system,
  (eq_steps * Δt, (eq_steps + prod_steps) * Δt),
  boundary_conditions,
  kb
)
@time prod_result = run_simulation(prod_simulation, simulator, dt=Δt);
prod_bodies = get_final_bodies(prod_result)
eq_time_range = [t * 1e12 for (i,t) ∈ enumerate(eq_result.solution.t) if (i - 1) % stride == 0]
prod_time_range = [t * 1e12 for (i,t) ∈ enumerate(prod_result.solution.t) if (i - 1) % stride == 0]
total_time_range = vcat(eq_time_range, prod_time_range)
plot(
  title="Temperature during Simulation",
  xlab="Time [ps]",
  ylab="Temperature [K]",
)
plot!(
  eq_time_range,
  t -> temperature(eq_result, t * 1e-12),
  label="Equilibrium Stage Temperature",
  color=1
)
plot!(
  prod_time_range,
  t -> temperature(prod_result, t * 1e-12),
  label="Production Stage Temperature",
  color=3
)
plot!(
  total_time_range,
  t -> reference_temp,
  label="Reference Temperature",
  color=2,
  linestyle=:dash
)
vline!(
  [eq_time_range[end]],
  label=false,
  color=:black,
  linestyle=:dash
)
plot(
  title="Energy during Simulation",
  xlab="Time [ps]",
  ylab="Energy [eV]",
  legend=:right
)
plot!(
  eq_time_range,
  t -> kinetic_energy(eq_result, t * 1e-12) / J_per_eV,
  label="Kinetic Energy",
  color=2
)
plot!(
  eq_time_range,
  t -> potential_energy(eq_result, t * 1e-12) / J_per_eV,
  label="Potential Energy",
  color=1
)
plot!(
  eq_time_range,
  t -> total_energy(eq_result, t * 1e-12) / J_per_eV,
  label="Total Energy",
  color=3
)
plot!(
  prod_time_range,
  t -> kinetic_energy(prod_result, t * 1e-12) / J_per_eV,
  label=false,
  color=2
)
plot!(
  prod_time_range,
  t -> potential_energy(prod_result, t * 1e-12) / J_per_eV,
  label=false,
  color=1
)
plot!(
  prod_time_range,
  t -> total_energy(prod_result, t * 1e-12) / J_per_eV,
  label=false,
  color=3
)
vline!(
  [eq_time_range[end]],
  label=false,
  color=:black,
  linestyle=:dash
)

using Profile
@Profile prod_result = run_simulation(prod_simulation, simulator, dt=Δt);
Juno.profiler()
Profile.clear()
```

Author: ChrisRackauckas

src/nbody_to_ode.jl

@@ -296,15 +296,17 @@ function DiffEqBase.ODEProblem(simulation::NBodySimulation{<:PotentialNBodySyste
 
     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);
+    ode_system! = let acceleration_functions = tuple(acceleration_functions...)
+        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
-            du[:, n + i] .= a
         end
     end
 
@@ -318,19 +320,20 @@ function DiffEqBase.SecondOrderODEProblem(simulation::NBodySimulation{<:Potentia
     acceleration_functions = gather_accelerations_for_potentials(simulation)
     simultaneous_acceleration = gather_simultaneous_acceleration(simulation)
 
-    function soode_system!(dv, v, u, p, t)
-        @inbounds for i = 1:n
-            a = MVector(0.0, 0.0, 0.0)
-            for acceleration! in acceleration_functions
-                acceleration!(a, u, v, t, i);
+    soode_system! = let acceleration_functions = tuple(acceleration_functions...), simultaneous_acceleration = tuple(simultaneous_acceleration...)
+        function soode_system!(dv, v, u, p, t)
+            @inbounds for i = 1:n
+                a = MVector(0.0, 0.0, 0.0)
+                for acceleration! in acceleration_functions
+                    acceleration!(a, u, v, t, i);
+                end
+                dv[:, i] .= a
+            end
+            for acceleration! in simultaneous_acceleration
+                acceleration!(dv, u, v, t);
             end
-            dv[:, i] .= a
-        end
-        for acceleration! in simultaneous_acceleration
-            acceleration!(dv, u, v, t);
         end
     end
-
     SecondOrderODEProblem(soode_system!, v0, u0, simulation.tspan)
 end
 
@@ -391,17 +394,19 @@ function DiffEqBase.SDEProblem(simulation::NBodySimulation{<:PotentialNBodySyste
 
     therm = simulation.thermostat
 
-    function deterministic_acceleration!(du, u, p, t)
-        du[:, 1:n] = @view u[:, n + 1:2n];
+    deterministic_acceleration! = let acceleration_functions = tuple(acceleration_functions...)
+        function deterministic_acceleration!(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, (@view u[:, 1:n]), (@view u[:, n + 1:end]), t, i);
+            @inbounds for i = 1:n
+                a = MVector(0.0, 0.0, 0.0)
+                for acceleration! in acceleration_functions
+                    acceleration!(a, (@view u[:, 1:n]), (@view u[:, n + 1:end]), t, i);
+                end
+                du[:, n + i] .= a
             end
-            du[:, n + i] .= a
+            @. du[:, n + 1:end] -= therm.γ*u[:, n + 1:end]
         end
-        @. du[:, n + 1:end] -= therm.γ*u[:, n + 1:end]
     end
 
     σ = sqrt(2*therm.γ*simulation.kb*therm.T/simulation.system.bodies[1].m)