Simulating a system of moving points

Question

This is an update of the code from my previous question.

I'm simulating a system of particles which interact with each other and move and rotate due to the interaction. There are 2000 particles in total but each particle interacts with other particles and their image (I'm using a periodic boundary condition) if particles and their images have distances larger thanl0. I just consider image of particles in 8 near boxes to the main box.

y[N],x[N] are the position of points in the plane; dx[n] and dy[N] their changes in each time step. c shows different time steps. e_x and e_y shows x and y component of direction of each particle respect to x and y-axis.

#include <iostream>
#include <math.h>
#include <vector>
#include <array>
#include <list>
#include <random>
#include <functional>
#include <fstream>
#include <string>
#include <sstream>
#include <algorithm>
#include <chrono>
using namespace std;

#define pi 3.14159265358979323846

// function declarations(these lines are defined to generate uniform random numbers)
double dist(double x1, double y1, double x2, double y2, double L);
double PBC(double x, double L);
uint64_t s[2] = { 0x41, 0x29837592 };
static inline uint64_t rotl(const uint64_t x, int k) {
    return (x << k) | (x >> (64 - k));
}

uint64_t next(void) {
    const uint64_t s0 = s[0];
    uint64_t s1 = s[1];
    const uint64_t result = s0 + s1;
    s1 ^= s0;
    s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
    s[1] = rotl(s1, 36); // c

    return result;
}

double uniform() {
    return next()*(1.0/18446744073709551616.0);
}

// main function
int main()
{
    unsigned short int N=2000, l0=2;
    double dt=0.001, a=0.001, b=0.002;
    long c_equilib=10;
    double L=10,
        x[N], // x coordinate of r
        y[N], // y coordinate of r
        theta[N]; // angle (direction)     
    double e_x[N], e_y[N];
    for (unsigned short int i=0; i<N; i++)
    {
        x[i]=uniform()*L; // in the range (0, L)
        y[i]=uniform()*L; // in the range (0, L)
        theta[i]=-pi+uniform()*2.0*pi;
        e_x[i]=cos(theta[i]);
        e_y[i]=sin(theta[i]);
    }

    // evolution of time
    long c = 0; //This is number of past time steps in while loop of time
    while (c < c_equilib)
    {
        double dx[N], dy[N],theta_new[N];
        c++;
        for(unsigned short int i = 0; i < N; i++)
        {
            // orientation (theta) change
            unsigned int count = 0;
            double S_theta=0.0, S_x=0.0, S_y=0.0, S;
            for(unsigned short int j=0; j<N; j++)
            {
                if (j==i) continue;
                for(int ky=-1; ky<=1; ky++)
                {
                    for(int kx=-1; kx<=1; kx++)
                    {
                        double r_x_ij=x[i]-(x[j]+kx*L),
                            r_y_ij=y[i]-(y[j]+ky*L),
                        r_ij=sqrt(pow(r_x_ij,2.0)+pow(r_y_ij,2.0));
                        if (r_ij >= l0)
                        { 
                            r_x_ij/=r_ij;
                            r_y_ij/=r_ij;
                            S_theta += (e_x[i]*e_x[j]+e_y[i]*e_y[j]) / pow(r_ij, 3) * ( 3.0* (r_x_ij*e_x[j]+r_y_ij*e_y[j]) *
                                                    pow(e_x[i]*r_y_ij - e_y[i]*r_x_ij, 2.0) -
                                                    e_x[i]*e_y[j] + e_y[i]*e_x[j] );
                            S = ( 3.0* pow(r_x_ij*e_x[j]+r_y_ij*e_y[j], 2.0) -1.0 ) / pow(r_ij, 2.0);
                            S_x += r_x_ij*S;
                            S_y += r_y_ij*S;
                        }
                    }
                }
            }
            theta_new[i] += dt*b*S_theta;
            dx[i] += dt*a*S_x;
            dy[i] += dt*a*S_y;
        }

        // update angles and positions
        for(unsigned short int i=0; i<N; i++)
        {
            theta[i]=theta_new[i];
            // keep theta in the range [-pi, pi]
            while(theta[i]>pi)
                theta[i]-=2*pi;
            while(theta[i]<-pi)
                theta[i]+=2*pi;

            e_x[i]=cos(theta[i]);
            e_y[i]=sin(theta[i]);
            x[i]+=dx[i];
            y[i]+=dy[i];

            // apply periodic boundary condition
            if (x[i] > L)
            {
                x[i]-= L;
            } else if( x[i] < 0) {
            x[i]+= L;
            }
            if (y[i] > L)
            {
            y[i]-=L;
            } else if( y[i] < 0) {
            y[i]+= L;
            }
        }  
    }
    return 0;
}

