Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: e94a7f85dd91
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 639ce0a999e6
Choose a head ref
  • 16 commits
  • 14 files changed
  • 1 contributor

Commits on May 31, 2020

  1. KSP resonance test.

    pleroy committed May 31, 2020
    Copy the full SHA
    58d0b2d View commit details
  2. KSP system test.

    pleroy committed May 31, 2020
    Copy the full SHA
    3412d9c View commit details

Commits on Jun 1, 2020

  1. New Logger::Set function.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    3662f4e View commit details
  2. Lunar orbit test.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    2211a56 View commit details
  3. If you test you'll find bugs.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    b475a2c View commit details
  4. Regression test support.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    df32831 View commit details
  5. Copy the full SHA
    9f1623d View commit details
  6. Solar system dynamics test.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    10d9c2a View commit details
  7. Trappist dynamics.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    d2d52e3 View commit details
  8. Molniya orbit.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    ac41bf1 View commit details
  9. Integrators and physics.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    36224ab View commit details
  10. Remove unused dependencies.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    393e44c View commit details
  11. Regression test notebook.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    00a0d49 View commit details
  12. wl file instead of nb.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    fa17fb1 View commit details
  13. Lint.

    pleroy committed Jun 1, 2020
    Copy the full SHA
    e632482 View commit details
  14. Merge pull request #2599 from pleroy/MathematicaLogger

    Use the mathematica::Logger everywhere to output Mathematica traces
    pleroy authored Jun 1, 2020
    Copy the full SHA
    639ce0a View commit details
31 changes: 16 additions & 15 deletions astronomy/ksp_resonance_test.cpp
Original file line number Diff line number Diff line change
@@ -5,7 +5,6 @@
#include <vector>

