C++

Forum

 
Programm crashes at runtime

Pages: 123
PhysicsIsFun (283)
Good evening again,

here is my complete N-particle simulation of the other thread.
The for-loop in my main() is supposed to do 200 000 iterations, but it only does around 2674 and then crashes at runtime without a proper error message.

What could this be? There should not be a point where I use too much memory, the large array in my calculate_force() function gets locally deleted each step.



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
 
  #include <iostream>
#include <cmath>
#include <fstream>
#include <vector>

using namespace std;

vector<vector <double>> readin_pos(){
    vector <vector <double>> pos;
    ifstream inFile;
    inFile.open("positions.txt");
    if(inFile.fail()){                  // Error test
        cerr<< "Error opening file"<<endl;
        exit(1);
    }
    vector <double> pos_xy;
    double pos_x;
    double pos_y;
    while(!inFile.eof()){
            inFile>> pos_x;
            inFile>> pos_y;
            pos_xy={pos_x, pos_y};
            pos.push_back(pos_xy);
    }
    inFile.close();
    pos.pop_back();
    return pos;
}
vector<vector <double>> readin_vel(){
    vector <vector <double>> vel;
    ifstream inFile;
    inFile.open("velocities.txt");
    if(inFile.fail()){                  // Error test
        cerr<< "Error opening file"<<endl;
        exit(1);
    }
    vector <double> vel_xy;
    double vel_x;
    double vel_y;
    while(!inFile.eof()){
            inFile>> vel_x;
            inFile>> vel_y;
            vel_xy={vel_x, vel_y};
            vel.push_back(vel_xy);
    }
    inFile.close();
    vel.pop_back();
    return vel;
}

void calculate_forces(vector <vector<double>> F, vector <vector<double>> pos, double R, double L, double Epot){ //also calculate E_pot
  double r;
  double rx;
  double ry;
  double f_abs;
  double fx;
  double fy;
  double total_fx;
  double total_fy;
  Epot=0;

  //copy of pos-vector, used to consider images of particles in neighbouring boxes if necessary
  vector <vector <double>> image_pos=pos;

  //matrix storing the 2d forces between all particles
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>[pos.size()];
    for(int j=0;j<pos.size();j++){  //initializing
        F_ij[i][j]={0,0};
    }
  }
  //filling the matrix
  for(int i=0;i<pos.size();i++){
    F_ij[i][i]={0,0};  //no force of particle on itself
    for(int j=0; j<i;j++){
        F_ij[i][j]={-F_ij[j][i].at(0),-F_ij[j][i].at(1)}; //Newton's 3rd law
    }
    for(int j=i+1; j<pos.size();j++){
        rx=image_pos.at(i).at(0)-image_pos.at(j).at(0); //distances between particles
        ry=image_pos.at(i).at(1)-image_pos.at(j).at(1);
        if(rx>=L/2){
                rx-=L;
                image_pos.at(j).at(0)+=L;   //imaging particle to neighbouring box
        }
        else if(rx<=-L/2){
                rx+=L;
                image_pos.at(j).at(0)-=L;
        }
        if(ry>=L/2){
                ry-=L;
                image_pos.at(j).at(1)+=L;
        }
        else if(ry<=-L/2){
                ry+=L;
                image_pos.at(j).at(1)-=L;
        }
        r=sqrt(pow(rx,2)+pow(ry,2));
        if(r<=R){                           //calculate forces
                f_abs=48*pow(r,-14)-24*pow(r,-8);
                fx=f_abs*(image_pos.at(i).at(0)-image_pos.at(j).at(0));
                fy=f_abs*(image_pos.at(i).at(1)-image_pos.at(j).at(1));
                F_ij[i][j]={fx, fy};

                Epot+=4*(pow(r,-12)-pow(r,-6))+1;
        }
    }

  }
//calculate total force vector F, holding the total forces for all particles
for(int i=0;i<pos.size();i++){
    total_fx=0;
    total_fy=0;
    for(int j=0; j<pos.size(); j++){
        total_fx+=F_ij[i][j].at(0);
        total_fy+=F_ij[i][j].at(1);
    }
    F.at(i)={total_fx, total_fy};
}
delete[] F_ij;
}
void update_pos(vector <vector<double>> pos, vector <vector<double>> vel, vector <vector<double>> F, double dt){
for (int i=0;i<pos.size();i++){
    pos.at(i).at(0)+=dt*vel.at(i).at(0)+(1./2)*pow(dt,2)*F.at(i).at(0);
    pos.at(i).at(1)+=dt*vel.at(i).at(1)+(1./2)*pow(dt,2)*F.at(i).at(1);
}
}
void update_vel(vector <vector<double>> vel, vector <vector<double>> F, vector <vector <double>> F_old, double dt){
for (int i=0; i<vel.size();i++){
    vel.at(i).at(0)=vel.at(i).at(0)+(1./2)*dt*(F.at(i).at(0)+F_old.at(i).at(0));
    vel.at(i).at(1)=vel.at(i).at(1)+(1./2)*dt*(F.at(i).at(1)+F_old.at(i).at(1));
}
}
void calculate_Ekin(vector <vector <double>> vel, double Ekin){
Ekin=0;
for(int i=0;i<vel.size();i++){
    Ekin+=(1./2)*(pow(vel.at(i).at(0),2)+pow(vel.at(i).at(1),2));
}
}