// calculates the distance between two points
double dist(double x1, double y1, double x2, double y2, double L)
{
    return sqrt(pow(PBC(x1-x2,L),2.0)+pow(PBC(y1-y2,L),2.0));
}

// periodic boundary condition
double PBC(double x, double L)
{
    if(fabs(x)>L/2.0)
        return L-fabs(x);
    else
        return fabs(x);
}

I'm looking to optimize the large loop of this C++ program. This is one part of my simulation (I have ignored some calculations inside the loop to make the question simpler). The main question is his part is extremely time-consuming. How can I make it faster? Second question: How should I choose discretization step size (dt)? could you please give me a good choice for it (comparing with other values in the question) This question is related to the other question which is linked to this one. But, according to comments in the previous question, I have added the runnable code here. As I was asked to post a new question for this purpose by the moderators.

Answer 1

I normally prefer to focus mostly on algorithmic and performance issues, and leave the coding style aspects for others to comment on. In this case, however, it turns out that your code is so hard to read that I can't really tell what your simulation is actually doing. Thus, I'm going to address some of the "superficial" issues that I can see below, in the hope that perhaps the third iteration of your code will actually be readable enough to let people analyze the simulation algorithm itself.

What's that weird code and where did it come from?

At the top of your source code, you have a bunch of "magic" code (specifically, the s array and the rotl() and next() functions) with nothing to indicate what it does or where you got it from.

I happen to be familiar enough with pseudorandom number generators to know one when I see one, and a bit of Googling revealed that your "magic code" is a straightforward copy-paste of the xoroshiro128+ PRNG reference implementation by David Blackman and Sebastiano Vigna. But there's no way to know that just by looking at your code. At a minimum, I'd suggest including a comment like this in your code:

// xoroshiro128+ pseudorandom number generator by Blackman and Vigna (2016)
// code based on http://xoroshiro.di.unimi.it/xoroshiro128plus.c

You should also consider splitting the PRNG code into a separate source file, and perhaps encapsulating it in a class (with the s array as a private member, initialized by the constructor) so that you could have multiple independent PRNG instances if you wanted.

using namespace std;

Many people consider this a bad practice. See this Stack Overflow thread for some arguments why. Consider removing it, and getting into habit of just writing out the std:: prefix where needed.

(If you feel that writing out the std:: prefix everywhere clutters up your code too much and hurts readability, you can also pull specific frequently used names from the std namespace into your current namespace or block scope with using-declarations like using std::string.)

Magic constants

In your uniform() function you have the "magic" constant 1.0/18446744073709551616.0. Again, I can guess from context that it's equal to 2⁻⁶⁴, and that it's used to convert the PRNG output from an unsigned 64-bit integer into a floating-point number between 0 and 1, but you really should at least have a comment in your code saying so.

(Note that, in C++17, you could write that constant using hexadecimal scientific notation simply as 0x1p-64. That syntax has been valid in C since C99, so in practice many compilers supported it in C++ too even before C++17. It probably still deserves a comment, but at least using that notation means you don't have to worry about one of the digits being mistyped.)

For that matter, the uniform() function itself should at least have a comment saying what it does. So should all other functions, too. The next person to read your code (who might be you, a year or ten years from now) will thank you for it.

Also, defining your own π constant with #define pi 3.14159265358979323746 seems like poor style to me. Your compiler almost certainly has one available already, and if not, you can define one yourself using a suitable trigonometric identity like:

static constexpr auto pi = acos(-1);

(Bonus exercise: I deliberately changed one of the digits of π while quoting your code in the previous paragraph. Can you tell which one? And are you sure your original definition is correct, and has enough digits? Can you tell that at a glance?)

More magic constants

At the top of your main() function, you have:

unsigned short int N=2000, l0=2;
double dt=0.001, a=0.001, b=0.002;
long c_equilib=10;
double L=10, // ...

None of those constants have any comments saying what they mean, and the variable names are almost entirely uninformative. I can at least guess that dt is probably a time interval (in some unspecified units), and the following lines suggest that N is the length of some arrays (apparently, based on your accompanying question text, the number of particles in your system), but the rest are pretty mysterious.

You should either rename your constants so that their names are self-explanatory, or add comments explaining what they stand for. Or preferably both. For physical constants, you should also note their units (e.g. does dt = 0.001 mean 0.001 seconds or 0.001 milliseconds or 0.001 days or what?).

For that matter, I can't even be sure that they really are constants, since you didn't include the const keyword. To know for sure, I'd have to go over every single line in your main() function to make sure none of those variables are changed anywhere. The compiler can do that. People reading your code really shouldn't need to.

BTW, l0 is a horrible name for a constant / variable, since in many fonts it's really hard to distinguish from the number 10. Change it.

Why so short?

You've declared your N constant and your loop counters as unsigned short int for no apparent reason. Changing their type to int (or unsigned int) would save you some typing and likely yield marginally more efficient code.

In general, there's rarely any reason to use short integer types for simple variables. They may be useful for saving space in arrays or other dynamic data structures, but in those cases you should normally use the explicitly sized (u)intN_t types anyway.

Don't put everything in main()

All your simulation code is inside your main() function, which is poor style. At least move it into a separate simulate() function instead, and call that from main().

(I know your real code probably doesn't have all the simulation code directly inside main(). At least I really hope it doesn't. But you should still fix that. Also, even when posting simplified code here, it would be nice if you could leave at least some of the output code in place. Right now, running your code just wastes a bunch of CPU cycles and then exists with no output.)

Inconsistent helper function placement

You've used forward declarations to put your dist() and PBC() functions below your main simulation code, but then you've stuck the next() and uniform() functions above it with no forward declaration. Why?

Unused variable

While simplifying your code, you've apparently left in an unsigned int count = 0 variable that isn't used for anything. I noticed because GCC warned me about it when I compiled your code. At least that was the only warning I got with -Wall, which is nice to see.

Comments, comments, comments

Really, your code needs a lot more of them. Add comments explaining what all your variables contain. Add comments explaining what each part of your code does, at a high level. Add comments explaining what your math formulas are supposed to be calculating.

Imagine you were showing your code to someone who had no idea what it's supposed to do (because that's exactly what you're doing here!). Do you really think most people reasonably familiar with C++ could figure out what it does at a glance? If not, you need more comments.

What is your simulation doing, anyway?

It would be a lot easier to follow your code if you actually told us how your particles are supposed to be interacting.

It could also be nice if you could include some kind of graphical output routines so that we could actually see the particles moving, but I can understand that this might result in too much boilerplate. But even just a verbal description would help a lot.

Answer 2

I'll second everything that @Ilmari_Karonen said. In addition, I have the following suggestions.

Avoid Parallel Arrays

Any time you find yourself using parallel arrays (as you have with x, y, theta, e_x, and e_y), that's a sign that you really need to have a data type. OOP is all about creating a data type and methods to act on that data type to separate it out from all the other code you write to make it easier to maintain.

So you should have a Point class (or use an existing one from a library of some sort). I recommend something like this:

class Point {
public:
    Point(double newX, double newY) : x(newX), y(newY) {}
    ~Point() {}
    double getX() const { return x; }
    double getY() const { return y; }
    void setX(double newX) { x = newX; }
    void setY(double newY) { y = newY; }
private:
    double x;
    double y;
};

You might then want the class to contain such useful methods as operator+(), operator-(), operator=(), operator==(), dot(const Point& x).

You might further create class like for a particle. Something like this:

class Particle {
    public:
        Particle(const Point& location, const Point& velocity);
        Particle(const Point& location, double speed, double angle);
        ~Particle();
    private:
        Point loc;
        Point vel;
        double theta;
};

Then you could create an array (or std::vector) of Particle objects. This would allow you to further simplify your code by using range-based loops to iterate over the collection.

Doing all of the above will make reading your code enormously easier.

Performance

As always, the best way to determine why your code is slow is to profile it. Anything else is just a guess. Having said that, I will take an educated guess at what might be the problem, as I've written similar code in the past and hit this problem.

You're using the pow() function to compute squares and cubes. The pow() function is very slow. You should pull out the terms that you want to square and cube and manually do it. Something like this:

double r_ij_squared = r_ij * r_ij;
double r_ij_cubed = r_ij_squared * r_ij;
double er_min = e_x[i] * r_y_ij - e_y[i] * r_x_ij;
double er_min_squared = er_min * er_min;
double er_plus = e_x[j] * r_x_ij + e_y[j] * r_y_ij;
double er_plus_squared = er_plus * er_plus;

You also have some common terms that you are continuously recalculating. It looks like you reuse r_x_ij * e_x[j] + r_y_ij * e_y[j] in the addition to S_theta, then again in the next line when setting S.

These lines can be simplified:

        while(theta[i]>pi)
            theta[i]-=2*pi;
        while(theta[i]<-pi)
            theta[i]+=2*pi;

This is basically a modulus operator. You can do something like:

theta[i] = fmod(theta[i], 2 * pi);
if (theta[i] < 0.0)
{
    theta[i] += 2 * pi;
}

I believe that accomplishes the same thing without a loop. If theta is 27000π, it will execute in the same amount of time as if theta was between 0 and 2π.