8
\$\begingroup\$

I have wrote the following code in C++ for the n-body problem. This code is sequential as later on I am planning to parallelize it using OpenMP. I want to know whether I have wrote the correct implementation, if there are some errors or bugs, or if this is an invalid approach.

#include <iostream>
#include <ctime>
#include <math.h>
using namespace std;
#define N 100
#define G 6.673e-11
#define TIMESTAMP 1e11
struct Particle{
    double rx, ry;//position components
    double vx, vy;//velocity components
    double fx, fy;//force components
    double mass;//mass of the particle

};
Particle Update(Particle p, double timestamp)
{
    p.vx += timestamp*p.fx / p.mass;
    p.vy += timestamp*p.fy / p.mass;
    p.rx += timestamp*p.vx;
    p.ry += timestamp*p.vy;
    return p;
}
void PrintParticle(Particle p)
{
    printf("rx == %f ry == %f vx == %f vy == %f mass == %f\n", p.rx,p.ry,p.vx,p.vy,p.mass);
}
//Reset the forces on particle
Particle ResetForce(Particle p)
{
    p.fx = 0.0;
    p.fy = 0.0;
    return p;
}
//Add force to particle a by particle b
Particle AddForce(Particle a,Particle b)
{
    double EPS = 3E4;      // softening parameter (just to avoid infinities)
    double dx = b.rx - a.rx;
    double dy = b.ry - a.ry;
    double dist = sqrt(dx*dx + dy*dy);
    double F = (G * a.mass * b.mass) / (dist*dist + EPS*EPS);
    a.fx += F * dx / dist;
    a.fy += F * dy / dist;
    return a;

}

int main()
{
    Particle particles[N];
    srand(time(NULL));
    //randomly generating N Particles
    for (int i = 0; i < N; i++){
        double rx = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].rx = rx;
        double ry = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].ry = ry;
        double vx = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].vx = vx;
        double vy = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].vy = vy;
        double mass = 1.98892e30*rand()*10 + 1e20;
        particles[i].mass = mass;

    }

    int numberofiterations = 10;
    int count = 0;
    while (count < numberofiterations){
        for (int i = 0; i < N; i++)
        {
            particles[i] = ResetForce(particles[i]);
            for (int j = 0; j < N; j++)
            {
                if (i != j)
                {
                    particles[i] = AddForce(particles[i], particles[j]);
                }

            }
        }
        //loop again to update the time stamp here
        for (int i = 0; i < N; i++)
        {
            particles[i] = Update(particles[i], TIMESTAMP);
        }
        for (int i = 0; i < N; i++)
        {
            PrintParticle(particles[i]);
        }
        count++;
    }




    return 0;
}
\$\endgroup\$
  • 2
    \$\begingroup\$ It's usually a bad idea to create a single-threaded algorithm and then adapt it to multi-threading. Unless you can apply some general pattern (like partitioning the problem space and parallelizing the different partitions), you will probably need something quite different to get performant and safe parallel execution. \$\endgroup\$ – Luaan Apr 20 '15 at 9:20
  • 1
    \$\begingroup\$ This is Code Review, so my point is a bit off-topic: You use the explicit Euler method to calculate the next step of your simulation. This is rather unstable and errors accumulate rapidly. You might want to take a look at a Runge-Kutta method. \$\endgroup\$ – M.Herzkamp Apr 20 '15 at 13:04
  • \$\begingroup\$ FYI: here is a nice talk about how to parallelize the n-body problem and this talk uses the n-body problem to demonstrate what you should look out for if you want to speed up your program on a modern processor. \$\endgroup\$ – MikeMB Apr 21 '15 at 22:19
  • \$\begingroup\$ when i run this code it doesn't create an image? \$\endgroup\$ – TMC Nov 24 '16 at 11:46
14
\$\begingroup\$

Here are some observations that may help you improve your code.

Make sure to #include all required headers

This program calls printf and srand but does not include the corresponding headers. Fix that by adding these lines:

#include <cstdio>
#include <cstdlib>

Use objects

You have a Particle structure and then separate functions that operate on Particle data. With only a slight syntax change, you would have a real object instead of C-style code written in C++. For example, you have a loop right now that says this:

for (int i = 0; i < N; i++)
{
    particles[i] = Update(particles[i], TIMESTAMP);
}

With objects, and with C++11 or newer, you could instead write this:

for (auto &p : particles)
    p.update(TIMESTAMP);

All that would be required would be to make Update a member function.

Prefer stream I/O to printf

The printf function was a capable workhorse for many years, but the C++ iostream library is better in a number of regards. Although it's often more typing for the programmer initially, it's better because it has better type checking, less possibility for runtime overhead, and it fits well with the rest of C++. So I would probably change your PrintParticle to this:

friend ostream& operator<<(ostream& out, const Particle &p) {
    return out << 
        "rx == " << p.rx_ <<
        " ry == " << p.ry_ <<
        " vx == " << p.vx_ <<
        " vy == " << p.vy_ <<
        " mass == " << p.mass_ << '\n';
}

Use const where practical