int main(){
double dt=0.0005;
double R=pow(2,1/6);
double L=14;
double Ekin;
double Epot;
vector <vector <double>> F_old;  //storage for old forces
vector <vector <double>> pos=readin_pos();
vector <vector <double>> vel=readin_vel();
vector <vector <double>> F(pos.size(), vector <double>(2,0));
calculate_forces(F, pos, R, L, Epot);
calculate_Ekin(vel, Ekin);

//simulation step - velocity-Verlet
for(int i=1; i<=200000;i++){
    update_pos(pos,vel,F,dt);
    F_old=F;
    calculate_forces(F,pos,R,L, Epot);
    update_vel(vel,F,F_old,dt);
    calculate_Ekin(vel, Ekin);
    cout <<i<<endl;
}
}


Regards
jlb (4973)
1
2
3
4
5
6
 
void calculate_Ekin(vector <vector <double>> vel, double Ekin){
   Ekin=0;
   for(int i=0;i<vel.size();i++){
      Ekin+=(1./2)*(pow(vel.at(i).at(0),2)+pow(vel.at(i).at(1),2));
   }
}


Why are you passing Ekin by value into this function just to change it's value to some silly new value.

What is the purpose of this function? You do realize that is is not returning anything to the calling function, right?


Aceix (1118)
Use a debugger to check the kind of signal sent when it crashes. Variable watches will also you more information.

Aceix.
lastchance (6980)
@PhysicsIsFun
Tracing the problem without an input file would be very difficult; you are using, after all, an explicit, rather than implicit, numerical method, so you will need a VERY small timestep and your method may simply be unstable. Your forces, after all, involve huge powers of r.

However, as far as C++ is concerned, many things stand out in your code.

- Why don't you use simple structs for your mathematical vectors like pos and the forces? They only have x and y components after all. It's no good saying you haven't learnt them in C++ yet: it would make for much easier coding if you got to grips with them.

- Reading files (line 19) - don't test for EOF: when you have read the last point you won't actually have set the EOF bit. If this is the reason for pop_back on line 26 then it is not a good one.

- Most of your functions will NOT RETURN ANYTHING. This is because in C++, by default, arguments are passed by VALUE, not by reference: i.e. unless you make them references (by using &) only a copy is sent to the routine; nothing is returned unless by the return value of the function. (If you are used to languages like Fortran the opposite is true, unfortunately.) Thus, for example,
calculate_forces (line 51) will NOT return your forces because F is not passed by reference, neither will update_pos (line 122) return your positions, nor update_vel (line 128) return updated velocities, for exactly the same reason. Nor will Ekin (line 134) return your kinetic energy.

- Line 144 will fail because of integer arithmetic: 1/6 will give 0.

- I may or may not right about this, but I think you may have a massive memory leak if line 120 is insufficient to delete a two-dimensional array. (I could be wrong here - it's not something I attempt very regularly, but I would expect it to be the reverse of the 'new' operations.) Fundamentally, however, you DON'T ACTUALLY NEED ARRAY F_ij - you could work out the total force on any particle and the total system potential energy in the main loop without it.


I really would advise you to strengthen your C++ before you try to tackle quite a complex problem. You need to know about:
- file reading and writing;
- structs (as a minimum);
- passing arguments by value or reference;
- use of arrays and handling of multidimensional arrays (preferably using vector).

Assuming you can get to grips with these first, then you need to check some basic physics. Why don't you print out total kinetic energy and potential energy at the end of each timestep? You have a conservative system (lines 77 and 105 - but I'm still not convinced by the sign of your potential energy), so kinetic energy + potential energy must remain constant: check it.


Last edited on
PhysicsIsFun (283)
@lastchance,
thanks for the detailed answer.

At first: I know that it is kind of stupid to write these large programms without solid basics in C++. The course is supposed to be for C++ beginners but they give zero programming tutorials and we have to do the exercises every week. The physics/programm alone is hard enough, but doing it on top of learning the language from the get go... it's a joke.

