Commit 69acc484 authored by Erik Strand's avatar Erik Strand

Start preparing an objective function

parent aba8edc2
......@@ -7,31 +7,58 @@ void display() { grapher.display(); }
void idle() { grapher.idle(); }
void mouse(int, int, int, int) { exit(0); }
//--------------------------------------------------------------------------------------------------
struct LaunchObjective {
LaunchObjective() {
masses_.resize(n_bodies_);
masses_[0] = 1e1; // planet
masses_[1] = 1e-1; // satellite
initial_conditions_.resize(2 * dim_ * n_bodies_);
initial_conditions_.segment<dim_>(0 * dim_) = Vector2<Scalar>(0, 0); // position of sun
initial_conditions_.segment<dim_>(1 * dim_) = Vector2<Scalar>(-1, 0); // position of earth
initial_conditions_.segment<dim_>(2 * dim_) = Vector2<Scalar>(0, 0); // velocity of sun
initial_conditions_.segment<dim_>(3 * dim_) = Vector2<Scalar>(0, 0); // velocity of earth
}
void operator()(VectorX<Scalar> const& velocity, Scalar& loss, VectorX<Scalar>& gradient) {
initial_conditions_.segment<dim_>(3 * dim_) = velocity;
nbody_.initialize(dim_, masses_, initial_conditions_, delta_t_);
for (uint32_t i = 0; i < n_steps_; ++i) {
nbody_.step();
}
std::cout << "state jacobian\n";
std::cout << nbody_.state_jacobian() << '\n';
}
constexpr static uint32_t dim_ = 2;
constexpr static uint32_t n_bodies_ = 2;
constexpr static Scalar sim_time_ = 1.4;
constexpr static Scalar delta_t_ = 1e-3;
constexpr static uint32_t n_steps_ = static_cast<uint32_t>(sim_time_ / delta_t_);
VectorX<Scalar> masses_;
VectorX<Scalar> initial_conditions_;
NBody nbody_;
};
//--------------------------------------------------------------------------------------------------
int main() {
constexpr uint32_t dim = 2;
uint32_t const n_bodies = 3;
VectorX<Scalar> masses(n_bodies);
masses[0] = 1.5; // sun
masses[1] = 1; // earth
masses[2] = 0.2; // moon
VectorX<Scalar> initial_conditions(dim * n_bodies * NBodyLaw::d);
Scalar const earth_v_y = -0.3;
Scalar const moon_v_x = 0.3;
Scalar const sun_v_y = -earth_v_y * masses[1] / masses[0];
Scalar const sun_v_x = -moon_v_x * masses[2] / masses[0];
initial_conditions.segment<dim>(0 * dim) = Vector2<Scalar>(0, 0); // position of sun
initial_conditions.segment<dim>(1 * dim) = Vector2<Scalar>(1, 0); // position of earth
initial_conditions.segment<dim>(2 * dim) = Vector2<Scalar>(0, 0.8); // position of moon
initial_conditions.segment<dim>(3 * dim) = Vector2<Scalar>(sun_v_x, sun_v_y); // velocity of sun
initial_conditions.segment<dim>(4 * dim) = Vector2<Scalar>(0, earth_v_y); // velocity of earth
initial_conditions.segment<dim>(5 * dim) = Vector2<Scalar>(moon_v_x, 0); // velocity of moon
NBody nbody;
nbody.initialize(dim, masses, initial_conditions);
LaunchObjective objective;
VectorX<Scalar> satellite_velocity;
satellite_velocity.resize(2);
satellite_velocity[0] = 0;
satellite_velocity[1] = 1;
Scalar loss;
VectorX<Scalar> gradient;
objective(satellite_velocity, loss, gradient);
/*
grapher.initialize(
&nbody,
BoundingBox<Vector2<Scalar>>(
......@@ -44,6 +71,7 @@ int main() {
glutMouseFunc(mouse);
glutIdleFunc(idle);
grapher.run();
*/
return 0;
}
......@@ -87,12 +87,13 @@ void NBodyLaw::operator()(
void NBody::initialize(
uint32_t dim,
VectorX<Scalar> const& masses,
VectorX<Scalar> const& state
VectorX<Scalar> const& state,
Scalar delta_t
) {
assert(2 * dim * masses.size() == state.size());
law_.set_masses(masses);
delta_t_ = 1e-4;
delta_t_ = delta_t;
// initialize state
t_ = 0;
......
......@@ -51,11 +51,17 @@ private:
class NBody : public SimulationInterface {
public:
NBody() {}
void initialize(uint32_t dim, VectorX<Scalar> const& masses, VectorX<Scalar> const& state);
void initialize(
uint32_t dim,
VectorX<Scalar> const& masses,
VectorX<Scalar> const& state,
Scalar delta_t
);
void step() final;
VectorX<Scalar> const& points() const final { return points_; }
VectorX<Scalar> const& state() const { return state_; }
MatrixX<Scalar> const& state_jacobian() const { return state_jacobian_; }
// These methods make NBody compatible with the differentiable integrators.
void operator()(Scalar, VectorX<Scalar> const& state) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment