Merge pull request #1606 from pleroy/1595

Don't lose accuracy when computing the times in the integrators

integrators/embedded_explicit_runge_kutta_nyström_integrator_body.hpp

@@ -211,7 +211,7 @@ Instance::Solve(Instant const& t_final) {
 
       // Runge-Kutta-Nyström iteration; fills |g|.
       for (int i = first_stage; i < stages; ++i) {
-        Instant const t_stage = t.value + c[i] * h;
+        Instant const t_stage = t.value + (t.error + c[i] * h);
         for (int k = 0; k < dimension; ++k) {
           Acceleration Σj_a_ij_g_jk{};
           for (int j = 0; j < i; ++j) {

integrators/embedded_explicit_runge_kutta_nyström_integrator_test.cpp

@@ -316,11 +316,11 @@ TEST_F(EmbeddedExplicitRungeKuttaNyströmIntegratorTest, Singularity) {
   auto const outcome = instance->Solve(t_final);
 
   EXPECT_EQ(termination_condition::VanishingStepSize, outcome.error());
-  EXPECT_EQ(130, solution.size());
+  EXPECT_EQ(132, solution.size());
   EXPECT_THAT(solution.back().time.value - t_initial,
-              AlmostEquals(t_singular - t_initial, 20));
+              AlmostEquals(t_singular - t_initial, 15));
   EXPECT_THAT(solution.back().positions.back().value,
-              AlmostEquals(specific_impulse * initial_mass / mass_flow, 711));
+              AlmostEquals(specific_impulse * initial_mass / mass_flow, 540));
 }
 
 TEST_F(EmbeddedExplicitRungeKuttaNyströmIntegratorTest, Restart) {

integrators/symplectic_runge_kutta_nyström_integrator_body.hpp

@@ -110,7 +110,7 @@ Instance::Solve(Instant const& t_final) {
       for (int k = 0; k < dimension; ++k) {
         q_stage[k] = q[k].value + Δq[k];
       }
-      equation.compute_acceleration(t.value + c[i] * h, q_stage, g);
+      equation.compute_acceleration(t.value + (t.error + c[i] * h), q_stage, g);
       for (int k = 0; k < dimension; ++k) {
         // exp(bᵢ h B)
         Δv[k] += h * b[i] * g[k];

journal/player_test.cpp

@@ -92,7 +92,7 @@ TEST_F(PlayerTest, DISABLED_Debug) {
   // An example of how journaling may be used for debugging.  You must set
   // |path| and fill the |method_in| and |method_out_return| protocol buffers.
   std::string path =
-      R"(\\venezia.mockingbirdnest.com\Namespaces\Public\Public Mockingbird\Principia\Crashes\1441\replication\JOURNAL.20170616-195725)";  // NOLINT
+      R"(\\venezia.mockingbirdnest.com\Namespaces\Public\Public Mockingbird\Principia\Crashes\1595\JOURNAL.20171014-075517)";  // NOLINT
   Player player(path);
   int count = 0;
   while (player.Play()) {
@@ -110,21 +110,20 @@ TEST_F(PlayerTest, DISABLED_Debug) {
   serialization::Method method_in;
   {
     auto* extension = method_in.MutableExtension(
-        serialization::AdvanceTime::extension);
+        serialization::UpdatePrediction::extension);
     auto* in = extension->mutable_in();
-    in->set_plugin(239018496);
-    in->set_t(55066.065156232558);
-    in->set_planetarium_rotation(-220.16314278280507);
+    in->set_plugin(25543983152);
+    in->set_vessel_guid("1ba916b3-ec0f-494a-a5bb-74956b745a58");
   }
   serialization::Method method_out_return;
   {
     auto* extension = method_out_return.MutableExtension(
-        serialization::AdvanceTime::extension);
+        serialization::UpdatePrediction::extension);
   }
   LOG(ERROR) << "Running unpaired method:\n" << method_in.DebugString();
-  CHECK(RunIfAppropriate<AdvanceTime>(method_in,
-                                      method_out_return,
-                                      player));
+  CHECK(RunIfAppropriate<UpdatePrediction>(method_in,
+                                           method_out_return,
+                                           player));
 #endif
 }
 

ksp_plugin_test/interface_test.cpp

@@ -582,7 +582,7 @@ TEST_F(InterfaceTest, DISABLED_DeserializePluginDebug) {
   PushDeserializer* deserializer = nullptr;
   Plugin const* plugin = nullptr;
   std::string const hexadecimal_plugin = ReadFromHexadecimalFile(
-      R"(P:\Public Mockingbird\Principia\Crashes\1422\persistent.proto.hex)");
+      R"(P:\Public Mockingbird\Principia\Crashes\1595\persistent.proto.hex)");
   principia__DeserializePlugin(
           hexadecimal_plugin.c_str(),
           hexadecimal_plugin.size(),

physics/continuous_trajectory_body.hpp

@@ -1,6 +1,8 @@
 
 #pragma once
 
+#include "physics/continuous_trajectory.hpp"
+
 #include <algorithm>
 #include <limits>
 #include <sstream>
@@ -10,7 +12,6 @@
 #include "astronomy/epoch.hpp"
 #include "glog/stl_logging.h"
 #include "numerics/ulp_distance.hpp"
-#include "physics/continuous_trajectory.hpp"
 #include "quantities/si.hpp"
 
 namespace principia {

physics/ephemeris_test.cpp

@@ -497,9 +497,10 @@ TEST_P(EphemerisTest, EarthProbe) {
   }
   // The solution is a line, so the rounding errors dominate.  Different
   // compilers result in different errors and thus different numbers of steps.
-  EXPECT_THAT(probe_positions.size(), AnyOf(Eq(434), Eq(419), Eq(-1)));
+  EXPECT_THAT(probe_positions.size(),
+              AnyOf(Eq(400), Eq(386), Eq(379)));
   EXPECT_THAT(probe_positions.back().coordinates().x,
-              AlmostEquals(1.00 * period * v_probe, 194, 207));
+              AlmostEquals(1.00 * period * v_probe, 194, 235));
   EXPECT_THAT(probe_positions.back().coordinates().y,
               Eq(q_probe));
 
@@ -650,9 +651,9 @@ TEST_P(EphemerisTest, EarthTwoProbes) {
   EXPECT_THAT(probe2_positions.size(), Eq(1001));
 #endif
   EXPECT_THAT(probe1_positions.back().coordinates().x,
-              AlmostEquals(1.00 * period * v_probe1, 40, 63));
+              AlmostEquals(1.00 * period * v_probe1, 40, 70));
   EXPECT_THAT(probe2_positions.back().coordinates().x,
-              AlmostEquals(1.00 * period * v_probe2, 1));
+              AlmostEquals(1.00 * period * v_probe2, 1, 2));
   EXPECT_THAT(probe1_positions.back().coordinates().y,
               Eq(q_probe1));
   EXPECT_THAT(probe2_positions.back().coordinates().y,