My assumptions:
-A void function will not return anything. Changing Ekin in this function is supposed to change the variable Ekin that got passed.
edit: I have just red the tutorial of this site: So I have to give the variable Ekin, e.g., via reference.

-The EOF stuff was taken from Youtube video. My reading functions at least work, I tested them.
-I have no clue what structs are :(

I did print out the energies, but I skipped it here because it was not so important (but I checked them now and it's all messed up. The energies don't change a bit, probably because of your point that my functions don't work^^")
The physics should be fine, though.
Last edited on
PhysicsIsFun (283)
I added the &-sign for the references in the function declarations, as well as a cout command for the calculate_Ekin function. Strangely, Ekin is set to 144 there, which is the size of my vel- and- pos vectors, i.e particle numbers.
I also included the ofstream commands for writing data.

I uploaded the txt. files for the read_in functions.
https://www.file-upload.net/download-12839491/velocities.txt.html
https://www.file-upload.net/download-12839494/positions.txt.html


Since it is night here in Germany, I am offline now, but I would be grateful for further help :)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
 
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>

using namespace std;

vector<vector <double>> readin_pos(){
    vector <vector <double>> pos;
    ifstream inFile;
    inFile.open("positions.txt");
    if(inFile.fail()){                  // Error test
        cerr<< "Error opening file"<<endl;
        exit(1);
    }
    vector <double> pos_xy;
    double pos_x;
    double pos_y;
    while(!inFile.eof()){
            inFile>> pos_x;
            inFile>> pos_y;
            pos_xy={pos_x, pos_y};
            pos.push_back(pos_xy);
    }
    inFile.close();
    pos.pop_back();
    return pos;
}
vector<vector <double>> readin_vel(){
    vector <vector <double>> vel;
    ifstream inFile;
    inFile.open("velocities.txt");
    if(inFile.fail()){                  // Error test
        cerr<< "Error opening file"<<endl;
        exit(1);
    }
    vector <double> vel_xy;
    double vel_x;
    double vel_y;
    while(!inFile.eof()){
            inFile>> vel_x;
            inFile>> vel_y;
            vel_xy={vel_x, vel_y};
            vel.push_back(vel_xy);
    }
    inFile.close();
    vel.pop_back();
    return vel;
}

void calculate_forces(vector <vector<double>>& F, vector <vector<double>>& pos, double R, double L, double& Epot){ //also calculate E_pot
  double r;
  double rx;
  double ry;
  double f_abs;
  double fx;
  double fy;
  double total_fx;
  double total_fy;
  Epot=0;

  //copy of pos-vector, used to consider images of particles in neighbouring boxes if necessary
  vector <vector <double>> image_pos=pos;

  //matrix storing the 2d forces between all particles
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>[pos.size()];
    for(int j=0;j<pos.size();j++){  //initializing
        F_ij[i][j]={0,0};
    }
  }
  //filling the matrix
  for(int i=0;i<pos.size();i++){
    F_ij[i][i]={0,0};  //no force of particle on itself
    for(int j=0; j<i;j++){
        F_ij[i][j]={-F_ij[j][i].at(0),-F_ij[j][i].at(1)}; //Newton's 3rd law
    }
    for(int j=i+1; j<pos.size();j++){
        rx=image_pos.at(i).at(0)-image_pos.at(j).at(0); //distances between particles
        ry=image_pos.at(i).at(1)-image_pos.at(j).at(1);
        if(rx>=L/2){
                rx-=L;
                image_pos.at(j).at(0)+=L;   //imaging particle to neighbouring box
        }
        else if(rx<=-L/2){
                rx+=L;
                image_pos.at(j).at(0)-=L;
        }
        if(ry>=L/2){
                ry-=L;
                image_pos.at(j).at(1)+=L;
        }
        else if(ry<=-L/2){
                ry+=L;
                image_pos.at(j).at(1)-=L;
        }
        r=sqrt(pow(rx,2)+pow(ry,2));
        if(r<=R){                           //calculate forces
                f_abs=48*pow(r,-14)-24*pow(r,-8);
                fx=f_abs*(image_pos.at(i).at(0)-image_pos.at(j).at(0));
                fy=f_abs*(image_pos.at(i).at(1)-image_pos.at(j).at(1));
                F_ij[i][j]={fx, fy};

                Epot+=4*(pow(r,-12)-pow(r,-6))+1;
        }
    }

  }
//calculate total force vector F, holding the total forces for all particles
for(int i=0;i<pos.size();i++){
    total_fx=0;
    total_fy=0;
    for(int j=0; j<pos.size(); j++){
        total_fx+=F_ij[i][j].at(0);
        total_fy+=F_ij[i][j].at(1);
    }
    F.at(i)={total_fx, total_fy};
}
delete[] F_ij;
}
void update_pos(vector <vector<double>>& pos, vector <vector<double>>& vel, vector <vector<double>>& F, double dt){
for (int i=0;i<pos.size();i++){
    pos.at(i).at(0)+=dt*vel.at(i).at(0)+(1./2)*pow(dt,2)*F.at(i).at(0);
    pos.at(i).at(1)+=dt*vel.at(i).at(1)+(1./2)*pow(dt,2)*F.at(i).at(1);
}
}
void update_vel(vector <vector<double>>& vel, vector <vector<double>>& F, vector <vector <double>>& F_old, double dt){
for (int i=0; i<vel.size();i++){
    vel.at(i).at(0)=vel.at(i).at(0)+(1./2)*dt*(F.at(i).at(0)+F_old.at(i).at(0));
    vel.at(i).at(1)=vel.at(i).at(1)+(1./2)*dt*(F.at(i).at(1)+F_old.at(i).at(1));
}
}
void calculate_Ekin(vector <vector <double>>& vel, double& Ekin){
Ekin=0;
for(int i=0;i<vel.size();i++){
    Ekin+=(1./2)*(pow(vel.at(i).at(0),2)+pow(vel.at(i).at(1),2));
}
cout <<Ekin<<endl;
}


