Odeint 与 openmp 的集成和并行计算

Question

我正在尝试使用 odeint 在 C++ 中实现三体问题集成，但是要进化，系统需要很多时间（系统集成仅 50 年大约需要 10 小时），所以我想尝试实现并行化脚本的版本。我看到可以使用 OpenMP 完成，并且当有 for 循环时，可以在 for 循环之前添加

 #pragma omp parallel for

。但是在这里我没有任何 for 循环。 我对并行计算和 OpenMP 一无所知。

这是代码：

#include <iostream>
#include <fstream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

// Define the state type as a Boost array of double
typedef std::array< double , 18 > state_type;

double G = 6.67e-11; // gravitational constant
double m1 = 1.8981e27, m2 = 2.40e22, m3 = 2.40e22; // masses of the giove and europa

double RR = 1.496e11; // Normalizing distance in km (= 1 AU)
double MM = 6e24; // Normalizing mass
double TT = 365 * 24 * 60 * 60.0; // nomalizing time (1 year)

double FF = (G * MM * MM) / (RR * RR); // Unit force
double EE = FF * RR; // unit energy   

double GG = (MM*G*TT*TT)/(RR*RR*RR);

double mn1 = m1/MM; // Normalized mass1
double mn2 = m2/MM; // Normalized mass2 
double mn3 = m3/MM; // Normalized mass3

double r12, r13, r23; // distance between bodies
double k = 10e4*RR; // spring constant
double L0 = 1e6/RR; // rest length of the spring

// Declare the energy variable 
double energy = 0;

// Set the output file
ofstream output("output_europa_beta_0.4_final.csv");

// Declare the total energy function
double total_energy( const state_type &x){
    
    double kinetic_energy;
    double potential_energy;

    r12 = sqrt( ( x[0] - x[6] ) * ( x[0] - x[6] ) + ( x[1] - x[7] ) * ( x[1] - x[7] ) + ( x[2] - x[8] ) * ( x[2] - x[8] ) );
    r13 = sqrt( ( x[0] - x[12] ) * ( x[0] - x[12] ) + ( x[1] - x[13] ) * ( x[1] - x[13] ) + ( x[2] - x[14] ) * ( x[2] - x[14] ) );
    r23 = sqrt( ( x[6] - x[12] ) * ( x[6] - x[12] ) + ( x[7] - x[13] ) * ( x[7] - x[13] ) + ( x[8] - x[14] ) * ( x[8] - x[14] ) );

    // kinetic energy
    kinetic_energy = 0.5 * mn1 * (x[3] * x[3] + x[4] * x[4] + x[5] * x[5]) +
                     0.5 * mn2 * (x[9] * x[9] + x[10] * x[10] + x[11] * x[11]) +
                     0.5 * mn3 * (x[15] * x[15] + x[16] * x[16] + x[17] * x[17]);
    
    // gravitational potential energy
    potential_energy  = -GG * mn1 * mn2 / r12 - GG * mn1 * mn3 / r13 - GG * mn2 * mn3 / r23;

    // elastic potential energy
    double dx = x[6] - x[12];
    double dy = x[7] - x[13];
    double dz = x[8] - x[14];
    potential_energy += 0.5 * k * (sqrt( dx * dx + dy * dy + dz * dz) - L0) * (sqrt( dx * dx + dy * dy + dz * dz) - L0);

    return kinetic_energy + potential_energy;
}

