mockingbirdnest · Aug 27, 2017 · Aug 28, 2017 · Aug 28, 2017 · Aug 31, 2017
Showing with 94 additions and 4 deletions.

+90 −0 benchmarks/ephemeris.cpp

+4 −4 physics/solar_system_body.hpp
diff --git a/benchmarks/ephemeris.cpp b/benchmarks/ephemeris.cpp
@@ -3,15 +3,18 @@
 
 #include <cmath>
 #include <limits>
+#include <list>
 #include <memory>
 #include <string>
 #include <vector>
 
 #include "astronomy/frames.hpp"
 #include "base/not_null.hpp"
+#include "base/thread_pool.hpp"
 #include "geometry/named_quantities.hpp"
 #include "geometry/quaternion.hpp"
 #include "geometry/rotation.hpp"
+#include "integrators/integrators.hpp"
 #include "integrators/embedded_explicit_runge_kutta_nyström_integrator.hpp"
 #include "integrators/symmetric_linear_multistep_integrator.hpp"
 #include "integrators/symplectic_runge_kutta_nyström_integrator.hpp"
@@ -36,28 +39,33 @@ using astronomy::ICRFJ2000Ecliptic;
 using astronomy::ICRFJ2000Equator;
 using astronomy::ICRFJ200EquatorialToEcliptic;
 using base::not_null;
+using base::ThreadPool;
 using geometry::Displacement;
 using geometry::Instant;
 using geometry::Position;
 using geometry::Quaternion;
 using geometry::Rotation;
 using geometry::Velocity;
+using integrators::Integrator;
 using integrators::DormandElMikkawyPrince1986RKN434FM;
 using integrators::McLachlanAtela1992Order5Optimal;
 using integrators::Quinlan1999Order8A;
 using integrators::QuinlanTremaine1990Order12;
 using quantities::DebugString;
+using quantities::Frequency;
 using quantities::Length;
 using quantities::Speed;
 using quantities::Sqrt;
 using quantities::Time;
 using quantities::astronomy::JulianYear;
 using quantities::bipm::NauticalMile;
 using quantities::si::AstronomicalUnit;
+using quantities::si::Hertz;
 using quantities::si::Kilo;
 using quantities::si::Metre;
 using quantities::si::Milli;
 using quantities::si::Minute;
+using quantities::si::Radian;
 using quantities::si::Second;
 using testing_utilities::SolarSystemFactory;
 
@@ -266,6 +274,82 @@ void EphemerisLEOProbeBenchmark(SolarSystemFactory::Accuracy const accuracy,
                  " nmi");
 }
 