int main(){
double dt=0.0005;
double R=pow(2,1./6);
double L=14;
double Ekin;
double Epot;
vector <vector <double>> F_old;  //storage for old forces
vector <vector <double>> pos=readin_pos();
vector <vector <double>> vel=readin_vel();
vector <vector <double>> F(pos.size(), vector <double>(2,0));
calculate_forces(F, pos, R, L, Epot);
calculate_Ekin(vel, Ekin);

//prepare output
ofstream E_kin("E_kin_dt="+to_string(dt)+".dat");
E_kin<< 0<<"  "<< Ekin<<endl;
ofstream E_pot("E_pot_dt="+to_string(dt)+".dat");
E_pot<< 0<<"  "<< Epot<<endl;
ofstream E_tot("E_tot_dt="+to_string(dt)+".dat");
E_tot<< 0<< "  "<<Ekin+Epot<<endl;


//simulation step - velocity-Verlet
for(int i=1; i<=3;i++){
    update_pos(pos,vel,F,dt);
    F_old=F;
    calculate_forces(F,pos,R,L, Epot);
    update_vel(vel,F,F_old,dt);
    calculate_Ekin(vel, Ekin);
    cout <<i<<endl;

    E_kin<< i*dt<<"  "<<Ekin<<endl;
    E_pot<< i*dt<<"  "<<Epot<<endl;
    E_tot<<i*dt<<"  "<<Ekin+Epot<<endl;
}

E_kin.close();
E_pot.close();
E_tot.close();
}

Last edited on
jlb (4973)
1
2
3
4
5
 