// Declare the force function
void three_body_force( const state_type &x , state_type &dxdt , double t )
{
    double beta = 0.4;

    r12 = sqrt( ( x[0] - x[6] ) * ( x[0] - x[6] ) + ( x[1] - x[7] ) * ( x[1] - x[7] ) + ( x[2] - x[8] ) * ( x[2] - x[8] ) );
    r13 = sqrt( ( x[0] - x[12] ) * ( x[0] - x[12] ) + ( x[1] - x[13] ) * ( x[1] - x[13] ) + ( x[2] - x[14] ) * ( x[2] - x[14] ) );
    r23 = sqrt( ( x[6] - x[12] ) * ( x[6] - x[12] ) + ( x[7] - x[13] ) * ( x[7] - x[13] ) + ( x[8] - x[14] ) * ( x[8] - x[14] ) );

    // First body's position and velocity
    dxdt[0] = x[3];
    dxdt[1] = x[4];
    dxdt[2] = x[5];

    // First body's acceleration
    dxdt[3] = ( GG * mn2 * ( x[6] - x[0] ) / ( r12 * r12 * r12 ) ) + ( GG * mn3 * ( x[12] - x[0] ) / ( r13 * r13 * r13 ) );
    dxdt[4] = ( GG * mn2 * ( x[7] - x[1] ) / ( r12 * r12 * r12 ) ) + ( GG * mn3 * ( x[13] - x[1] ) / ( r13 * r13 * r13 ) );
    dxdt[5] = ( GG * mn2 * ( x[8] - x[2] ) / ( r12 * r12 * r12 ) ) + ( GG * mn3 * ( x[14] - x[2] ) / ( r13 * r13 * r13 ) );

    // Second body's position and velocity
    dxdt[6] = x[9];
    dxdt[7] = x[10];
    dxdt[8] = x[11];

    // Second body's acceleration
    dxdt[9] = ( GG* mn1 * ( x[0] - x[6] ) / ( r12 * r12 * r12 ) ) + ( GG * mn3 * ( x[12] - x[6] ) / ( r23 * r23 * r23 ) );
    dxdt[10] = ( GG * mn1 * ( x[1] - x[7] ) / ( r12 * r12 * r12 ) ) + ( GG * mn3 * ( x[13] - x[7] ) / ( r23 * r23 * r23 ) );
    dxdt[11] = ( GG * mn1 * ( x[2] - x[8] ) / ( r12 * r12 * r12 ) ) + ( GG * mn3 * ( x[14] - x[8] ) / ( r23 * r23 * r23 ) );

    // Third body's position and velocity
    dxdt[12] = x[15];
    dxdt[13] = x[16];
    dxdt[14] = x[17];

    // Third body's acceleration
    dxdt[15] = ( GG * mn1 * ( x[0] - x[12] ) / ( r13 * r13 * r13 ) ) + ( GG * mn2 * ( x[6] - x[12] ) / ( r23 * r23 * r23 ) );
    dxdt[16] = ( GG * mn1 * ( x[1] - x[13] ) / ( r13 * r13 * r13 ) ) + ( GG * mn2 * ( x[7] - x[13] ) / ( r23 * r23 * r23 ) );
    dxdt[17] = ( GG * mn1 * ( x[2] - x[14] ) / ( r13 * r13 * r13 ) ) + ( GG * mn2 * ( x[8] - x[14] ) / ( r23 * r23 * r23 ) );

    // Spring force between second and third bodies
    double dx = x[6] - x[12];
    double dy = x[7] - x[13];
    double dz = x[8] - x[14];
    
    double dvx = x[9] - x[15];
    double dvy = x[10] - x[16];
    double dvz = x[11] - x[17];
    
    double spring_force = - k * ( sqrt( dx * dx + dy * dy + dz * dz ) - L0 ) - beta * sqrt(dvx * dvx + dvy * dvy + dvz * dvz);

    dxdt[9] += dx * spring_force / r23; 
    dxdt[10] += dy * spring_force / r23;
    dxdt[11] += dz * spring_force / r23;
    dxdt[15] += -dx * spring_force / r23;
    dxdt[16] += -dy * spring_force / r23; 
    dxdt[17] += -dz * spring_force / r23;

}

// Declare the function that compute the distance of the first body from the line between the second and the third
double distance_from_line(const state_type &x)
{
    double x1 = x[0], y1 = x[1], z1 = x[2]; // position of first body
    double x2 = x[6], y2 = x[7], z2 = x[8]; // position of second body
    double x3 = x[12], y3 = x[13], z3 = x[14]; // position of third body

    double a = (y2 - y3) * z1 - (z2 - z3) * y1 + (z2 * y3 - y2 * z3);
    double b = (z2 - z3) * x1 - (x2 - x3) * z1 + (x2 * z3 - z2 * x3);
    double c = (x2 - x3) * y1 - (y2 - y3) * x1 + (y2 * x3 - x2 * y3);
    double d = sqrt(a * a + b * b + c * c) / sqrt((y2 - y3) * (y2 - y3) + (z2 - z3) * (z2 - z3) + (x2 - x3) * (x2 - x3));
    return d;
}
int i=0;
template<class Stepper>
struct observer {
    Stepper &m_stepper;
    observer(Stepper &stepper) : m_stepper(stepper) {}
    void operator()(const state_type &x, double t) {
        energy = total_energy(x);
        
        // Compute the distance of the first body from the line between the second and the third
        double distance = distance_from_line(x);

        // Write the current state to the file
        if (i%7000000==0|| i==0){
        output  << t << "," << x[0] << "," << x[1] << "," << x[2] << "," << x[3] << "," << x[4] << "," << x[5] << "," << x[6] << "," << x[7] << "," << x[8] << "," << x[9] << "," << x[10] << "," << x[11] << "," << x[12] << "," << x[13] << "," << x[14] << "," << x[15] << "," << x[16] << "," << x[17] << "," << energy << "," << distance <<endl;
        }
        i++;
    }
};


int main()
{
    // Set up the initial conditions
    state_type x = {{ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , -6.647e8/RR , 0.0 , 0.0 , 0.0 , (-13871)*TT/RR , 0.0, -6.657e8/RR , 0.0 , 0.0 , 0.0 , (-13872)*TT/RR, 0.0 }}; // vere distanze giove europa
   
    double dt = 0.06;

    double tend = 50;

    output << "t,x1,y1,z1,vx1,vy1,vz1,x2,y2,z2,vx2,vy2,vz2,x3,y3,z3,vx3,vy3,vz3,energy,distance" << endl;

    typedef runge_kutta_dopri5< state_type > stepper_type;
    auto controlled_stepper = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type());
    observer<decltype(controlled_stepper)> obs(controlled_stepper);

    
    integrate_adaptive( controlled_stepper, three_body_force, x , 0.0 , tend , dt, obs );

    // Close the output file
    output.close();

    return 0;
}

你能帮我减少集成时间并告诉我如何实现这段代码的并行化版本吗？ 非常感谢。