#include "astronomy/stabilize_ksp.hpp"
#include "base/file.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "integrators/methods.hpp"
@@ -26,9 +25,9 @@ namespace principia {
using base::Error;
using base::make_not_null_unique;
using base::not_null;
using base::OFStream;
using geometry::AngularVelocity;
using geometry::BarycentreCalculator;
using geometry::Displacement;
using geometry::Frame;
using geometry::Inertial;
using geometry::Instant;
@@ -126,37 +125,39 @@ class KSPResonanceTest : public ::testing::Test {
Instant const& t_min,
Instant const& t_max,
std::string const& name) {
// Mathematica tends to be slow when dealing with quantities, so we give
// everything in SI units.
std::vector<double> times;
// Indexed chronologically, then by body.
std::vector<std::vector<Vector<double, KSP>>> barycentric_positions;
mathematica::Logger logger(TEMP_DIR / (name + ".generated.wl"),
/*make_unique=*/false);

for (Instant t = t_min; t <= t_max; t += Δt) {
auto const position = [&ephemeris, t](not_null<MassiveBody const*> body) {
return ephemeris.trajectory(body)->EvaluatePosition(t);
};

times.emplace_back((t - solar_system_.epoch()) / Second);

BarycentreCalculator<Position<KSP>, GravitationalParameter>
jool_system_barycentre;
for (not_null<MassiveBody const*> const body : jool_system_) {
jool_system_barycentre.Add(position(body),
body->gravitational_parameter());
}
barycentric_positions.emplace_back();

// Indexed by body.
std::vector<Displacement<KSP>> barycentric_positions;
for (not_null<MassiveBody const*> const body : jool_system_) {
// TODO(egg): when our dynamic frames support that, it would make sense
// to use a nonrotating dynamic frame centred at the barycentre of the
// Jool system, instead of computing the barycentre and difference
// ourselves.
barycentric_positions.back().emplace_back(
(position(body) - jool_system_barycentre.Get()) / Metre);
barycentric_positions.emplace_back(
position(body) - jool_system_barycentre.Get());
}

logger.Append(name + "t",
t - solar_system_.epoch(),
mathematica::ExpressIn(Second));
logger.Append(name + "q",
barycentric_positions,
mathematica::ExpressIn(Metre));
}
OFStream file(TEMP_DIR / (name + ".generated.wl"));
file << mathematica::Assign(name + "q", barycentric_positions);
file << mathematica::Assign(name + "t", times);
}

// Compute and log the measured periods of the moons.
200 changes: 80 additions & 120 deletions astronomy/ksp_system_test.cpp
Original file line number Diff line number Diff line change
@@ -7,8 +7,8 @@
#include <tuple>
#include <vector>

#include "absl/strings/ascii.h"
#include "astronomy/stabilize_ksp.hpp"
#include "base/file.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "mathematica/mathematica.hpp"
@@ -25,8 +25,8 @@
namespace principia {

using base::not_null;
using base::OFStream;
using geometry::BarycentreCalculator;
using geometry::Displacement;
using geometry::Frame;
using geometry::Inertial;
using geometry::Instant;
@@ -62,6 +62,8 @@ using quantities::si::Milli;
using quantities::si::Minute;
using quantities::si::Second;
using ::testing::Lt;
using ::testing::Matcher;
using ::testing::_;

namespace astronomy {

@@ -145,10 +147,12 @@ class KSPSystemTest : public ::testing::Test, protected KSPSystem {
jool_system_{jool_, laythe_, vall_, tylo_, pol_, bop_},
joolian_moons_{laythe_, vall_, tylo_, pol_, bop_} {}

void FillPositions(Ephemeris<KSP> const& ephemeris,
Instant const& initial_time,
Time const& duration,
std::vector<std::vector<Vector<double, KSP>>>& container) {
void LogPositions(Ephemeris<KSP> const& ephemeris,
Instant const& initial_time,
Time const& duration,
Matcher<Length> const& matcher,
std::string const& name,
mathematica::Logger& logger) {
for (Instant t = initial_time;
t < initial_time + duration;
t += 45 * Minute) {
@@ -162,11 +166,13 @@ class KSPSystemTest : public ::testing::Test, protected KSPSystem {
jool_system_barycentre.Add(position(body),
body->gravitational_parameter());
}
container.emplace_back();
std::vector<Displacement<KSP>> barycentric_positions;
for (not_null<MassiveBody const*> const body : jool_system_) {
container.back().emplace_back(
(position(body) - jool_system_barycentre.Get()) / Metre);
barycentric_positions.emplace_back(
position(body) - jool_system_barycentre.Get());
EXPECT_THAT(barycentric_positions.back().Norm(), matcher);
}
logger.Append(name, barycentric_positions, mathematica::ExpressIn(Metre));
}
}

@@ -198,6 +204,8 @@ class KSPSystemTest : public ::testing::Test, protected KSPSystem {
#if !defined(_DEBUG)
TEST_F(KSPSystemTest, KerbalSystem) {
google::LogToStderr();
mathematica::Logger logger(TEMP_DIR / "ksp_system.generated.wl",
/*make_unique=*/false);

#if 0
auto const a_century_hence = solar_system_.epoch() + 100 * JulianYear;
@@ -209,29 +217,11 @@ TEST_F(KSPSystemTest, KerbalSystem) {
ephemeris_->Prolong(a_century_hence);
LOG(INFO) << "Integration done";

std::map<not_null<MassiveBody const*>, std::vector<double>>
extremal_separations_in_m;
std::map<not_null<MassiveBody const*>, std::vector<double>> times_in_s;
std::map<not_null<MassiveBody const*>, Length> last_separations;
std::map<not_null<MassiveBody const*>, Sign> last_separation_changes;

Instant t = solar_system_.epoch();

// Stock elements.
std::vector<double> bop_eccentricities;
std::vector<double> bop_inclinations_in_degrees;
std::vector<double> bop_nodes_in_degrees;
std::vector<double> bop_arguments_of_periapsis_in_degrees;

// Elements around the barycentre of Jool, Laythe, and Vall.
std::vector<double> bop_jacobi_eccentricities;
std::vector<double> bop_jacobi_nodes_in_degrees;
std::vector<double> bop_jacobi_inclinations_in_degrees;
std::vector<double> bop_jacobi_arguments_of_periapsis_in_degrees;

std::vector<double> tylo_bop_separations_in_m;
std::vector<double> pol_bop_separations_in_m;

for (not_null<MassiveBody const*> const moon : joolian_moons_) {
last_separation_changes.emplace(moon, Sign::Positive());
}
@@ -249,22 +239,27 @@ TEST_F(KSPSystemTest, KerbalSystem) {
auto const jool_position = position(jool_);

for (not_null<MassiveBody const*> const moon : joolian_moons_) {
std::string const moon_name = absl::AsciiStrToLower(moon->name());
Length const separation = (jool_position - position(moon)).Norm();
Sign const separation_change = Sign(separation - last_separations[moon]);
if (separation_change != last_separation_changes.at(moon)) {
extremal_separations_in_m[moon].emplace_back(last_separations[moon] /
Metre);
times_in_s[moon].emplace_back((t - 1 * Hour - solar_system_.epoch()) /
Second);
logger.Append(moon_name + "Separations",
last_separations[moon],
mathematica::ExpressIn(Metre));
logger.Append(moon_name + "Times",
t - 1 * Hour - solar_system_.epoch(),
mathematica::ExpressIn(Second));
}
last_separations[moon] = separation;
last_separation_changes.at(moon) = separation_change;
}

tylo_bop_separations_in_m.emplace_back(
(position(tylo_) - position(bop_)).Norm() / Metre);
pol_bop_separations_in_m.emplace_back(
(position(pol_) - position(bop_)).Norm() / Metre);
logger.Append("tyloBop",
(position(tylo_) - position(bop_)).Norm(),
mathematica::ExpressIn(Metre));
logger.Append("polBop",
(position(pol_) - position(bop_)).Norm(),
mathematica::ExpressIn(Metre));

{
// KSP's osculating elements.
@@ -273,13 +268,16 @@ TEST_F(KSPSystemTest, KerbalSystem) {
MasslessBody(),
degrees_of_freedom(bop_) - degrees_of_freedom(jool_),
t).elements_at_epoch();
bop_eccentricities.emplace_back(*bop_elements.eccentricity);
bop_inclinations_in_degrees.emplace_back(bop_elements.inclination /
Degree);
bop_nodes_in_degrees.emplace_back(
bop_elements.longitude_of_ascending_node / Degree);
bop_arguments_of_periapsis_in_degrees.emplace_back(
*bop_elements.argument_of_periapsis / Degree);
logger.Append("bopEccentricities", *bop_elements.eccentricity);
logger.Append("bopInclinations",
bop_elements.inclination,
mathematica::ExpressIn(Degree));
logger.Append("bopNodes",
bop_elements.longitude_of_ascending_node,
mathematica::ExpressIn(Degree));
logger.Append("bopArguments",
*bop_elements.argument_of_periapsis,
mathematica::ExpressIn(Degree));
}

{
@@ -296,70 +294,32 @@ TEST_F(KSPSystemTest, KerbalSystem) {
*bop_,
degrees_of_freedom(bop_) - innermost_jool_system.Get(),
t).elements_at_epoch();
bop_jacobi_eccentricities.emplace_back(*bop_jacobi_elements.eccentricity);
bop_jacobi_inclinations_in_degrees.emplace_back(
bop_jacobi_elements.inclination / Degree);
bop_jacobi_nodes_in_degrees.emplace_back(
bop_jacobi_elements.longitude_of_ascending_node / Degree);
bop_jacobi_arguments_of_periapsis_in_degrees.emplace_back(
*bop_jacobi_elements.argument_of_periapsis / Degree);
logger.Append("bopJacobiEccentricities",
*bop_jacobi_elements.eccentricity);
logger.Append("bopJacobiInclinations",
bop_jacobi_elements.inclination,
mathematica::ExpressIn(Degree));
logger.Append("bopJacobiNodes",
bop_jacobi_elements.longitude_of_ascending_node,
mathematica::ExpressIn(Degree));
logger.Append("bopJacobiArguments",
*bop_jacobi_elements.argument_of_periapsis,
mathematica::ExpressIn(Degree));
}
}

std::vector<std::vector<Vector<double, KSP>>> barycentric_positions_1_year;
FillPositions(*ephemeris_,
LogPositions(*ephemeris_,
solar_system_.epoch(),
1 * JulianYear,
Lt(3e8 * Metre),
"barycentricPositions1",
logger);
LogPositions(*ephemeris_,
solar_system_.epoch(),
1 * JulianYear,
barycentric_positions_1_year);
std::vector<std::vector<Vector<double, KSP>>> barycentric_positions_2_year;
FillPositions(*ephemeris_,
solar_system_.epoch(),
2 * JulianYear,
barycentric_positions_2_year);

for (auto const& body_positions : barycentric_positions_1_year) {
for (auto const& body_position : body_positions) {
EXPECT_THAT(body_position.Norm(), Lt(3e8));
}
}

OFStream file(TEMP_DIR / "ksp_system.generated.wl");
file << mathematica::Assign("laytheTimes", times_in_s[laythe_]);
file << mathematica::Assign("vallTimes", times_in_s[vall_]);
file << mathematica::Assign("tyloTimes", times_in_s[tylo_]);
file << mathematica::Assign("polTimes", times_in_s[pol_]);
file << mathematica::Assign("bopTimes", times_in_s[bop_]);
file << mathematica::Assign("laytheSeparations",
extremal_separations_in_m[laythe_]);
file << mathematica::Assign("vallSeparations",
extremal_separations_in_m[vall_]);
file << mathematica::Assign("tyloSeparations",
extremal_separations_in_m[tylo_]);
file << mathematica::Assign("polSeparations",
extremal_separations_in_m[pol_]);
file << mathematica::Assign("bopSeparations",
extremal_separations_in_m[bop_]);

file << mathematica::Assign("bopEccentricities", bop_eccentricities);
file << mathematica::Assign("bopInclinations", bop_inclinations_in_degrees);
file << mathematica::Assign("bopNodes", bop_nodes_in_degrees);
file << mathematica::Assign("bopArguments",
bop_arguments_of_periapsis_in_degrees);
file << mathematica::Assign("bopJacobiEccentricities",
bop_jacobi_eccentricities);
file << mathematica::Assign("bopJacobiInclinations",
bop_jacobi_inclinations_in_degrees);
file << mathematica::Assign("bopJacobiNodes", bop_jacobi_nodes_in_degrees);
file << mathematica::Assign("bopJacobiArguments",
bop_jacobi_arguments_of_periapsis_in_degrees);

file << mathematica::Assign("tyloBop", tylo_bop_separations_in_m);
file << mathematica::Assign("polBop", pol_bop_separations_in_m);

file << mathematica::Assign("barycentricPositions1",
barycentric_positions_1_year);
file << mathematica::Assign("barycentricPositions2",
barycentric_positions_2_year);
2 * JulianYear,
_,
"barycentricPositions2",
logger);
}
#endif

@@ -375,8 +335,13 @@ class KSPSystemConvergenceTest
protected KSPSystem {
public:
static void SetUpTestCase() {
file_ = OFStream(SOLUTION_DIR / "mathematica" /
"ksp_system_convergence.generated.wl");
logger_ = new mathematica::Logger(
SOLUTION_DIR / "mathematica" / "ksp_system_convergence.generated.wl",
/*make_unique=*/false);
}

static void TearDownTestCase() {
delete logger_;
}

protected:
@@ -394,14 +359,15 @@ class KSPSystemConvergenceTest
return GetParam().first_step_in_seconds;
}

static OFStream file_;
static mathematica::Logger* logger_;
};

OFStream KSPSystemConvergenceTest::file_;
mathematica::Logger* KSPSystemConvergenceTest::logger_ = nullptr;

// This takes 2 minutes to run.
TEST_P(KSPSystemConvergenceTest, DISABLED_Convergence) {
google::LogToStderr();

Time const integration_duration = 1 * JulianYear;

std::map<std::string, std::vector<DegreesOfFreedom<KSP>>>
@@ -440,12 +406,10 @@ TEST_P(KSPSystemConvergenceTest, DISABLED_Convergence) {
}
}

using MathematicaEntry = std::tuple<Time, Length, std::string, Time>;
using MathematicaEntries = std::vector<MathematicaEntry>;

std::vector<Length> position_errors(iterations() - 1);
std::vector<std::string> worst_body(iterations() - 1);
MathematicaEntries mathematica_entries;
std::string const test_name(
::testing::UnitTest::GetInstance()->current_test_info()->name());
for (int i = 0; i < iterations() - 1; ++i) {
for (auto const& [name, errors] : name_to_errors) {
if (position_errors[i] < errors[i].displacement().Norm()) {
@@ -454,17 +418,13 @@ TEST_P(KSPSystemConvergenceTest, DISABLED_Convergence) {
position_errors[i] = std::max(position_errors[i],
errors[i].displacement().Norm());
}
mathematica_entries.push_back({steps[i + 1],
position_errors[i],
worst_body[i],
durations[i + 1].count() * Second});
logger_->Append(std::string("ppaKSPSystemConvergence") +
test_name[test_name.size() - 1],
std::tuple(steps[i + 1],
position_errors[i],
worst_body[i],
durations[i + 1].count() * Second));
}

std::string const test_name(
::testing::UnitTest::GetInstance()->current_test_info()->name());
file_ << mathematica::Assign(
std::string("ppaKSPSystemConvergence") + test_name[test_name.size() - 1],
mathematica::ToMathematica(mathematica_entries));
}

INSTANTIATE_TEST_CASE_P(
Loading