void calculate_forces(vector <vector<double>>& F, vector <vector<double>>& pos, double R, double L, double& Epot)  //also calculate E_pot
{
  ...
    //matrix storing the 2d forces between all particles
    vector <double> **F_ij = new vector <double> *[pos.size()];


What is with this try at manual memory allocation? The std::vector was designed to get rid of the manual memory allocation.

PhysicsIsFun (283)
So you would do a 3d-vector instead? Do you think the runtime crash might come from there?


Now with the added &-operators, I realized that my code seems to work for a smaller iteration number.
The total energy of 144 is correct, since the initial conditions were chosen such that each particle (of N=144) holds 1 energy unit. In the beginning, there is no potential energy, so all the energy is stored as kinetic energy. The total energy also stays conserved.

The numerical algorithm is stable.


Still, after around 2000 iterations, the programm crashes at runtime.
ne555 (10692)
here's an idea, post the code that does give you the runtime error. (¿where are those 2000 iterations?)
also, a small explanation about the purpose of your program and the expected output may help.
lastchance (6980)
@PhysicsIsFun,

Could you try to explain (or better, provide a reference) what you are trying to do with your "imaging" in lines 82-97. This will not achieve periodicity, and I could easily make this crash.

Suppose that you had just two particles at
(-7,-7) and (+7,+7)

Then lines 80 and 81 will calculate (give or take the sign).
(rx,ry) = (14,14)

With L=14, your lines 82-97 will then turn this into
(rx,ry) = (0,0)

You then calculate the magnitude of displacement as
r = 0
and when you try to calculate
f_abs=48*pow(r,-14)-24*pow(r,-8);
this will blow up because of the negative powers of r.

I found this out because I couldn't download your particle and velocity files so (after rewriting your code somewhat) I made up some numbers and this was the first that I tried because of your value of L. (Actually I tried 4 particles, distributed one at the centre of each quadrant). It blew up immediately ... but didn't do so when I removed your imaging displacements.

So there is a lot wrong with what you are doing in lines 82-97. They can easily make the program crash. There is also no purpose to the array image_pos: it is not needed. Please confirm exactly what you are trying to do here. It will not impose periodicity (because there is nothing to keep your particles in a finite box). If you want to bounce the particles off the walls then simply reflect them and reverse their velocities at the end of a timestep. (With a little bit of effort and application of change-of-momentum you could even work out the wall pressure as a bonus).

Last edited on
jlb (4973)
So you would do a 3d-vector instead?

Probably not, however I would first need to see the data. I would guess that all of those 2d and 3d vectors could probably be replaced by a vector of a structure that holds the data for each element. But as I said this is only a guess since I haven't seen the actual data.

But by looking at the following function I will say that the vector<vector<double>> should not be required:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
vector<vector <double>> readin_pos(){
    vector <vector <double>> pos;
    ifstream inFile;
    inFile.open("positions.txt");
    if(inFile.fail()){                  // Error test
        cerr<< "Error opening file"<<endl;
        exit(1);
    }
    vector <double> pos_xy;
    double pos_x;
    double pos_y;
    while(!inFile.eof()){
            inFile>> pos_x;
            inFile>> pos_y;
            pos_xy={pos_x, pos_y};
            pos.push_back(pos_xy);
    }
    inFile.close();
    pos.pop_back();
    return pos;
}


Something more like (and please note the comments):
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 
struct Point
{
    double x;
    double y;
};

vector<Point> readin_pos()
{
    vector<Point> position;

    ifstream inFile("positions.txt");

    if(!inFile)// Error test
    {
        cerr << "Error opening file" << endl;
        // Have the calling function handle an empty vector as an error condition. 
        return(position);        // Try to avoid using exit() in a C++ program.
    }

    double pos_x;
    double pos_y;
    
    // Prefer not to use EOF for a read loop, instead use the read operation.
    while(inFile >> pos_x >> pos_y) 
    {
        position.push_back({pos_x, pos_y});
    }

    return position;
}


The reason to avoid exit() is that this C function doesn't properly call the class destructors, which could lead to data corruption.

Still, after around 2000 iterations, the programm crashes at runtime.

Where exactly does it crash? Your debugger should be able to tell you exactly where it detects the problem and allow you to view the data at the time of the crash.

Last edited on
lastchance (6980)
@PhysicsIsFun,

I had a go at rewriting some elements of your code, so some of the below might be useful.

I had to make up some sample data, because I can't read yours.

I do know that I can create initial data that would deliberately crash this code, based on your "imaging" techniques. If I avoid those, then it runs continuously. It is not perfectly conservative (because there is nothing in the code to enforce that): there is a gradual loss of total energy.

If you rearranged your main loop you could do away with the extra F_ij completely.

Please explain what you are trying to do with your "imaging" lines, as per my previous post.

If you want to test the code with some real physics then replace your force system with something for which there is an exact solution (e.g. attractive inverse-square force law between just two particles: there should be a solution which is just circular motion about the mutual centre of mass - e.g. earth and moon).


(Apologies to @jlb - I have just read his post in more detail and appreciate that I have just used exit(1) here.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
 
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <vector>
using namespace std;


struct Point                                     // I think this really would make your code more understandable
{
   double x = 0;                                 // Initialised on creation
   double y = 0;
};


vector<Point> readFile( string filename )        // Use a single reusable routine for reading positions or velocities
{
   vector<Point> p;
   ifstream inFile( filename );   if ( !inFile ) { cerr << "Error opening file\n";  exit( 1 ); }
   double x, y;
   while ( inFile >> x >> y ) p.push_back( { x, y } );
   inFile.close();
   return p;
}


void calculate_forces( vector<Point> &F, const vector<Point> &pos, double R, double L, double &Epot )
{ 
  double r, rx, ry;
  double f_abs;
  int size = pos.size();

  Epot = 0;

  // Matrix storing the 2d forces between all particles     // **** This is NOT needed if you rearrange code
  vector< vector<Point> > F_ij( size, vector<Point>( size ) );

  // Filling the matrix
  for ( int i = 0; i < size; i++ )
  {
     F_ij[i][i] = { 0, 0 };                                 // No force of particle on itself

     for( int j = 0; j < i; j++ )
     {
        F_ij[i][j].x = -F_ij[j][i].x;                       // Newton's 3rd law
        F_ij[i][j].y = -F_ij[j][i].y;
     }

     for ( int j = i + 1; j < size; j++ )
     {
        rx = pos[i].x - pos[j].x;                           // Displacements between particles
        ry = pos[i].y - pos[j].y; 

        //*************************
        // Use of image particles *
        // This is highly dubious *
        //*************************
        if ( rx >=  L / 2 )      rx -= L;
        else if ( rx <= -L / 2 ) rx += L;
        if ( ry >=  L / 2 )      ry -= L;
        else if ( ry <= -L / 2 ) ry += L;

        r = sqrt( rx * rx + ry * ry );
        if ( r <= R )
        {                              // Calculate forces and potential energy
           f_abs = 4.0 * ( 12.0 * pow( r, -14 ) - 6.0 * pow( r, -8 ) );
           F_ij[i][j] = { f_abs * rx, f_abs * ry };
           Epot += 4.0 * (        pow( r, -12 ) -       pow( r, -6 ) ) + 1.0;
        }
     }
   }


   // Calculate total force vector F, holding the total forces for all particles
   // *******************************************************************
   // * You could do this in the main loop, avoiding a separate array F *
   // *******************************************************************
   for ( int i= 0; i < size; i++ )
   {
      F[i] = { 0.0, 0.0 };
      for ( int j = 0; j < pos.size(); j++ )
      {
         F[i].x += F_ij[i][j].x;
         F[i].y += F_ij[i][j].y;
      }
   }

   // F_ij will be removed automatically as it goes out of scope; you don't need to use delete
}


void update_pos( vector<Point> &pos, const vector<Point> &vel, const vector<Point> &F, double dt )
{
   for ( int i = 0; i < pos.size(); i++ )
   {
      pos[i].x += dt * vel[i].x + 0.5 * dt * dt * F[i].x;
      pos[i].y += dt * vel[i].y + 0.5 * dt * dt * F[i].y;
   }
}


void update_vel( vector<Point> &vel, const vector<Point> &F, const vector<Point> &F_old, double dt )
{
   for ( int i = 0; i < vel.size(); i++ )
   {
      vel[i].x = vel[i].x + 0.5 * dt * ( F[i].x + F_old[i].x );
      vel[i].y = vel[i].y + 0.5 * dt * ( F[i].y + F_old[i].y );
   }
}


void calculate_Ekin( const vector<Point> &vel, double &Ekin )      // This would be better as a value-returning function
{
   Ekin = 0;
   for ( int i = 0; i < vel.size(); i++ ) Ekin += 0.5 * ( vel[i].x * vel[i].x + vel[i].y * vel[i].y );
}


int main()
{
   double dt  = 0.0005;
   int    ndt = 100;
   double R   = pow( 2, 1.0 / 6.0  );
   double L   = 14;

   double Ekin, Epot;

   vector<Point> F_old;   
   vector<Point> pos = readFile( "positions.txt" );
   vector<Point> vel = readFile( "velocities.txt" );
   vector<Point> F( pos.size(), { 2, 0 } );

   calculate_forces( F, pos, R, L, Epot );
   calculate_Ekin( vel, Ekin );

   // Simulation step: velocity-Verlet
   for ( int i = 1; i <= ndt; i++ )
   {
       update_pos( pos, vel, F, dt );
       F_old = F;
       calculate_forces( F, pos, R, L, Epot );
       update_vel( vel, F, F_old, dt);
       calculate_Ekin( vel, Ekin );
       cout << i << " " << Ekin << " " << Epot << " " << Ekin + Epot << endl;
   }
}

Last edited on
PhysicsIsFun (283)
Hey guys, thanks for your efforts!

@ne555
The code crashes in the for loop in the main() programm. There is no error message (using windows). The output should be: E_tot= const, E_pot increasing (in the beginning), E_kin decreasing (in the beginning).

@jlb I promise I will learn about structures before I write my next program :))
Thanks, I will also go through your code examples.

