Skip to content

Commit

Permalink
new constant values now not affecting the tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@darioizzo
darioizzo committed Sep 20, 2024
1 parent e8f153d commit 0e8406b
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 27 deletions.
2 changes: 1 addition & 1 deletion include/kep3/core_astro/constants.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ inline constexpr double CAVENDISH = 6.67430e-11; // Cavendish consta
inline constexpr double MU_SUN = 1.32712440041279419e20; // Sun's gravitational parameter (m^3/s^2 kg) - DE440
inline constexpr double MU_EARTH = 398600435507000.0; // Earth's gravitational parameter (m^3/s^2 kg) - DE440
inline constexpr double MU_MOON = 4902800118000.0; // Moon's gravitational parameter (m^3/s^2 kg) - DE440
inline constexpr double EARTH_VELOCITY = 29784.691831696804; // Earth's velocity. (m/s)
inline constexpr double EARTH_VELOCITY = 29784.69183430911; // Earth's velocity. (m/s)
inline constexpr double EARTH_J2 = 1.08262668E-03;
inline constexpr double EARTH_RADIUS = 6378137; // Earth's radius (m)
inline constexpr double DEG2RAD = (pi / 180.0);
Expand Down
18 changes: 11 additions & 7 deletions test/leg_sims_flanagan_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,20 +246,23 @@ TEST_CASE("compute_mismatch_constraints_test")

TEST_CASE("mismatch_constraints_test2")
{
// We test the correctness of the compute_mismatch_constraints computations against a ground truth (computed with a different program)
// The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep ones.
double AU_OLD = 149597870700.0;
double EV_OLD = 29784.691831696804;
double MU_OLD = 1.32712440018e20;
// We test the correctness of the compute_mismatch_constraints computations against a ground truth (computed with a
// different program)
std::array<std::array<double, 3>, 2> rvs{
{{1 * kep3::AU, 0.1 * kep3::AU, -0.1 * kep3::AU},
{0.2 * kep3::EARTH_VELOCITY, 1 * kep3::EARTH_VELOCITY, -0.2 * kep3::EARTH_VELOCITY}}};
{{1 * AU_OLD, 0.1 * AU_OLD, -0.1 * AU_OLD}, {0.2 * EV_OLD, 1 * EV_OLD, -0.2 * EV_OLD}}};

std::array<std::array<double, 3>, 2> rvf{
{{1.2 * kep3::AU, -0.1 * kep3::AU, 0.1 * kep3::AU},
{-0.2 * kep3::EARTH_VELOCITY, 1.023 * kep3::EARTH_VELOCITY, -0.44 * kep3::EARTH_VELOCITY}}};
{{1.2 * AU_OLD, -0.1 * AU_OLD, 0.1 * AU_OLD}, {-0.2 * EV_OLD, 1.023 * EV_OLD, -0.44 * EV_OLD}}};

double ms = 1500.;
double mf = 1300.;
std::vector<double> throttles
= {0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24};
kep3::leg::sims_flanagan sf(rvs, ms, throttles, rvf, mf, 324.0 * kep3::DAY2SEC, 0.12, 100, kep3::MU_SUN, 0.6);
kep3::leg::sims_flanagan sf(rvs, ms, throttles, rvf, mf, 324.0 * kep3::DAY2SEC, 0.12, 100, MU_OLD, 0.6);
auto retval = sf.compute_mismatch_constraints();
std::vector<double> ground_truth
= {-1.9701274809621304e+11, 4.6965044246848071e+11, -1.5007523306033661e+11, -2.9975151466948650e+04,
Expand Down Expand Up @@ -292,7 +295,8 @@ TEST_CASE("grad_test")
auto grad_a = udp.gradient(x);
auto xgrad = xt::adapt(grad, {1u + 7u + 5u, 17u});
auto xgrad_a = xt::adapt(grad_a, {1u + 7u + 5u, 17u});
REQUIRE(xt::linalg::norm(xgrad - xgrad_a) < 1e-8); // With the high fidelity gradient this is still the best we can achieve
REQUIRE(xt::linalg::norm(xgrad - xgrad_a)
< 1e-8); // With the high fidelity gradient this is still the best we can achieve
}

TEST_CASE("serialization_test")
Expand Down
30 changes: 20 additions & 10 deletions test/stark_problem_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,29 +34,39 @@ TEST_CASE("construct")

