Skip to content

Commit

Permalink
Merge pull request #2574 from eggrobin/angular-momentum-accounting
Browse files Browse the repository at this point in the history
A new approach to conserving angular momentum
eggrobin authored May 16, 2020
2 parents 1d9f60a + b175f10 commit 828547f
Showing 7 changed files with 172 additions and 495 deletions.
6 changes: 5 additions & 1 deletion geometry/rotation_body.hpp
Original file line number Diff line number Diff line change
@@ -198,7 +198,11 @@ Bivector<double, FromFrame> Rotation<FromFrame, ToFrame>::RotationAxis() const {

template<typename FromFrame, typename ToFrame>
Angle Rotation<FromFrame, ToFrame>::RotationAngle() const {
return 2 * ArcTan(quaternion_.imaginary_part().Norm(),
// NOTE(egg): We intentionally use the single-parameter arctangent, as we want
// a result in [-π, π], thus an arctangent in [-π/2, π/2].
// The two-parameter arctangent, with a positive numerator, would lie in
// [0, π], and thus the result in [0, 2π].
return 2 * ArcTan(quaternion_.imaginary_part().Norm() /
quaternion_.real_part());
}

10 changes: 3 additions & 7 deletions ksp_plugin/interface_vessel.cpp
Original file line number Diff line number Diff line change
@@ -192,14 +192,10 @@ XYZ __cdecl principia__VesselVelocity(Plugin const* const plugin,
}

void __cdecl principia__SetAngularMomentumConservation(
bool const correct_orientation,
bool const correct_angular_velocity,
bool const thresholding) {
bool const conserve_angular_momentum) {
journal::Method<journal::SetAngularMomentumConservation> m(
{correct_orientation, correct_angular_velocity, thresholding});
ksp_plugin::PileUp::correct_orientation = correct_orientation;
ksp_plugin::PileUp::correct_angular_velocity = correct_angular_velocity;
ksp_plugin::PileUp::thresholding = thresholding;
{conserve_angular_momentum});
ksp_plugin::PileUp::conserve_angular_momentum = conserve_angular_momentum;
return m.Return();
}

444 changes: 154 additions & 290 deletions ksp_plugin/pile_up.cpp

Large diffs are not rendered by default.

23 changes: 1 addition & 22 deletions ksp_plugin/pile_up.hpp
Original file line number Diff line number Diff line change
@@ -16,7 +16,6 @@
#include "geometry/grassmann.hpp"
#include "geometry/named_quantities.hpp"
#include "integrators/integrators.hpp"
#include "numerics/pid.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/ephemeris.hpp"
#include "physics/euler_solver.hpp"
@@ -47,7 +46,6 @@ using geometry::NonRotating;
using geometry::RigidTransformation;
using geometry::Vector;
using integrators::Integrator;
using numerics::PID;
using physics::DiscreteTrajectory;
using physics::DegreesOfFreedom;
using physics::Ephemeris;
@@ -104,9 +102,7 @@ class PileUp {
virtual ~PileUp();

std::string trace;
static bool correct_orientation;
static bool correct_angular_velocity;
static bool thresholding;
static bool conserve_angular_momentum;

// This class is moveable.
PileUp(PileUp&& pile_up) = default;
@@ -190,16 +186,6 @@ class PileUp {
template<AppendToPartTrajectory append_to_part_trajectory>
void AppendToPart(DiscreteTrajectory<Barycentric>::Iterator it) const;

// Computes a rigid motion that adjusts the pile-up to make L_apparent match
// L_actual.
static RigidMotion<ApparentPileUp, NonRotatingPileUp>
ComputeAngularMomentumCorrection(
Time const& Δt,
Bivector<AngularMomentum, ApparentPileUp> const& L_apparent,
Bivector<AngularMomentum, NonRotatingPileUp> const& L_actual,
InertiaTensor<ApparentPileUp> const& inertia_tensor,
std::string& trace);

// Wrapped in a |unique_ptr| to be moveable.
not_null<std::unique_ptr<absl::Mutex>> lock_;

@@ -255,13 +241,6 @@ class PileUp {
// Called in the destructor.
std::function<void()> deletion_callback_;

// A PID used to smoothen the value of the apparent angular momentum obtained
// from KSP.
PID<Bivector<AngularMomentum, ApparentPileUp>,
Bivector<AngularMomentum, ApparentPileUp>,
/*horizon=*/25,
/*finite_difference_order=*/5> apparent_angular_momentum_controller_;

friend class TestablePileUp;
};

19 changes: 5 additions & 14 deletions ksp_plugin_adapter/main_window.cs
Original file line number Diff line number Diff line change
@@ -239,17 +239,10 @@ protected override void RenderWindow(int window_id) {

private void RenderKSPFeatures() {
if (show_2519_debugging_ui) {
correct_orientation = UnityEngine.GUILayout.Toggle(
correct_orientation,
"Correct orientation");
correct_angular_velocity = UnityEngine.GUILayout.Toggle(
correct_angular_velocity,
"Correct angular velocity");
thresholding = UnityEngine.GUILayout.Toggle(
thresholding,
"Fubini orientation correction");
Interface.SetAngularMomentumConservation(
correct_orientation, correct_angular_velocity, thresholding);
conserve_angular_momentum = UnityEngine.GUILayout.Toggle(
conserve_angular_momentum,
"Conserve angular momentum");
Interface.SetAngularMomentumConservation(conserve_angular_momentum);
string trace = null;
if (FlightGlobals.ActiveVessel != null &&
plugin.HasVessel(FlightGlobals.ActiveVessel.id.ToString())) {
@@ -515,9 +508,7 @@ private void RenderToggleableSection(string name,
// These flags exist to facilitate investigation of #2519.
// They must not be serialized: their non-default values can lead to absurd
// behaviour.
private static bool correct_orientation = true;
private static bool correct_angular_velocity = true;
private static bool thresholding = false;
private static bool conserve_angular_momentum = true;
private static readonly bool show_2519_debugging_ui = true;

private static readonly double[] prediction_length_tolerances_ =
161 changes: 3 additions & 158 deletions ksp_plugin_test/pile_up_test.cpp
Original file line number Diff line number Diff line change
@@ -89,7 +89,6 @@ class TestablePileUp : public PileUp {
using PileUp::DeformPileUpIfNeeded;
using PileUp::AdvanceTime;
using PileUp::NudgeParts;
using PileUp::ComputeAngularMomentumCorrection;

Mass const& mass() const {
return mass_;
@@ -253,49 +252,6 @@ class PileUpTest : public testing::Test {
EXPECT_THAT(pile_up.apparent_part_rigid_motion(), IsEmpty());
}

void CheckAngularMomentumCorrection(
InertiaTensor<Vessel> const& vessel_inertia_tensor,
AngularVelocity<Apparent> const& vessel_angular_velocity,
Bivector<AngularMomentum, NonRotatingPileUp> const& L_actual,
Matcher<Bivector<AngularMomentum, ApparentPileUp>> const&
L_apparent_matcher,
Matcher<Bivector<AngularMomentum, CorrectedPileUp>> const&
L_corrected_matcher) {
// Prepare the motion and inertia of a vessel that yields the L_apparent
// that we want to exercise.
RigidMotion<Vessel, Apparent> const apparent_vessel_rigid_motion(
RigidTransformation<Vessel, Apparent>::Identity(),
vessel_angular_velocity,
Velocity<Apparent>());
MechanicalSystem<Apparent, ApparentPileUp> apparent_system;
apparent_system.AddRigidBody(apparent_vessel_rigid_motion,
/*mass=*/1 * Kilogram,
vessel_inertia_tensor);

Bivector<AngularMomentum, ApparentPileUp> const L_apparent =
apparent_system.AngularMomentum();
EXPECT_THAT(L_apparent, L_apparent_matcher);
InertiaTensor<ApparentPileUp> const inertia_tensor =
apparent_system.InertiaTensor();

std::string trace;
RigidMotion<ApparentPileUp, NonRotatingPileUp> const correction =
TestablePileUp::ComputeAngularMomentumCorrection(
/*Δt=*/0.02 * Second, L_apparent, L_actual, inertia_tensor, trace);
LOG(ERROR) << trace;

MechanicalSystem<NonRotatingPileUp, CorrectedPileUp> corrected_system;
RigidMotion<Vessel, NonRotatingPileUp> const corrected_vessel_rigid_motion =
correction * apparent_system.LinearMotion().Inverse() *
apparent_vessel_rigid_motion;
corrected_system.AddRigidBody(corrected_vessel_rigid_motion,
/*mass=*/1 * Kilogram,
vessel_inertia_tensor);
Bivector<AngularMomentum, CorrectedPileUp> const L_corrected =
corrected_system.AngularMomentum();
EXPECT_THAT(L_corrected, L_corrected_matcher);
}

MockFunction<void()> deletion_callback_;

PartId const part_id1_ = 111;
@@ -322,120 +278,7 @@ class PileUpTest : public testing::Test {
Part p2_;
};

TEST_F(PileUpTest, AngularMomentum) {
{
InertiaTensor<Vessel> const vessel_inertia_tensor(
si::Unit<MomentOfInertia> * R3x3Matrix<double>({1, 0, 0},
{0, 1, 0},
{0, 0, 2}));
AngularVelocity<Apparent> const vessel_angular_velocity(
{0.2 * Radian / Second,
0.1 * Radian / Second,
0.95 * Radian / Second});
Bivector<AngularMomentum, NonRotatingPileUp> const L_actual(
{0 * si::Unit<AngularMomentum>,
0 * si::Unit<AngularMomentum>,
2 * si::Unit<AngularMomentum>});
CheckAngularMomentumCorrection(
vessel_inertia_tensor,
vessel_angular_velocity,
L_actual,
/*L_apparent_matcher=*/
AlmostEquals(Bivector<AngularMomentum, ApparentPileUp>(
{0.2 * si::Unit<AngularMomentum>,
0.1 * si::Unit<AngularMomentum>,
1.9 * si::Unit<AngularMomentum>}),
0),
/*L_corrected_matcher=*/
Componentwise(AlmostEquals(L_actual.coordinates().x, 0),
VanishesBefore(L_actual.Norm(), 0),
AlmostEquals(L_actual.coordinates().z, 0)));
}
{
InertiaTensor<Vessel> const vessel_inertia_tensor(
si::Unit<MomentOfInertia> * R3x3Matrix<double>({1, 0, 0},
{0, 2, 0},
{0, 0, 4}));
AngularVelocity<Apparent> const vessel_angular_velocity(
{0.2 * Radian / Second,
0.05 * Radian / Second,
0.475 * Radian / Second});
Bivector<AngularMomentum, NonRotatingPileUp> const L_actual(
{0 * si::Unit<AngularMomentum>,
0 * si::Unit<AngularMomentum>,
2 * si::Unit<AngularMomentum>});
CheckAngularMomentumCorrection(
vessel_inertia_tensor,
vessel_angular_velocity,
L_actual,
/*L_apparent_matcher=*/
AlmostEquals(Bivector<AngularMomentum, ApparentPileUp>(
{0.2 * si::Unit<AngularMomentum>,
0.1 * si::Unit<AngularMomentum>,
1.9 * si::Unit<AngularMomentum>}),
2),
/*L_corrected_matcher=*/
Componentwise(VanishesBefore(L_actual.Norm(), 0),
VanishesBefore(L_actual.Norm(), 0),
AlmostEquals(L_actual.coordinates().z, 0)));
}
{
InertiaTensor<Vessel> const vessel_inertia_tensor(
si::Unit<MomentOfInertia> * R3x3Matrix<double>({1, 0, 0},
{0, 2, 0},
{0, 0, 4}));
AngularVelocity<Apparent> const vessel_angular_velocity(
{0.02 * Radian / Second,
0.005 * Radian / Second,
0.4975 * Radian / Second});
Bivector<AngularMomentum, NonRotatingPileUp> const L_actual(
{0 * si::Unit<AngularMomentum>,
0 * si::Unit<AngularMomentum>,
2 * si::Unit<AngularMomentum>});
CheckAngularMomentumCorrection(
vessel_inertia_tensor,
vessel_angular_velocity,
L_actual,
/*L_apparent_matcher=*/
AlmostEquals(Bivector<AngularMomentum, ApparentPileUp>(
{0.02 * si::Unit<AngularMomentum>,
0.01 * si::Unit<AngularMomentum>,
1.99 * si::Unit<AngularMomentum>}),
1),
/*L_corrected_matcher=*/
Componentwise(VanishesBefore(L_actual.Norm(), 0),
VanishesBefore(L_actual.Norm(), 0),
AlmostEquals(L_actual.coordinates().z, 1)));
}
{
InertiaTensor<Vessel> const vessel_inertia_tensor(
si::Unit<MomentOfInertia> * R3x3Matrix<double>({1, 0, 0},
{0, 2, 0},
{0, 0, 4}));
AngularVelocity<Apparent> const vessel_angular_velocity(
{3.1 * Radian / Second,
0.45 * Radian / Second,
1.275 * Radian / Second});
Bivector<AngularMomentum, NonRotatingPileUp> const L_actual(
{3 * si::Unit<AngularMomentum>,
1 * si::Unit<AngularMomentum>,
5 * si::Unit<AngularMomentum>});
CheckAngularMomentumCorrection(
vessel_inertia_tensor,
vessel_angular_velocity,
L_actual,
/*L_apparent_matcher=*/
AlmostEquals(Bivector<AngularMomentum, ApparentPileUp>(
{3.1 * si::Unit<AngularMomentum>,
0.9 * si::Unit<AngularMomentum>,
5.1 * si::Unit<AngularMomentum>}),
1),
/*L_corrected_matcher=*/
Componentwise(AlmostEquals(L_actual.coordinates().x, 2),
AlmostEquals(L_actual.coordinates().y, 2),
AlmostEquals(L_actual.coordinates().z, 2)));
}
}
#if 0

// Exercises the entire lifecycle of a |PileUp| that is subject to an intrinsic
// force.
@@ -733,6 +576,8 @@ TEST_F(PileUpTest, LifecycleWithoutIntrinsicForce) {
46)));
}

#endif

TEST_F(PileUpTest, MidStepIntrinsicForce) {
// An empty ephemeris; the parameters don't matter, since there are no bodies
// to integrate.
4 changes: 1 addition & 3 deletions serialization/journal.proto
Original file line number Diff line number Diff line change
@@ -219,9 +219,7 @@ message SetAngularMomentumConservation {
optional SetAngularMomentumConservation extension = 5999;
}
message In {
required bool correct_orientation = 1;
required bool correct_angular_velocity = 2;
required bool thresholding = 3;
required bool conserve_angular_momentum = 1;
}
optional In in = 1;
}

0 comments on commit 828547f

Please sign in to comment.