@lastchance
I try to explain my thoughts with the image vector.
Suppose the particle, let's call it i, is on the right hand side of the box. Other particles only exerts forces on i if they are within a radius R. Now suppose another particle j is very close to the lefthand side of the box. The difference in position vectors of i and j is larger than R. However, due to periodic boundary conditions, the particle i feels particle j through the walls. It is like the force comes from an image of particle j from the other side of the wall on the righthand side where particle i is.

I did not want to change the positions of the particles here out of fear the vector would change even in the main() function, but as I have learned now, if I give pos by value to calculate_forces(), there should'nt be a problem.
Anyway, I copied all my positions in this image vector. If particle j is in another quadrant of the box than my particle i, I calculated the image position of particle j in the imaginative neighbouring box. For example, if particle j is right next to the lefthand side of the box and on the same height as particle i, the image of j would be directly on the righthand side of the righthand border of the box. The distance between i and the images is then always the closest distance between two particles, so the decision wether the distance is <R has to be made with respect to these images.

That's the (physically correct) idea, but it seems like, according to you, the code does not reproduce that.

It's strange that you can't download the files, I just checked them and I can download them. The positions of the particles are on square lattice points, so the particles are distributed equidistantly over the L x L box. The initial velocities are distributed randomly, such that the total energy is 144 units (= number of particles).
Also, it should be noted that my box goes from [0,L] x [0,L].
The energy must be conserved in this system.