TEST_CASE("propagate")
{
auto s = kep3::stark_problem(kep3::MU_SUN, 3000 * kep3::G0, 1e-16);
// The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep ones.
double AU_OLD = 149597870700.0;
double EV_OLD = 29784.691831696804;
double MU_OLD = 1.32712440018e20;

auto s = kep3::stark_problem(MU_OLD, 3000 * kep3::G0, 1e-16);
// Mass is zero
REQUIRE_THROWS_AS(s.propagate({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
REQUIRE_THROWS_AS(s.propagate({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
// Radius is zero
REQUIRE_THROWS_AS(s.propagate({0., 0., 0., 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
REQUIRE_THROWS_AS(s.propagate({0., 0., 0., 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
// Thrust is zero
REQUIRE_NOTHROW(s.propagate({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC));
REQUIRE_NOTHROW(s.propagate({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC));

auto res = s.propagate({kep3::AU, 1000, 1000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
auto res = s.propagate({AU_OLD, 1000, 1000, 0.123, EV_OLD, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
std::array<double, 7> const ground_truth{86447318382.73862, -109196478678.49327, 9774837.693084948, 25535.89383875655, 18933.302734531168, -47.00078590023982, 954.2200884785473};
REQUIRE(L_infinity_norm_rel(res, ground_truth) <=1e-13);
}
TEST_CASE("propagate_var")
{
auto s = kep3::stark_problem(kep3::MU_SUN, 3000 * kep3::G0, 1e-16);
// The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep ones.
double AU_OLD = 149597870700.0;
double EV_OLD = 29784.691831696804;
double MU_OLD = 1.32712440018e20;

auto s = kep3::stark_problem(MU_OLD, 3000 * kep3::G0, 1e-16);
// Mass is zero
REQUIRE_THROWS_AS(s.propagate_var({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
REQUIRE_THROWS_AS(s.propagate_var({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
// Radius is zero
REQUIRE_THROWS_AS(s.propagate_var({0., 0., 0., 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
REQUIRE_THROWS_AS(s.propagate_var({0., 0., 0., 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
// Thrust is zero
REQUIRE_THROWS_AS(s.propagate_var({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC), std::domain_error);
REQUIRE_THROWS_AS(s.propagate_var({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC), std::domain_error);

auto res = s.propagate_var({kep3::AU, 1000, 1000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
auto res = s.propagate_var({AU_OLD, 1000, 1000, 0.123, EV_OLD, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
std::array<double, 7> const ground_truth{86447318382.73862, -109196478678.49327, 9774837.693084948, 25535.89383875655, 18933.302734531168, -47.00078590023982, 954.2200884785473};
REQUIRE(L_infinity_norm_rel(std::get<0>(res), ground_truth) <=1e-13);
}
Expand Down
25 changes: 16 additions & 9 deletions test/ta_stark_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@ using heyoka::taylor_adaptive;
using heyoka::taylor_outcome;

using kep3::ta::get_ta_stark;
using kep3::ta::get_ta_stark_var;
using kep3::ta::get_ta_stark_cache_dim;
using kep3::ta::get_ta_stark_var;
using kep3::ta::get_ta_stark_var_cache_dim;

using kep3_tests::L_infinity_norm_rel;
Expand All @@ -37,9 +37,9 @@ TEST_CASE("caches")
auto ta_cached = get_ta_stark(1e-16);
REQUIRE(get_ta_stark_cache_dim() == 1u);
ta_cached = get_ta_stark(1e-16);
REQUIRE(get_ta_stark_cache_dim() ==1u);
REQUIRE(get_ta_stark_cache_dim() == 1u);
ta_cached = get_ta_stark(1e-8);
REQUIRE(get_ta_stark_cache_dim() ==2u);
REQUIRE(get_ta_stark_cache_dim() == 2u);

// The variational integrator.
REQUIRE(get_ta_stark_var_cache_dim() == 0u);
Expand All @@ -48,7 +48,7 @@ TEST_CASE("caches")
ta_cached = get_ta_stark_var(1e-16);
REQUIRE(get_ta_stark_var_cache_dim() == 1u);
ta_cached = get_ta_stark_var(1e-8);
REQUIRE(get_ta_stark_var_cache_dim() == 2u);
REQUIRE(get_ta_stark_var_cache_dim() == 2u);
}

TEST_CASE("dynamics")
Expand All @@ -72,12 +72,15 @@ TEST_CASE("dynamics")
}
{
// We test a generic case.
// The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep
// ones.
double MU_OLD = 1.32712440018e20;
taylor_adaptive<double> ta(ta_cached); // making a copy as to be able to modify the object.
ta.set_time(0.);
std::vector<double> ic{164557657760.1, 179517444829.19998, 47871318621.12, 32763.159550000004,
23827.7524, 9531.10096, 1200.};
std::copy(ic.begin(), ic.end(), ta.get_state_data());
std::vector<double> pars{kep3::MU_SUN, 3000.*kep3::G0, 0.05, 0., 0.032};
std::vector<double> pars{MU_OLD, 3000. * kep3::G0, 0.05, 0., 0.032};
std::copy(pars.begin(), pars.end(), ta.get_pars_data());
auto out = ta.propagate_until(3888000.0);
std::vector<double> const ground_truth
Expand Down Expand Up @@ -107,24 +110,28 @@ TEST_CASE("variational_dynamics")
std::copy(pars.begin(), pars.end(), ta.get_pars_data());
auto out = ta.propagate_until(2. * kep3::pi);
REQUIRE(std::get<0>(out) == taylor_outcome::time_limit);
REQUIRE(L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7), ic) <= 1e-13);
REQUIRE(L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7), ic)
<= 1e-13);
}
{
// We test a generic case.
// The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep
// ones.
double MU_OLD = 1.32712440018e20;
taylor_adaptive<double> ta(ta_cached); // making a copy as to be able to modify the object.
ta.set_time(0.);
std::vector<double> ic{164557657760.1, 179517444829.19998, 47871318621.12, 32763.159550000004,
23827.7524, 9531.10096, 1200.};
std::copy(ic.begin(), ic.end(), ta.get_state_data());
std::vector<double> pars{kep3::MU_SUN, 3000.*kep3::G0, 0.05, 0., 0.032};
std::vector<double> pars{MU_OLD, 3000. * kep3::G0, 0.05, 0., 0.032};
std::copy(pars.begin(), pars.end(), ta.get_pars_data());
auto out = ta.propagate_until(3888000.0);
std::vector<double> const ground_truth
= {284296823432.0578, 263961690798.0665, 82814214381.94377, 29341.50292902231,
20219.294034700008, 8592.028822618351, 1192.1548315009777};
REQUIRE(std::get<0>(out) == taylor_outcome::time_limit);
REQUIRE(kep3_tests::L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7), ground_truth)
REQUIRE(kep3_tests::L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7),
ground_truth)
<= 1e-13);
}
}

0 comments on commit 0e8406b

Please sign in to comment.