Skip to content

Commit

Permalink
Merge pull request #3518 from eggrobin/no-vectors
Browse files Browse the repository at this point in the history 
Get rid of the obnoxious vectors
  • Loading branch information
@eggrobin
eggrobin authored Jan 29, 2023
2 parents 3b6b829 + 18d2fd1 commit fefdd4f
Show file tree
Hide file tree
Showing 12 changed files with 89 additions and 206 deletions.
67 changes: 14 additions & 53 deletions integrators/embedded_explicit_runge_kutta_integrator_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,28 +106,10 @@ Solve(typename ODE::IndependentVariable const& s_final) {
DependentVariableDerivatives f;
DependentVariableDerivatives last_f;
std::vector<DependentVariableDifferences> k(stages_);
for (auto& k_stage : k) {
for_all_of(ŷ, Δŷ, k_stage).loop([](auto const& ŷ, auto& Δŷ, auto& k_stage) {
int const dimension = ŷ.size();
Δŷ.resize(dimension);
k_stage.resize(dimension);
});
}

for_all_of(ŷ, error_estimate, f, last_f, y_stage)
.loop([](auto const& ŷ,
auto& error_estimate,
auto& f,
auto& last_f,
auto& y_stage) {
int const dimension = ŷ.size();
error_estimate.resize(dimension);
f.resize(dimension);
last_f.resize(dimension);
for (auto const& ŷₗ : ŷ) {
y_stage.push_back(ŷₗ.value);
}
});
for_all_of(ŷ, y_stage).loop([](auto const& ŷ, auto& y_stage) {
y_stage = ŷ.value;
});

bool at_end = false;
double tolerance_to_error_ratio;
Expand Down Expand Up @@ -202,34 +184,21 @@ Solve(typename ODE::IndependentVariable const& s_final) {
// TODO(phl): Should dimension |Σⱼ_aᵢⱼ_kⱼ| in the not FSAL case.
DependentVariableDifferences Σⱼ_aᵢⱼ_kⱼ{};
for (int j = 0; j < i; ++j) {
for_all_of(k[j], ŷ, y_stage, Σⱼ_aᵢⱼ_kⱼ)
.loop([&a, i, j](auto const& kⱼ,
auto const& ŷ,
auto& y_stage,
auto& Σⱼ_aᵢⱼ_kⱼ) {
int const dimension = ŷ.size();
Σⱼ_aᵢⱼ_kⱼ.resize(dimension);
for (int l = 0; l < dimension; ++l) {
Σⱼ_aᵢⱼ_kⱼ[l] += a(i, j) * kⱼ[l];
}
for_all_of(k[j], Σⱼ_aᵢⱼ_kⱼ)
.loop([&a, i, j](auto const& kⱼ, auto& Σⱼ_aᵢⱼ_kⱼ) {
Σⱼ_aᵢⱼ_kⱼ += a(i, j) * kⱼ;
});
}
for_all_of(ŷ, Σⱼ_aᵢⱼ_kⱼ, y_stage)
.loop([](auto const& ŷ, auto const& Σⱼ_aᵢⱼ_kⱼ, auto& y_stage) {
int const dimension = ŷ.size();
for (int l = 0; l < dimension; ++l) {
y_stage[l] = ŷ[l].value + Σⱼ_aᵢⱼ_kⱼ[l];
}
y_stage = ŷ.value + Σⱼ_aᵢⱼ_kⱼ;
});

termination_condition::UpdateWithAbort(
equation.compute_derivative(s_stage, y_stage, f), step_status);
}
for_all_of(f, k[i]).loop([h](auto const& f, auto& kᵢ) {
int const dimension = f.size();
for (int l = 0; l < dimension; ++l) {
kᵢ[l] = h * f[l];
}
kᵢ = h * f;
});
}

Expand All @@ -244,16 +213,11 @@ Solve(typename ODE::IndependentVariable const& s_final) {
auto& Σᵢ_b̂ᵢ_kᵢ,
auto& Σᵢ_bᵢ_kᵢ,
auto& error_estimate) {
int const dimension = ŷ.size();
Σᵢ_b̂ᵢ_kᵢ.resize(dimension);
Σᵢ_bᵢ_kᵢ.resize(dimension);
for (int l = 0; l < dimension; ++l) {
Σᵢ_b̂ᵢ_kᵢ[l] += b̂[i] * kᵢ[l];
Σᵢ_bᵢ_kᵢ[l] += b[i] * kᵢ[l];
Δŷ[l] = Σᵢ_b̂ᵢ_kᵢ[l];
auto const Δyₗ = Σᵢ_bᵢ_kᵢ[l];
error_estimate[l] = Δyₗ - Δŷ[l];
}
Σᵢ_b̂ᵢ_kᵢ += b̂[i] * kᵢ;
Σᵢ_bᵢ_kᵢ += b[i] * kᵢ;
Δŷ = Σᵢ_b̂ᵢ_kᵢ;
auto const Δy = Σᵢ_bᵢ_kᵢ;
error_estimate = Δy - Δŷ;
});
}
tolerance_to_error_ratio =
Expand All @@ -275,10 +239,7 @@ Solve(typename ODE::IndependentVariable const& s_final) {
// Increment the solution with the high-order approximation.
s.Increment(h);
for_all_of(Δŷ, ŷ).loop([](auto const& Δŷ, auto& ŷ) {
int const dimension = ŷ.size();
for (int l = 0; l < dimension; ++l) {
ŷ[l].Increment(Δŷ[l]);
}
ŷ.Increment(Δŷ);
});

RETURN_IF_STOPPED;
Expand Down
44 changes: 19 additions & 25 deletions integrators/embedded_explicit_runge_kutta_integrator_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ double HarmonicOscillatorToleranceRatio(
Speed const& v_tolerance,
std::function<void(bool tolerable)> callback) {
auto const& [position_error, velocity_error] = error;
double const r = std::min(q_tolerance / Abs(position_error[0]),
v_tolerance / Abs(velocity_error[0]));
double const r = std::min(q_tolerance / Abs(position_error),
v_tolerance / Abs(velocity_error));
callback(r > 1.0);
return r;
}
Expand Down Expand Up @@ -140,9 +140,9 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest,
}
{
auto const& [positions, velocities] = solution.back().y;
EXPECT_THAT(AbsoluteError(x_initial, positions[0].value),
EXPECT_THAT(AbsoluteError(x_initial, positions.value),
IsNear(7.5e-2_(1) * Metre));
EXPECT_THAT(AbsoluteError(v_initial, velocities[0].value),
EXPECT_THAT(AbsoluteError(v_initial, velocities.value),
IsNear(6.8e-2_(1) * Metre / Second));
EXPECT_EQ(t_final, solution.back().s.value);
EXPECT_EQ(steps_forward, solution.size());
Expand Down Expand Up @@ -175,9 +175,9 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest,
}
{
auto const& [positions, velocities] = solution.back().y;
EXPECT_THAT(AbsoluteError(x_initial, positions[0].value),
EXPECT_THAT(AbsoluteError(x_initial, positions.value),
IsNear(2.1e-1_(1) * Metre));
EXPECT_THAT(AbsoluteError(v_initial, velocities[0].value),
EXPECT_THAT(AbsoluteError(v_initial, velocities.value),
IsNear(6.8e-2_(1) * Metre / Second));
EXPECT_EQ(t_initial, solution.back().s.value);
EXPECT_EQ(steps_backward, solution.size() - steps_forward);
Expand Down Expand Up @@ -241,12 +241,12 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest, MaxSteps) {
StatusIs(termination_condition::ReachedMaximalStepCount));
EXPECT_THAT(AbsoluteError(
x_initial * Cos(ω * (solution.back().s.value - t_initial)),
positions[0].value),
positions.value),
IsNear(5.9e-2_(1) * Metre));
EXPECT_THAT(AbsoluteError(
-v_amplitude *
Sin(ω * (solution.back().s.value - t_initial)),
velocities[0].value),
velocities.value),
IsNear(5.6e-2_(1) * Metre / Second));
EXPECT_THAT(solution.back().s.value, Lt(t_final));
EXPECT_EQ(40, solution.size());
Expand All @@ -268,9 +268,9 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest, MaxSteps) {
auto const outcome = instance->Solve(t_final);
auto const& [positions, velocities] = solution.back().y;
EXPECT_THAT(outcome, StatusIs(termination_condition::Done));
EXPECT_THAT(AbsoluteError(x_initial, positions[0].value),
EXPECT_THAT(AbsoluteError(x_initial, positions.value),
IsNear(7.5e-2_(1) * Metre));
EXPECT_THAT(AbsoluteError(v_initial, velocities[0].value),
EXPECT_THAT(AbsoluteError(v_initial, velocities.value),
IsNear(6.8e-2_(1) * Metre / Second));
EXPECT_EQ(t_final, solution.back().s.value);
EXPECT_EQ(steps_forward, solution.size());
Expand Down Expand Up @@ -308,8 +308,8 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest, Singularity) {
ODE::DependentVariableDerivatives& dependent_variable_derivatives) {
auto const& [q, v] = dependent_variables;
auto& [qʹ, vʹ] = dependent_variable_derivatives;
qʹ[0] = v[0];
vʹ[0] = mass_flow * specific_impulse / mass(t);
qʹ = v;
vʹ = mass_flow * specific_impulse / mass(t);
return absl::OkStatus();
};
InitialValueProblem<ODE> problem;
Expand All @@ -326,8 +326,8 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest, Singularity) {
ODE::State const& /*state*/,
ODE::State::Error const& error) {
auto const& [position_error, velocity_error] = error;
return std::min(length_tolerance / Abs(position_error[0]),
speed_tolerance / Abs(velocity_error[0]));
return std::min(length_tolerance / Abs(position_error),
speed_tolerance / Abs(velocity_error));
};

AdaptiveStepSizeIntegrator<ODE> const& integrator =
Expand All @@ -340,12 +340,12 @@ TEST_F(EmbeddedExplicitRungeKuttaIntegratorTest, Singularity) {
parameters);
auto const outcome = instance->Solve(t_final);

auto const& [positions, velocities] = solution.back().y;
auto const& [position, velocity] = solution.back().y;
EXPECT_THAT(outcome, StatusIs(termination_condition::VanishingStepSize));
EXPECT_EQ(101, solution.size());
EXPECT_THAT(solution.back().s.value - t_initial,
AlmostEquals(t_singular - t_initial, 4));
EXPECT_THAT(positions.back().value,
EXPECT_THAT(position.value,
AlmostEquals(specific_impulse * initial_mass / mass_flow, 155));
}

Expand Down Expand Up @@ -514,16 +514,10 @@ void PrintTo(
typename internal_embedded_explicit_runge_kutta_integrator::ODE::
State const& state,
std::ostream* const out) {
auto const& [positions, velocities] = state.y;
auto const& [position, velocity] = state.y;
*out << "\nTime: " << state.s << "\n";
*out << "Positions:\n";
for (int i = 0; i < positions.size(); ++i) {
*out << " " << i << ": " << positions[i] << "\n";
}
*out << "Velocities:\n";
for (int i = 0; i < velocities.size(); ++i) {
*out << " " << i << ": " << velocities[i] << "\n";
}
*out << "Position: " << position << "\n";
*out << "Velocity: " << velocity << "\n";
}

} // namespace internal_ordinary_differential_equations
Expand Down
42 changes: 8 additions & 34 deletions integrators/explicit_linear_multistep_integrator_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,23 +63,12 @@ ExplicitLinearMultistepIntegrator<Method, ODE_>::Instance::Solve(

// Current stage of the integrator.
DependentVariables y_stage;
for_all_of(y, y_stage).loop([](auto const& y, auto& y_stage) {
int const dimension = y.size();
y_stage.resize(dimension);
});

absl::Status status;

while (h <= (s_final - s.value) - s.error) {
DoubleDependentVariables Σⱼ_minus_αⱼ_yⱼ{};
DependentVariableDerivatives Σⱼ_βⱼ_numerator_fⱼ{};
for_all_of(y_stage, Σⱼ_minus_αⱼ_yⱼ, Σⱼ_βⱼ_numerator_fⱼ)
.loop(
[](auto const& y, auto& Σⱼ_minus_αⱼ_yⱼ, auto& Σⱼ_βⱼ_numerator_fⱼ) {
int const dimension = y.size();
Σⱼ_minus_αⱼ_yⱼ.resize(dimension);
Σⱼ_βⱼ_numerator_fⱼ.resize(dimension);
});

// Scan the previous steps from the most recent to the oldest. That's how
// the Adams-Bashforth formulæ are typically written.
Expand All @@ -100,11 +89,8 @@ ExplicitLinearMultistepIntegrator<Method, ODE_>::Instance::Solve(
auto const& fⱼ,
auto& Σⱼ_minus_αⱼ_yⱼ,
auto& Σⱼ_βⱼ_numerator_fⱼ) {
int const dimension = yⱼ.size();
for (int i = 0; i < dimension; ++i) {
Σⱼ_minus_αⱼ_yⱼ[i] -= Scale(αⱼ, yⱼ[i]);
Σⱼ_βⱼ_numerator_fⱼ[i] += βⱼ_numerator * fⱼ[i];
}
Σⱼ_minus_αⱼ_yⱼ -= Scale(αⱼ, yⱼ);
Σⱼ_βⱼ_numerator_fⱼ += βⱼ_numerator * fⱼ;
});
}

Expand All @@ -118,11 +104,7 @@ ExplicitLinearMultistepIntegrator<Method, ODE_>::Instance::Solve(
for_all_of(Σⱼ_βⱼ_numerator_fⱼ, Σⱼ_minus_αⱼ_yⱼ)
.loop([h, β_denominator](auto const& Σⱼ_βⱼ_numerator_fⱼ,
auto& Σⱼ_minus_αⱼ_yⱼ) {
int const dimension = Σⱼ_βⱼ_numerator_fⱼ.size();
for (int i = 0; i < dimension; ++i) {
Σⱼ_minus_αⱼ_yⱼ[i].Increment(h * Σⱼ_βⱼ_numerator_fⱼ[i] /
β_denominator);
}
Σⱼ_minus_αⱼ_yⱼ.Increment(h * Σⱼ_βⱼ_numerator_fⱼ / β_denominator);
});

auto& yₙ₊₁ = Σⱼ_minus_αⱼ_yⱼ;
Expand All @@ -131,12 +113,8 @@ ExplicitLinearMultistepIntegrator<Method, ODE_>::Instance::Solve(
auto& y,
auto& y_stage,
auto& current_step_yʹ) {
DCHECK_EQ(y_stage.size(), yₙ₊₁.size());
current_step_yʹ.resize(y_stage.size());
for (int i = 0; i < y_stage.size(); ++i) {
y_stage[i] = yₙ₊₁[i].value;
y[i] = yₙ₊₁[i];
}
y_stage = yₙ₊₁.value;
y = yₙ₊₁;
});
current_step.y = yₙ₊₁;
termination_condition::UpdateWithAbort(
Expand Down Expand Up @@ -177,13 +155,9 @@ FillStepFromState(ODE const& equation,
step.s = state.s;
step.y = state.y;
typename ODE::DependentVariables y;
for_all_of(state.y, y, step.yʹ)
.loop([](auto const& state_y, auto& y, auto& step_yʹ) {
step_yʹ.resize(state_y.size());
for (auto const& state_yᵢ : state_y) {
y.push_back(state_yᵢ.value);
}
});
for_all_of(state.y, y).loop([](auto const& state_y, auto& y) {
y = state_y.value;
});
// Ignore the status here. We are merely computing yʹ to store it, not to
// advance an integrator.
equation.compute_derivative(step.s.value, y, step.yʹ).IgnoreError();
Expand Down
8 changes: 4 additions & 4 deletions integrators/explicit_linear_multistep_integrator_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,8 +183,8 @@ TEST_P(ExplicitLinearMultistepIntegratorTest, Convergence) {
ODE::DependentVariableDerivatives& dependent_variable_derivatives) {
auto const& [q, v] = dependent_variables;
auto& [qʹ, vʹ] = dependent_variable_derivatives;
qʹ[0] = v[0];
vʹ[0] = mass_flow * specific_impulse / mass(t);
qʹ = v;
vʹ = mass_flow * specific_impulse / mass(t);
return absl::OkStatus();
};
InitialValueProblem<ODE> problem;
Expand All @@ -205,8 +205,8 @@ TEST_P(ExplicitLinearMultistepIntegratorTest, Convergence) {
integrator.NewInstance(problem, append_state, step);
EXPECT_OK(instance->Solve(t_final));
Time const t = final_state.s.value - t_initial;
Length const& q = std::get<0>(final_state.y)[0].value;
Speed const& v = std::get<1>(final_state.y)[0].value;
Length const& q = std::get<0>(final_state.y).value;
Speed const& v = std::get<1>(final_state.y).value;
double const log_q_error = std::log10(RelativeError(
q,
specific_impulse *
Expand Down
Loading

0 comments on commit fefdd4f

Please sign in to comment.