I will go through your code now and try to understand your modifications.
Last edited on
lastchance (6980)
@PhysicsIsFun

I note what you are saying about the image particles and what you are trying to do. However, I don't think there is anything in your code that guarantees to keep your particles in a box
-L/2 <= x <= L/2
-L/2 <= y <= L/2
and without that, it is easy to produce numerical examples that will make your equations blow up. A possible workaround might be to force r to have a (small) minimum positive value before you work out forces.

If you wrote down the equations with calculus then, certainly, total energy would be conserved. However, your update_pos and update_vel routines apply numerical approximations. Whilst energy changes are small, they can't be guaranteed to conserve energy. Believe me: numerical dissipation is a near-universal problem in simulation codes!
PhysicsIsFun (283)
I just reread my post and realized: You are right, I take into account the periodicity for the forces, but never do I specifiy what happens if a particle hits a wall!
So there really was a model error in my code, not only c++ problems. Thanks mate!

Of course, with numerical algorithms, conservation is not ensured 100%. But the velocity-Verlet algorithm on this kind of problem is supposed to be very well behaved.
PhysicsIsFun (283)
Now, back to business.
I have changed my update_position() function so that periodic boundaries are ensured.
1
2
3
4
5
6
7
8
9
10
 
void update_pos(vector <vector<double>>& pos, vector <vector<double>>& vel, vector <vector<double>>& F, double dt, double L){
for (int i=0;i<pos.size();i++){
    pos.at(i).at(0)+=dt*vel.at(i).at(0)+(1./2)*pow(dt,2)*F.at(i).at(0);
    if(pos.at(i).at(0)<0) pos.at(i).at(0)+=L; //periodic boundaries
    else if(pos.at(i).at(0)>=L)pos.at(i).at(0)-=L;
    pos.at(i).at(1)+=dt*vel.at(i).at(1)+(1./2)*pow(dt,2)*F.at(i).at(1);
    if(pos.at(i).at(1)<0) pos.at(i).at(1)+=L;
    else if(pos.at(i).at(1)>=L)pos.at(i).at(1)-=L;
}
}


Otherwise I did not change anything. I really want to understand where the error comes from, so I did not implement the other simplifying stuff you recommended (e.g. structure).
I am still working on making the debugger of my IDE work, but printing out the iteration count of the main loop, it suddenly slows down at around 2783, and stops working at 2786.
How does this sound like to you? Is there not enough memory due to saving the variables wrongly?

edit: I printed out a "success" after each successful running of a subfunction in the main loop.
It turns out the last thing the programm did right was the F_old=F command, so the bug seems to be in the calculate_force() function.

edit²: more accurate: The programm seems to crash in the calculate_force() function during the creation of the F_ij array, so in this part:
1
2
3
4
5
6
7
 
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>[pos.size()];
    for(int j=0;j<pos.size();j++){  //initializing
        F_ij[i][j]={0,0};
    }
  }

pos.size() stays the same the whole time.

Sometimes, there is an error message like: "std::bad_alloc()" or such. So some wrong memory handling here?
Last edited on
lastchance (6980)
The whole idea of using vector<> ... is that you DON'T have to use new!

