From b80b7bf3a0554a9ff9addc77ce6312e955d2113f Mon Sep 17 00:00:00 2001 From: Robin Leroy Date: Thu, 1 Nov 2018 00:54:53 +0100 Subject: [PATCH] better accuracy parameters --- physics/ephemeris.hpp | 5 +- physics/ephemeris_body.hpp | 11 ++-- physics/geopotential.hpp | 61 ++++++++++---------- physics/geopotential_body.hpp | 98 ++++++++++++++++----------------- physics/geopotential_test.cpp | 8 +-- serialization/physics.proto | 2 +- suppress_useless_warnings.props | 2 +- 7 files changed, 94 insertions(+), 93 deletions(-) diff --git a/physics/ephemeris.hpp b/physics/ephemeris.hpp index 7dbfa2a30d..2202deb646 100644 --- a/physics/ephemeris.hpp +++ b/physics/ephemeris.hpp @@ -75,7 +75,7 @@ class Ephemeris { // Implicit for compatibility. AccuracyParameters(Length const& fitting_tolerance); // NOLINT AccuracyParameters(Length const& fitting_tolerance, - serialization::Numerics::Mode geopotential_mode); + double geopotential_tolerance); void WriteToMessage( not_null const @@ -85,8 +85,7 @@ class Ephemeris { private: Length fitting_tolerance_; - serialization::Numerics::Mode geopotential_mode_ = - serialization::Numerics::PRECISE; + double geopotential_tolerance_ = 0; friend class Ephemeris; }; diff --git a/physics/ephemeris_body.hpp b/physics/ephemeris_body.hpp index 1e9f98e23b..d4d050ba43 100644 --- a/physics/ephemeris_body.hpp +++ b/physics/ephemeris_body.hpp @@ -121,16 +121,16 @@ Ephemeris::AccuracyParameters::AccuracyParameters( template Ephemeris::AccuracyParameters::AccuracyParameters( Length const& fitting_tolerance, - serialization::Numerics::Mode const geopotential_mode) + double const geopotential_tolerance) : fitting_tolerance_(fitting_tolerance), - geopotential_mode_(geopotential_mode) {} + geopotential_tolerance_(geopotential_tolerance) {} template void Ephemeris::AccuracyParameters::WriteToMessage( not_null const message) const { fitting_tolerance_.WriteToMessage(message->mutable_fitting_tolerance()); - message->set_geopotential_mode(geopotential_mode_); + message->set_geopotential_tolerance(geopotential_tolerance_); } template @@ -139,7 +139,7 @@ Ephemeris::AccuracyParameters::ReadFromMessage( serialization::Ephemeris::AccuracyParameters const& message) { return AccuracyParameters( Length::ReadFromMessage(message.fitting_tolerance()), - message.geopotential_mode()); + message.geopotential_tolerance()); } template @@ -302,7 +302,8 @@ Ephemeris::Ephemeris( if (body->is_oblate()) { geopotentials_.emplace( geopotentials_.cbegin(), - dynamic_cast_not_null const*>(body.get())); + dynamic_cast_not_null const*>(body.get()), + accuracy_parameters_.geopotential_tolerance_); // Inserting at the beginning of the vectors is O(N). bodies_.insert(bodies_.begin(), std::move(body)); trajectories_.insert(trajectories_.begin(), trajectory); diff --git a/physics/geopotential.hpp b/physics/geopotential.hpp index 186966ee9d..d5e279bf4a 100644 --- a/physics/geopotential.hpp +++ b/physics/geopotential.hpp @@ -25,6 +25,37 @@ using quantities::Length; using quantities::Quotient; using quantities::Square; +// Specification of the damping of a spherical harmonic, acting as a radial +// multiplier on the potential: +// V_damped = σ(|r|) V(r). +struct HarmonicDamping final { + explicit HarmonicDamping() = default; + explicit HarmonicDamping(Length const& inner_threshold); + + // Above this threshold, the contribution to the potential from this + // harmonic is 0, i.e., σ = 0. + Length outer_threshold = Infinity(); + // Below this threshold, the contribution to the potential from this + // harmonic is undamped, σ = 1. + // inner_threshold = outer_threshold / 3. + Length inner_threshold = Infinity(); + // For r in [outer_threshold, inner_threshold], σ is a polynomial with the + // following coefficients in monomial basis. + // The constant term is always 0, and is thus ignored in the evaluation. + Derivatives sigmoid_coefficients; + + // Sets σℜ_over_r and grad_σℜ according to σ as defined by |*this|. + template + void ComputeDampedRadialQuantities( + Length const& r_norm, + Square const& r², + Vector const& r_normalized, + Inverse> const& ℜ_over_r, + Inverse> const& ℜʹ, + Inverse>& σℜ_over_r, + Vector>, Frame>& grad_σℜ) const; +}; + // Representation of the geopotential model of an oblate body. template class Geopotential { @@ -77,36 +108,6 @@ class Geopotential { template struct AllDegrees; - // Specification of the damping of a spherical harmonic, acting as - // a radial multiplier on the potential: - // V_damped = σ(|r|) V(r). - struct HarmonicDamping { - explicit HarmonicDamping() = default; - explicit HarmonicDamping(Length const& inner_threshold); - - // Above this threshold, the contribution to the potential from this - // harmonic is 0, i.e., σ = 0. - Length const outer_threshold = Infinity(); - // Below this threshold, the contribution to the potential from this - // harmonic is undamped, σ = 1. - // inner_threshold = outer_threshold / 3. - Length const inner_threshold = Infinity(); - // For r in [outer_threshold, inner_threshold], σ is a polynomial with the - // following coefficients in monomial basis. - // The constant term is always 0, and is thus ignored in the evaluation. - Derivatives const sigmoid_coefficients; - - // Sets σℜ_over_r and grad_σℜ according to σ as defined by |*this|. - void ComputeDampedRadialQuantities( - Length const& r_norm, - Square const& r², - Vector const& r_normalized, - Inverse> const& ℜ_over_r, - Inverse> const& ℜʹ, - Inverse>& σℜ_over_r, - Vector>, Frame>& grad_σℜ) const; - }; - // If z is a unit vector along the axis of rotation, and r a vector from the // center of |body_| to some point in space, the acceleration computed here // is: diff --git a/physics/geopotential_body.hpp b/physics/geopotential_body.hpp index 13bea82585..0847b6f7e5 100644 --- a/physics/geopotential_body.hpp +++ b/physics/geopotential_body.hpp @@ -38,6 +38,50 @@ using quantities::Sin; using quantities::SIUnit; // The notation in this file follows documentation/Geopotential.pdf. + +inline HarmonicDamping::HarmonicDamping( + Length const& inner_threshold) + : outer_threshold(inner_threshold * 3), + inner_threshold(outer_threshold), + sigmoid_coefficients{0, + 9 / (4 * inner_threshold), + -3 / (2 * Pow<2>(inner_threshold)), + 1 / (4 * Pow<3>(inner_threshold))} {} + +template +void HarmonicDamping::ComputeDampedRadialQuantities( + Length const& r_norm, + Square const& r², + Vector const& r_normalized, + Inverse> const& ℜ_over_r, + Inverse> const& ℜʹ, + Inverse>& σℜ_over_r, + Vector>, Frame>& grad_σℜ) const { + Length const& s1 = outer_threshold; + Length const& s0 = inner_threshold; + if (r_norm <= s0) { + // Below the inner threshold, σ = 1. + σℜ_over_r = ℜ_over_r; + grad_σℜ = ℜʹ * r_normalized; + } else { + auto const& c = sigmoid_coefficients; + Derivative const c1 = std::get<1>(c); + Derivative const c2 = std::get<2>(c); + Derivative const c3 = std::get<3>(c); + auto const r³ = r² * r_norm; + double const c3r³ = c3 * r³; + double const c2r² = c2 * r²; + double const c1r = c1 * r_norm; + double const σ = c3r³ + c2r² + c1r; + double const σʹr = 3 * c3r³ + 2 * c2r² + c1r; + + σℜ_over_r = σ * ℜ_over_r; + // Writing this as σ′ℜ + ℜ′σ rather than ℜ∇σ + σ∇ℜ turns some vector + // operations into scalar ones. + grad_σℜ = (σʹr * ℜ_over_r + ℜʹ * σ) * r_normalized; + } +} + template template struct Geopotential::Precomputations { @@ -243,50 +287,6 @@ auto Geopotential::DegreeNOrderM::Acceleration( } } -template -Geopotential::HarmonicDamping::HarmonicDamping( - Length const& inner_threshold) - : outer_threshold(inner_threshold * 3), - inner_threshold(outer_threshold), - sigmoid_coefficients{0, - 9 / (4 * inner_threshold), - -3 / (2 * Pow<2>(inner_threshold)), - 1 / (4 * Pow<3>(inner_threshold))} {} - -template -void Geopotential::HarmonicDamping::ComputeDampedRadialQuantities( - Length const& r_norm, - Square const& r², - Vector const& r_normalized, - Inverse> const& ℜ_over_r, - Inverse> const& ℜʹ, - Inverse>& σℜ_over_r, - Vector>, Frame>& grad_σℜ) const { - Length const& s1 = outer_threshold; - Length const& s0 = inner_threshold; - if (r_norm <= s0) { - // Below the inner threshold, σ = 1. - σℜ_over_r = ℜ_over_r; - grad_σℜ = ℜʹ * r_normalized; - } else { - auto const& c = sigmoid_coefficients; - Derivative const c1 = std::get<1>(c); - Derivative const c2 = std::get<2>(c); - Derivative const c3 = std::get<3>(c); - auto const r³ = r² * r_norm; - double const c3r³ = c3 * r³; - double const c2r² = c2 * r²; - double const c1r = c1 * r_norm; - double const σ = c3r³ + c2r² + c1r; - double const σʹr = 3 * c3r³ + 2 * c2r² + c1r; - - σℜ_over_r = σ * ℜ_over_r; - // Writing this as σ′ℜ + ℜ′σ rather than ℜ∇σ + σ∇ℜ turns some vector - // operations into scalar ones. - grad_σℜ = (σʹr * ℜ_over_r + ℜʹ * σ) * r_normalized; - } -} - template template auto Geopotential:: @@ -492,11 +492,11 @@ Geopotential::Geopotential(not_null const*> body, tesseral_damping_ = HarmonicDamping(Length{}); Length tesseral_threshold; bool is_tesseral = false; - degree_damping_.push_back(); - degree_damping_.push_back(); + degree_damping_.emplace_back(); + degree_damping_.emplace_back(); for (int n = 2; n <= body_->geopotential_degree(); ++n) { Length degree_n_threshold = degree_damping_[n - 1].inner_threshold; - for (m = 0; m <= n; ++m) { + for (int m = 0; m <= n; ++m) { double const max_abs_Pnm = MaxAbsNormalizedAssociatedLegendreFunction[n][m]; double const Cnm = body->cos()[n][m]; @@ -504,8 +504,8 @@ Geopotential::Geopotential(not_null const*> body, // TODO(egg): write a rootn. Length const r = body->reference_radius() * - std::exp( - (max_abs_Pnm * (1 + n) * Sqrt(Pow<2>(Cnm) + Pow<2>(Snm))) / ε, + std::pow( + (max_abs_Pnm * (n + 1) * Sqrt(Pow<2>(Cnm) + Pow<2>(Snm))) / ε, 1.0 / n); if (m == 0 || is_tesseral) { degree_n_threshold = std::max(r, degree_n_threshold); diff --git a/physics/geopotential_test.cpp b/physics/geopotential_test.cpp index e50cd9a96f..cca9e40f22 100644 --- a/physics/geopotential_test.cpp +++ b/physics/geopotential_test.cpp @@ -97,7 +97,7 @@ TEST_F(GeopotentialTest, J2) { OblateBody(massive_body_parameters_, rotating_body_parameters_, OblateBody::Parameters(/*j2=*/6, 1 * Metre)); - Geopotential const geopotential(&body); + Geopotential const geopotential(&body, /*tolerance=*/0); // The acceleration at a point located on the axis is along the axis. { @@ -169,7 +169,7 @@ TEST_F(GeopotentialTest, C22S22) { rotating_body_parameters_, OblateBody::Parameters::ReadFromMessage( message, 1 * Metre)); - Geopotential const geopotential(&body); + Geopotential const geopotential(&body, /*tolerance=*/0); // The acceleration at a point located on the axis is along the axis for the // (2, 2) harmonics. @@ -224,7 +224,7 @@ TEST_F(GeopotentialTest, J3) { rotating_body_parameters_, OblateBody::Parameters::ReadFromMessage( message, 1 * Metre)); - Geopotential const geopotential(&body); + Geopotential const geopotential(&body, /*tolerance=*/0); // The acceleration at a point located on the axis is along the axis. { @@ -279,7 +279,7 @@ TEST_F(GeopotentialTest, TestVector) { rotating_body_parameters, OblateBody::Parameters::ReadFromMessage( earth_message.geopotential(), earth_reference_radius)); - Geopotential const geopotential(&earth); + Geopotential const geopotential(&earth, /*tolerance=*/0); // This test vector is from Kuga & Carrara, "Fortran- and C-codes for higher // order and degree geopotential computation", diff --git a/serialization/physics.proto b/serialization/physics.proto index 692aad5fb5..87b29a974a 100644 --- a/serialization/physics.proto +++ b/serialization/physics.proto @@ -98,7 +98,7 @@ message Ephemeris { // Added in Ἐρατοσθένης. message AccuracyParameters { required Quantity fitting_tolerance = 1; - required Numerics.Mode geopotential_mode = 2; + required double geopotential_tolerance = 2; } message AdaptiveStepParameters { required AdaptiveStepSizeIntegrator integrator = 1; diff --git a/suppress_useless_warnings.props b/suppress_useless_warnings.props index b0a26dbf53..1f4f15266c 100644 --- a/suppress_useless_warnings.props +++ b/suppress_useless_warnings.props @@ -5,7 +5,7 @@ - 4018;4244;4251;4267;4373;4503;4506;4554;4566;4624;4800;4996;%(DisableSpecificWarnings) + 4018;4244;4251;4267;4373;4503;4506;4554;4566;4624;4723;4800;4996;%(DisableSpecificWarnings)