There are a number of things in your code that look and act like constants, such as numberofiterations in main(). Declare them as const for a program more resistant to error. Also, move the constants, such as N into main rather than having them at file scope.

Create a function rather than repeating code

The functions which create the particle use the same rand()-based function several times. Instead of repeating it, make it into a function:

double xrand() {
    return 1e18*exp(-1.8)*(.5 - rand());
}

Don't abuse using namespace std

Putting using namespace std at the top of every program is a bad habit that you'd do well to avoid. It isn't necessarily wrong to use, but be aware of when you absolutely shouldn't do it (such as in header files).

Consider using a better random number generator

If you are using a compiler that supports at least C++11, consider using a better random number generator. In particular, instead of rand, you might want to look at std::uniform_real_distribution and friends in the <random> header.

Consider altering your algorithm

At the moment, the force on each particle from every other particle is calculated exhaustively, but Newton tells us that the force actually acts equally (but opposite) on both particles, so your loop could be simplified and calculations speeded if you take advantage of that fact to only calculate force vectors once and then apply them to both.

Double check your constant values

It seems to me that 3e4 is a somewhat large value for EPS. Additionally, in the calculation of the random mass, the math is currently 1.98892e30*rand()*10 + 1e20, but why not simply write 1.98892e31*rand() + 1e20 instead?

Eliminate return 0 at the end of main

When a C++ program reaches the end of main the compiler will automatically generate code to return 0, so there is no reason to put return 0; explicitly at the end of main.

\$\endgroup\$
  • 2
    \$\begingroup\$ return 0; might be optional but it is confusing that main is the only function in the program that behaves differently in this regard. \$\endgroup\$ – usr Apr 19 '15 at 15:35
  • 3
    \$\begingroup\$ It's confusing after 100 years and it has no value leaving it out. \$\endgroup\$ – usr Apr 19 '15 at 16:16
  • 3
    \$\begingroup\$ @usr: I use the existance of a return at the end of main to indicate if the application reports errors or not. If your application can and does return alterntive values to indicate error situations then a return value is useful other wise it has no meaning so let the compiler generate the correct code for you. \$\endgroup\$ – Martin York Apr 19 '15 at 16:17
  • 1
    \$\begingroup\$ @andrepd: That is a myth (that they are vastly different) and has been covered many times on SO. Just decouple the C++ streams from the C streams std::ios_base::sync_with_stdio(false) and it works as efficiently (slightly slower). see printf more than 5 times faster than std::cout? (Note its not). Also this has been discussed on comp.lang.c++ where the consensus seems to be "be robust". \$\endgroup\$ – Martin York Apr 20 '15 at 18:52
  • 1
    \$\begingroup\$ @andrepd: There is not a 25% difference in printf() and iostreams when it comes to printing floats (on my machine I am getting less than 5%). But yes, if you just want to dump floats and there are billions of them fine use printf (but make sure that the number is significant enough that it is worth it). But if you want more functionality then use streams; it provides much more functionality under the hood (localization/translation etc) in addition to type to checking. \$\endgroup\$ – Martin York Apr 22 '15 at 16:43
10
\$\begingroup\$

All points by @Edvard are good. Here are a couple more.

Encapsulation

All your functions mutate a Particle structure. So why not just make this methods on your class?

struct Particle{
    double rx, ry;//position components
    double vx, vy;//velocity components
    double fx, fy;//force components
    double mass;//mass of the particle

};
Particle Update(Particle p, double timestamp);
void PrintParticle(Particle p);
Particle ResetForce(Particle p);
Particle AddForce(Particle a,Particle b);

I would have written this as:

classe Particle
{
    // Members are all private.
    // You mutate and access them via methods.
    double rx, ry;//position components
    double vx, vy;//velocity components
    double fx, fy;//force components
    double mass;//mass of the particle

    public:

        Particle(double rx, double ry, double vx, double vy, double fx, double fy, double mass)
            :rx(rx) ,ry(ry)
            ,vx(vx) ,vy(vy)
            ,fx(fx) ,fy(fy)
            ,mass(mass)
        {} 

        // If these methods return `void` or `Particle&` depends
        // on context and usage patterns. I have left then as close
        // to the original functions as possible.
        Particle& update(double timestamp);
        void printParticle(std::ostream& str);
        Particle& resetForce();
        Particle& addForce(Particle b);

        friend std::ostream& operator<<(std::ostream& str, Particle const& data)
        {
            printParticle(str);
            return str;
        }
};

Don't use macros

#define N 100
#define G 6.673e-11
#define TIMESTAMP 1e11

Macros don't obey scope rules and have no type information (they are simple text substitutions). Prefer to use const variables (or constexpr in C++11);

constexpr int    N         = 100;
constexpr double G         = 6.673e-11;
constexpr double TIMESTAMP = 1e11;

Prefer Standard containers

Arrays are limited. They don't work well with a lot of other constructs in the language as they quickly degrade into pointers just by looking at them funnily (them people start worrying about memory management).

Particle particles[N];