+void BM_EphemerisMultithreadingBenchmark(benchmark::State& state) {
+  auto const at_спутник_1_launch =
+      SolarSystemFactory::AtСпутник1Launch(
+          SolarSystemFactory::Accuracy::AllBodiesAndOblateness);
+  Instant const epoch = at_спутник_1_launch->epoch();
+  auto const ephemeris =
+      at_спутник_1_launch->MakeEphemeris(1 * Milli(Metre),
+                                         EphemerisParameters());
+  std::string const& earth_name =
+      SolarSystemFactory::name(SolarSystemFactory::Earth);
+  auto const earth_massive_body =
+      at_спутник_1_launch->massive_body(*ephemeris, earth_name);
+  auto const earth_degrees_of_freedom =
+      at_спутник_1_launch->degrees_of_freedom(earth_name);
+
+  MasslessBody probe;
+  std::list<DiscreteTrajectory<ICRFJ2000Equator>> trajectories;
+  for (int i = 0; i < state.range_x(); ++i) {
+    KeplerianElements<ICRFJ2000Equator> elements;
+    elements.eccentricity = 0;
+    elements.semimajor_axis = (i + 1) * 10'000 * Kilo(Metre);
+    elements.inclination = 0 * Radian;
+    elements.longitude_of_ascending_node = 0 * Radian;
+    elements.argument_of_periapsis = 0 * Radian;
+    elements.true_anomaly = 0 * Radian;
+    KeplerOrbit<ICRFJ2000Equator> const orbit(
+        *earth_massive_body, probe, elements, epoch);
+    trajectories.emplace_back();
+    auto& trajectory = trajectories.back();
+    trajectory.Append(epoch,
+                      earth_degrees_of_freedom + orbit.StateVectors(epoch));
+  }
+
+  ThreadPool<void> pool(/*pool_size=*/state.range_y());
+  static constexpr int warp_factor = 6E6;
+  static constexpr Frequency refresh_frequency = 50 * Hertz;
+  static constexpr Time step = warp_factor / refresh_frequency;
+  Instant final_time = epoch;
+  while (state.KeepRunning()) {
+    state.PauseTiming();
+    std::vector<not_null<std::unique_ptr<Integrator<
+        Ephemeris<ICRFJ2000Equator>::NewtonianMotionEquation>::Instance>>>
+        instances;
+    for (auto& trajectory : trajectories) {
+      instances.push_back(ephemeris->NewInstance(
+          {&trajectory},
+          Ephemeris<ICRFJ2000Equator>::NoIntrinsicAccelerations,
+          Ephemeris<ICRFJ2000Equator>::FixedStepParameters(
+              Quinlan1999Order8A<Position<ICRFJ2000Equator>>(),
+              /*step=*/10 * Second)));
+    }
+    final_time += step;
+    state.ResumeTiming();
+
+    std::vector<std::future<void>> futures;
+    for (auto& instance : instances) {
+      futures.push_back(pool.Add([&ephemeris, &instance, final_time]() {
+        ephemeris->FlowWithFixedStep(final_time, *instance);
+      }));
+    }
+    for (auto const& future : futures) {
+      future.wait();
+    }
+  }
+
+  std::stringstream ss;
+  for (auto const& trajectory : trajectories) {
+    Length const earth_distance =
+        (at_спутник_1_launch->trajectory(*ephemeris, earth_name).
+             EvaluatePosition(final_time) -
+         trajectory.last().degrees_of_freedom().position()).Norm();
+    ss << earth_distance << " ";
+  }
+  state.SetLabel(ss.str());
+}
+
 }  // namespace
 
 void BM_EphemerisSolarSystemMajorBodiesOnly(benchmark::State& state) {
@@ -383,6 +467,12 @@ void FlowEphemerisWithFixedStepSRKN(
   ephemeris.FlowWithFixedStep(t, *instance);
 }
 
+BENCHMARK(BM_EphemerisMultithreadingBenchmark)
+    ->ArgPair(3, 1)
+    ->ArgPair(3, 2)
+    ->ArgPair(3, 3)
+    ->ArgPair(3, 4)
+    ->ArgPair(3, 5);
 BENCHMARK(BM_EphemerisSolarSystemMajorBodiesOnly)->Arg(-3);
 BENCHMARK(BM_EphemerisSolarSystemMinorAndMajorBodies)->Arg(-3);
 BENCHMARK(BM_EphemerisSolarSystemAllBodiesAndOblateness)->Arg(-3);

diff --git a/physics/solar_system_body.hpp b/physics/solar_system_body.hpp
@@ -207,8 +207,8 @@ Length SolarSystem<Frame>::mean_radius(std::string const& name) const {
 
 template<typename Frame>
 not_null<MassiveBody const*> SolarSystem<Frame>::massive_body(
-    Ephemeris<Frame> const & ephemeris,
-    std::string const & name) const {
+    Ephemeris<Frame> const& ephemeris,
+    std::string const& name) const {
   return ephemeris.bodies()[index(name)];
 }
 
@@ -223,8 +223,8 @@ not_null<RotatingBody<Frame> const*> SolarSystem<Frame>::rotating_body(
 
 template<typename Frame>
 ContinuousTrajectory<Frame> const& SolarSystem<Frame>::trajectory(
-    Ephemeris<Frame> const & ephemeris,
-    std::string const & name) const {
+    Ephemeris<Frame> const& ephemeris,
+    std::string const& name) const {
   MassiveBody const* const body = ephemeris.bodies()[index(name)];
   return *ephemeris.trajectory(body);
 }