And for every new there would have to be a corresponding delete. "Corresponding", as in "precisely equal numbers of". So I have little doubt that repeated calls to a routine with imbalanced new and delete will inevitably run out of memory. Your original code had a massive memory leak: you deleted the pointer to the array of rows, but not the individual row pointers: this memory (an entire array's worth for every timestep!) is lost for the remainder of your program run. But you shouldn't be using new anyway. You should be using vectors (and structs) alone.

I think you really should go over to using the much simpler code elements that you have been advised. Few people will have any interest in looking at unnecessarily contorted code. You are asking a lot of people to correct code in a form that they would never dream of writing.
Last edited on
PhysicsIsFun (283)
I understand that it's asked a lot if this code is really impractical. But before I change the code, I just want to see where my error lies, so that I at least learnt from experience that one should use vectors.

The whole idea of using vector<> ... is that you DON'T have to use new!

But in this case, I have an array of vectors, so the 'new' is for the array, right?

And for every new there would have to be a corresponding delete. "Corresponding", as in "precisely equal numbers of". So I have little doubt that repeated calls to a routine with imbalanced new and delete will inevitably run out of memory. Your original code had a massive memory leak: you deleted the pointer to the array of rows, but not the individual row pointers: this memory (an entire array's worth for every timestep!) is lost for the remainder of your program run.

But didn't you say I don't need to delete[ ] in a subfunction because the arrays go out of scope anyway?

How about this then? I only entered (hopefully fitting) delete[] statements.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
void calculate_forces(vector <vector<double>>& F, vector <vector<double>> pos, double R, double L, double& Epot, int i){ //also calculate E_pot
  double r;
  double rx;
  double ry;
  double f_abs;
  double fx;
  double fy;
  double total_fx;
  double total_fy;
  Epot=0;

  //copy of pos-vector, used to consider images of particles in neighbouring boxes if necessary
  vector <vector <double>> image_pos=pos;

  //matrix storing the 2d forces between all particles
//  cout<<"pos.size for i= "<<i<<" : "<<pos.size()<<endl;
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>[pos.size()];
    for(int j=0;j<pos.size();j++){  //initializing
        F_ij[i][j]={0,0};
    }
  }
//  cout <<"creating F_ij possible for i= "<<i<<endl;
  //filling the matrix
  for(int i=0;i<pos.size();i++){
    F_ij[i][i]={0,0};  //no force of particle on itself
    for(int j=0; j<i;j++){
        F_ij[i][j]={-F_ij[j][i].at(0),-F_ij[j][i].at(1)}; //Newton's 3rd law
    }
    for(int j=i+1; j<pos.size();j++){
        rx=image_pos.at(i).at(0)-image_pos.at(j).at(0); //distances between particles
        ry=image_pos.at(i).at(1)-image_pos.at(j).at(1);
        if(rx>=L/2){
                rx-=L;
                image_pos.at(j).at(0)+=L;   //imaging particle to neighbouring box
        }
        else if(rx<=-L/2){
                rx+=L;
                image_pos.at(j).at(0)-=L;
        }
        if(ry>=L/2){
                ry-=L;
                image_pos.at(j).at(1)+=L;
        }
        else if(ry<=-L/2){
                ry+=L;
                image_pos.at(j).at(1)-=L;
        }
        r=sqrt(pow(rx,2)+pow(ry,2));
        if(r<=R){                           //calculate forces
                f_abs=48*pow(r,-14)-24*pow(r,-8);
                fx=f_abs*(image_pos.at(i).at(0)-image_pos.at(j).at(0));
                fy=f_abs*(image_pos.at(i).at(1)-image_pos.at(j).at(1));
                F_ij[i][j]={fx, fy};

                Epot+=4*(pow(r,-12)-pow(r,-6))+1;
        }
    }


  }
//calculate total force vector F, holding the total forces for all particles
for(int i=0;i<pos.size();i++){
    total_fx=0;
    total_fy=0;
    for(int j=0; j<pos.size(); j++){
        total_fx+=F_ij[i][j].at(0);
        total_fy+=F_ij[i][j].at(1);
    }
    F.at(i)={total_fx, total_fy};
}
    for(int i=0;i<pos.size(); i++){
        delete[] F_ij[i];
    }
    delete[] F_ij;

//cout <<"update forces worked for i= "<<i<<endl;
}


Code keeps running now. However it takes a giant amount of time to reach 200 000. Is this to be expected from such a problem? Program is still running...
Last edited on
lastchance (6980)
But didn't you say I don't need to delete[ ] in a subfunction because the arrays go out of scope anyway?


My code DIDN'T use 'new' ... so it DIDN'T need delete. Look carefully!
Your (EDIT: original) code DID use 'new' ... so it DID need delete.

Your (EDIT: original) code didn't have (anything like) enough deletes. It had a thumping great memory hole. EDIT: That appears to be partially fixed.


You are computing O(N^2) particle interaction forces per timestep and doing 200000 timesteps, so yes, it will take some time.

I'd parallelise it!


Last edited on
PhysicsIsFun (283)
My code DIDN'T use 'new' ... so it DIDN'T need delete. Look carefully!
Your (EDIT: original) code DID use 'new' ... so it DID need delete.

Your (EDIT: original) code didn't have (anything like) enough deletes. It had a thumping great memory hole. EDIT: That appears to be partially fixed.

I see, sorry!

Program is done: 8882.232 s, so almost 2,5h.

If I am home, I take a look at the energy plots.
Pages: 123