It is better to use modern constructs like std::vector<> or std::array<> as these object maintain there type information correctly.

std::array<Particle, N>  particles;   // fixed size array of particles.
                                      // useful if you they can be 
                                      // build by default contructor

Alternatively:

 std::vector<Particle>   particles;   // Dynamically sized array.
 particles.reserve(N);                // Useful when built dynamically
                                      // from a loop as you can call
                                      // the constructor specifically.

Prefer for(;;) over while() when counter is used

int count = 0;
while (count < numberofiterations){
  // STUFF
  count++;
}

// Easier to read and reason (as the increment is obvious)
for(int count = 0; count < numberofiterations; ++count) {
    // STUFF
}

Standard algorithms

Look at the standard algorithms and use them when you can. Alternatively use the for() syntax for looping over containers.

Refactoring the main now looks like this:

int main()
{
    srand(time(NULL));

    std::vector<Particle> particles;
    particles.reserve(N);

    //randomly generating N Particles
    for (int i = 0; i < N; i++)
    {
        particles.emplace_back(getRand(), getRand(), 
                               getRand(), getRand(),
                               getRand(), getRand(),
                               getRandMass());
    }

    int constexpr numberofiterations = 10;

    for(int count = 0; count < numberofiterations; ++count) {
        for (int i = 0; i < N; ++i) {
            particles[i].resetForce();
            for (int j = 0; j < N; ++j) {
                if (i != j) {
                    particles[i].addForce(particles[j]);
                }
            }
        }
        //loop again to update the time stamp here
        for (auto& particle: particles) {
            particle.update(TIMESTAMP);
        }
        std::copy(std::begin(particles), std::end(particles),
                  std::istream_iterator<Particle>(std::cout));
    }
}
\$\endgroup\$
  • \$\begingroup\$ The rewrite I did looks almost identical to this except that I also made N and TIMESTAMP local to main, and I missed the opportunity to use emplace_back as you have. Nice! \$\endgroup\$ – Edward Apr 19 '15 at 19:31
  • \$\begingroup\$ he seems to want a normal/Gaussian distribution (bell curve), so std::uniform_real_distribution is not ideal -- try std::normal_distribution instead. en.cppreference.com/w/cpp/numeric/random/normal_distribution \$\endgroup\$ – Snowbody Apr 20 '15 at 17:37
  • \$\begingroup\$ @Snowbody: I did not cover random number generation as Edward had already got that one down. But yes; pick an appropriate distribution. \$\endgroup\$ – Martin York Apr 20 '15 at 18:40
3
\$\begingroup\$

All the advice from @Edvard and @Loki Astari are good. However, one thing I'd also recommend you on doing is to use dimensionless units. That way you have two advantages:

  • It is easier to generalize the system;
  • You are less prone to loss of precision.

@Snowbody (fairly) asked if such system as for \$G=1\$. Yes, but not only. One should set the whole system of units.

Say, if you have the following energy:

\$ U = G \sum_{ij} \frac{m_i m_j}{r_{ij}}\$

You can define dimensionless units as

\$ m' = m/m_0 \$;

\$ r' = r/r_0 \$;

\$ U' = U/E_0 \$.

Thus,

\$U' = \frac{U}{E_0} = \frac{G}{E0} \sum_{ij} \frac{m_i m_j}{r_{ij}} \$.

Replacing \$ r\$ and \$ m \$ in \$ U' \$.

\$U' = \frac{U}{E_0} = \frac{G}{E_0} \sum_{ij} \frac{m_i' m_j' m_0^2}{r_{ij'} r_0} \$.

\$U' = \frac{G m_0^2}{E_0r_0} \sum_{ij} \frac{m_i' m_j'}{r_{ij}'} \$.

Now comes the cool part. We want the first term to unity. So, we say that

\$ \frac{G m_0^2}{E_0r_0} =1 \$

or

\$ \frac{G m_0^2}{r_0} = E_0 \$

\$U' =\sum_{ij} \frac{m_i' m_j'}{r_{ij}} \$.

Now one knows how to link the real values with the simulations. Say one wants to simulate the earth-moon. Then \$ m_0 \$ can be the mass of the Earth. Then, in the program, \$ m_{Earth} = 1.0\$ and \$ m_{Moon} = 0.012300 \$*. One can do the same thing for \$ r_0 \$.

* Moon mass in Earth units according to wikipedia

\$\endgroup\$
  • \$\begingroup\$ Could you expand on what you mean by "dimensionless units" or include a link to same? I believe you mean choosing units such that G=1, right? \$\endgroup\$ – Snowbody Apr 20 '15 at 16:05
  • \$\begingroup\$ Yes, but not only. What you want is a whole system without dimension. I'll edit the main part to add that. \$\endgroup\$ – rmk236 Apr 21 '15 at 16:48
0
\$\begingroup\$
  • why are you storing the force component instead of the acceleration? Do you expect particles' masses to change? Division is usually more expensive than multiplication, so do it as little as possible.
  • Why are you calculating distance squared, then square rooting it, then squaring it again? The compiler may be able to optimize this one though.
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.