Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Crash when warping at high speed with short history #2165

Closed
pleroy opened this issue May 11, 2019 · 0 comments · Fixed by #2168
Closed

Crash when warping at high speed with short history #2165

pleroy opened this issue May 11, 2019 · 0 comments · Fixed by #2168
Labels
Milestone

Comments

@pleroy
Copy link
Member

pleroy commented May 11, 2019

F0511 20:11:52.193073 4916 continuous_trajectory_body.hpp:173] Check failed: t_min_locked() <= time (+4.28789530227277544e+05 s vs. +4.27813530227277544e+05 s)
Stack:

    @   000007FED18D915B  	principia::physics::internal_ephemeris::Ephemeris<principia::geometry::internal_frame::Frame<enum principia::serialization::Frame_PluginTag,7,1> >::ComputeGravitationalAccelerationByMassiveBodyOnMasslessBodies<0> [0x000007FED18D915A+90] (c:\users\phl\projects\github\principia\principia\physics\ephemeris_body.hpp:1016)
    @   000007FED196E670  	principia::physics::internal_ephemeris::Ephemeris<principia::geometry::internal_frame::Frame<enum principia::serialization::Frame_PluginTag,7,1> >::ComputeMasslessBodiesGravitationalAccelerations [0x000007FED196E66F+447] (c:\users\phl\projects\github\principia\principia\physics\ephemeris_body.hpp:1121)
    @   000007FED190F409  	<lambda_238bcc44e5612def8e4e635268f83d36>::operator() [0x000007FED190F408+72] (c:\users\phl\projects\github\principia\principia\physics\ephemeris_body.hpp:471)
    @   000007FED19DE4BC  	std::_Func_impl_no_alloc<<lambda_238bcc44e5612def8e4e635268f83d36>,principia::base::Status,principia::geometry::internal_point::Point<principia::quantities::internal_quantities::Quantity<principia::quantities::internal_dimensions::Dimensions<0,0,1,0,0,0,0 [0x000007FED19DE4BB+27] (c:\program files (x86)\microsoft visual studio\preview\enterprise\vc\tools\msvc\14.16.27023\include\functional:16707566)
    @   000007FED22CFF8E  	principia::integrators::internal_embedded_explicit_runge_kutta_nyström_integrator::EmbeddedExplicitRungeKuttaNyströmIntegrator<principia::integrators::methods::Dormand???????Prince1986RKN434FM,principia::geometry::internal_point::Point<principia::geometry [0x000007FED22CFF8D+1933] (c:\users\phl\projects\github\principia\principia\integrators\embedded_explicit_runge_kutta_nyström_integrator_body.hpp:190)
    @   000007FED18E07E4  	principia::physics::internal_ephemeris::Ephemeris<principia::geometry::internal_frame::Frame<enum principia::serialization::Frame_PluginTag,7,1> >::FlowODEWithAdaptiveStep<principia::integrators::internal_ordinary_differential_equations::SpecialSecondOrde [0x000007FED18E07E3+2163] (c:\users\phl\projects\github\principia\principia\physics\ephemeris_body.hpp:1199)
    @   000007FED1979262  	principia::physics::internal_ephemeris::Ephemeris<principia::geometry::internal_frame::Frame<enum principia::serialization::Frame_PluginTag,7,1> >::FlowWithAdaptiveStep [0x000007FED1979261+161] (c:\users\phl\projects\github\principia\principia\physics\ephemeris_body.hpp:481)
    @   000007FED239AB12  	principia::ksp_plugin::internal_vessel::Vessel::FlowPrognostication [0x000007FED239AB11+321] (c:\users\phl\projects\github\principia\principia\ksp_plugin\vessel.cpp:504)
    @   000007FED239CDBC  	principia::ksp_plugin::internal_vessel::Vessel::RepeatedlyFlowPrognostication [0x000007FED239CDBB+843] (c:\users\phl\projects\github\principia\principia\ksp_plugin\vessel.cpp:477)
    @   000007FED2391BCC  	std::_LaunchPad<std::unique_ptr<std::tuple<std::_Binder<std::_Unforced,void (__cdecl principia::base::ThreadPool<principia::base::Status>::*)(void),principia::base::ThreadPool<principia::base::Status> *> >,std::default_delete<std::tuple<std::_Binder<std:: [0x000007FED2391BCB+107] (c:\program files (x86)\microsoft visual studio\preview\enterprise\vc\tools\msvc\14.16.27023\include\thr\xthread:230)
    @   000007FED2327129  	std::_Pad::_Call_func [0x000007FED2327128+8] (c:\program files (x86)\microsoft visual studio\preview\enterprise\vc\tools\msvc\14.16.27023\include\thr\xthread:208)
    @   000007FED954CD70  	o__realloc_base [0x000007FED954CD6F+95]
    @   00000000771C555D  	BaseThreadInitThunk [0x00000000771C555C+12]
    @   000000007742385D  	RtlUserThreadStart [0x000000007742385C+28]

It seems that the parameters passed to the prognosticator have a first_time which is before the t_min protected by the guard.

Decoded stack trace:

void Vessel::RepeatedlyFlowPrognostication() {
std::optional<PrognosticatorParameters> previous_prognosticator_parameters;
for (;;) {
// No point in going faster than 50 Hz.
std::chrono::steady_clock::time_point const wakeup_time =
std::chrono::steady_clock::now() + std::chrono::milliseconds(20);
// The thread is only started after the parameters have been set, so we
// should always find parameters here.
std::optional<PrognosticatorParameters> prognosticator_parameters;
{
absl::ReaderMutexLock l(&prognosticator_lock_);
CHECK(prognosticator_parameters_);
prognosticator_parameters = prognosticator_parameters_;
}
if (prognosticator_parameters->shutdown) {
break;
}
// Do not reflow if the parameters have not changed: the same causes would
// produce the same effects.
if (!previous_prognosticator_parameters ||
*previous_prognosticator_parameters != prognosticator_parameters) {
std::unique_ptr<DiscreteTrajectory<Barycentric>> prognostication;
Status const status = FlowPrognostication(*prognosticator_parameters,
std::unique_ptr<DiscreteTrajectory<Barycentric>>& prognostication) {
// This function may be run in a separate thread, and bad things will happen
// if |EventuallyForgetBefore| runs in parallel with it: the ephemeris'
// |t_min| may end up being after the time that |FlowWithAdaptiveStep| tries
// to integrate, causing various check failures.
Ephemeris<Barycentric>::Guard g(ephemeris_);
prognostication = std::make_unique<DiscreteTrajectory<Barycentric>>();
prognostication->Append(
prognosticator_parameters.first_time,
prognosticator_parameters.first_degrees_of_freedom);
Status status;
status = ephemeris_->FlowWithAdaptiveStep(
bool const last_point_only) {
auto compute_acceleration = [this, &intrinsic_acceleration](
Instant const& t,
std::vector<Position<Frame>> const& positions,
std::vector<Vector<Acceleration, Frame>>& accelerations) {
Error const error =
ComputeMasslessBodiesGravitationalAccelerations(t,
positions,
accelerations);
if (intrinsic_acceleration != nullptr) {
accelerations[0] += intrinsic_acceleration(t);
}
return error == Error::OK ? Status::OK :
CollisionDetected();
};
return FlowODEWithAdaptiveStep<NewtonianMotionEquation>(
if (trajectory_last_time == t) {
return Status::OK;
}
std::vector<not_null<DiscreteTrajectory<Frame>*>> const trajectories =
{trajectory};
// The |min| is here to prevent us from spending too much time computing the
// ephemeris. The |max| is here to ensure that we always try to integrate
// forward. We use |last_state_.time.value| because this is always finite,
// contrary to |t_max()|, which is -∞ when |empty()|.
Instant const t_final =
std::min(std::max(instance_time() +
max_ephemeris_steps * fixed_step_parameters_.step(),
trajectory_last_time + fixed_step_parameters_.step()),
t);
Prolong(t_final);
IntegrationProblem<ODE> problem;
problem.equation.compute_acceleration = std::move(compute_acceleration);
auto const trajectory_last = trajectory->last();
auto const last_degrees_of_freedom = trajectory_last.degrees_of_freedom();
problem.initial_state = {{last_degrees_of_freedom.position()},
{last_degrees_of_freedom.velocity()},
trajectory_last.time()};
typename AdaptiveStepSizeIntegrator<ODE>::Parameters const
integrator_parameters(
/*first_time_step=*/t_final - problem.initial_state.time.value,
/*safety_factor=*/0.9,
parameters.max_steps_,
/*last_step_is_exact=*/true);
CHECK_GT(integrator_parameters.first_time_step, 0 * Second)
<< "Flow back to the future: " << t_final
<< " <= " << problem.initial_state.time.value;
auto const tolerance_to_error_ratio =
std::bind(&Ephemeris<Frame>::ToleranceToErrorRatio,
std::cref(parameters.length_integration_tolerance_),
std::cref(parameters.speed_integration_tolerance_),
_1, _2);
typename AdaptiveStepSizeIntegrator<ODE>::AppendState append_state;
typename ODE::SystemState last_state;
if (last_point_only) {
append_state = [&last_state](typename ODE::SystemState const& state) {
last_state = state;
};
} else {
append_state = std::bind(
&Ephemeris::AppendMasslessBodiesState, _1, std::cref(trajectories));
}
auto const instance =
parameters.integrator_->NewInstance(problem,
append_state,
tolerance_to_error_ratio,
integrator_parameters);
auto status = instance->Solve(t_final);
// We probably don't care if the vessel gets too close to the singularity, as
// we only use this integrator for the future. So we swallow the error. Note
Instance::Solve(Instant const& t_final) {
using Displacement = typename ODE::Displacement;
using Velocity = typename ODE::Velocity;
using Acceleration = typename ODE::Acceleration;
auto const& a = integrator_.a_;
auto const& b̂ = integrator_.b̂_;
auto const& b̂ʹ = integrator_.b̂ʹ_;
auto const& b = integrator_.b_;
auto const& bʹ = integrator_.bʹ_;
auto const& c = integrator_.c_;
auto& append_state = this->append_state_;
auto& current_state = this->current_state_;
auto& first_use = this->first_use_;
auto& parameters = this->parameters_;
auto const& equation = this->equation_;
// |current_state| gets updated as the integration progresses to allow
// restartability.
// State before the last, truncated step.
std::optional<typename ODE::SystemState> final_state;
// Argument checks.
int const dimension = current_state.positions.size();
Sign const integration_direction = Sign(parameters.first_time_step);
if (integration_direction.Positive()) {
// Integrating forward.
CHECK_LT(current_state.time.value, t_final);
} else {
// Integrating backward.
CHECK_GT(current_state.time.value, t_final);
}
CHECK(first_use || !parameters.last_step_is_exact)
<< "Cannot reuse an instance where the last step is exact";
first_use = false;
// Time step. Updated as the integration progresses to allow restartability.
Time& h = this->time_step_;
// Current time. This is a non-const reference whose purpose is to make the
// equations more readable.
DoublePrecision<Instant>& t = current_state.time;
// Position increment (high-order).
std::vector<Displacement> Δq̂(dimension);
// Velocity increment (high-order).
std::vector<Velocity> Δv̂(dimension);
// Current position. This is a non-const reference whose purpose is to make
// the equations more readable.
std::vector<DoublePrecision<Position>>& q̂ = current_state.positions;
// Current velocity. This is a non-const reference whose purpose is to make
// the equations more readable.
std::vector<DoublePrecision<Velocity>>& v̂ = current_state.velocities;
// Difference between the low- and high-order approximations.
typename ODE::SystemStateError error_estimate;
error_estimate.position_error.resize(dimension);
error_estimate.velocity_error.resize(dimension);
// Current Runge-Kutta-Nyström stage.
std::vector<Position> q_stage(dimension);
// Accelerations at each stage.
// TODO(egg): this is a rectangular container, use something more appropriate.
std::vector<std::vector<Acceleration>> g(stages_);
for (auto& g_stage : g) {
g_stage.resize(dimension);
}
bool at_end = false;
double tolerance_to_error_ratio;
// The first stage of the Runge-Kutta-Nyström iteration. In the FSAL case,
// |first_stage == 1| after the first step, since the first RHS evaluation has
// already occurred in the previous step. In the non-FSAL case and in the
// first step of the FSAL case, |first_stage == 0|.
int first_stage = 0;
// The number of steps already performed.
std::int64_t step_count = 0;
Status status;
Status step_status;
// No step size control on the first step. If this instance is being
// restarted we already have a value of |h| suitable for the next step, based
// on the computation of |tolerance_to_error_ratio_| during the last
// invocation.
goto runge_kutta_nyström_step;
while (!at_end) {
// Compute the next step with decreasing step sizes until the error is
// tolerable.
do {
// Reset the status as any error returned by a force computation for a
// rejected step is now moot.
step_status = Status::OK;
// Adapt step size.
// TODO(egg): find out whether there's a smarter way to compute that root,
// especially since we make the order compile-time.
h *= parameters.safety_factor *
std::pow(tolerance_to_error_ratio, 1.0 / (lower_order + 1));
// TODO(egg): should we check whether it vanishes in double precision
// instead?
if (t.value + (t.error + h) == t.value) {
return Status(termination_condition::VanishingStepSize,
"At time " + DebugString(t.value) +
", step size is effectively zero. Singularity or "
"stiff system suspected.");
}
runge_kutta_nyström_step:
// Termination condition.
if (parameters.last_step_is_exact) {
Time const time_to_end = (t_final - t.value) - t.error;
at_end = integration_direction * h >=
integration_direction * time_to_end;
if (at_end) {
// The chosen step size will overshoot. Clip it to just reach the
// end, and terminate if the step is accepted. Note that while this
// step size is a good approximation, there is no guarantee that it
// won't over/undershoot, so we still need to special case the very
// last stage below.
h = time_to_end;
final_state = current_state;
}
}
auto const h² = h * h;
// Runge-Kutta-Nyström iteration; fills |g|.
for (int i = first_stage; i < stages_; ++i) {
Instant const t_stage =
(parameters.last_step_is_exact && at_end && c[i] == 1.0)
? t_final
: 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) {
Σj_a_ij_g_jk += a[i][j] * g[j][k];
}
q_stage[k] = q̂[k].value + h * c[i] * v̂[k].value + h² * Σj_a_ij_g_jk;
}
step_status.Update(
std::vector<Vector<Acceleration, Frame>>& accelerations) {
Error const error =
ComputeMasslessBodiesGravitationalAccelerations(t,
accelerations.assign(accelerations.size(), Vector<Acceleration, Frame>());
Error error = Error::OK;
// Locking ensures that we see a consistent state of all the trajectories.
absl::ReaderMutexLock l(&lock_);
for (std::size_t b1 = 0; b1 < number_of_oblate_bodies_; ++b1) {
MassiveBody const& body1 = *bodies_[b1];
error |= ComputeGravitationalAccelerationByMassiveBodyOnMasslessBodies<
/*body1_is_oblate=*/true>(
t,
body1, b1,
positions,
accelerations);
}
for (std::size_t b1 = number_of_oblate_bodies_;
b1 < number_of_oblate_bodies_ +
number_of_spherical_bodies_;
++b1) {
MassiveBody const& body1 = *bodies_[b1];
error |= ComputeGravitationalAccelerationByMassiveBodyOnMasslessBodies<
/*body1_is_oblate=*/false>(
t,
GravitationalParameter const& μ1 = body1.gravitational_parameter();
auto const& trajectory1 = *trajectories_[b1];
Position<Frame> const position1 = trajectory1.EvaluatePosition(t);
Length const body1_collision_radius =
mean_radius_tolerance * body1.mean_radius();

template<typename Frame>
Position<Frame> ContinuousTrajectory<Frame>::EvaluatePosition(
Instant const& time) const {
absl::ReaderMutexLock l(&lock_);
CHECK_LE(t_min_locked(), time);

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant