Programm crashes at runtime
there are a bunch of 'tiny' tweaks you could apply to this code for example integer power issue (pow is significantly slower than X*X for a square and this is generally true for all small integer powers). These may add up given the run time, things that can normally be ignored might shave off 5-10 min!
@jonnin is absolutely right - those calculations with pow() are very expensive, and completely unnecessary. If you write pow(x,a) the computer will probably work this out as exp(a.log(x)) and that combines two very expensive calculations.
In general, write small-integer powers out as simple multiplies: it is many times faster. You should write x * x and never pow(x,2), for example.
Take a look at the force and potential-energy calculations:
Here you use both powers and an unnecessary sqrt evaluation: all of those powers are even powers of r. Provided you have already calculated RSQ=R*R then you will vastly reduce the computational time by something like:
This gets rid of all the expensive pow() calls.
You can remove the power calls in your calculation of kinetic energy as well - you only need velocity squared: do a simple multiply for each component, not a call to pow().
I've also commented before that you simply do not need the separate array F_ij (which you dynamically construct and then delete at each timestep). It is perfectly possible to rearrange your code to evaluate the force on the ith particle without having to separately calculate these interactions.
In general, write small-integer powers out as simple multiplies: it is many times faster. You should write x * x and never pow(x,2), for example.
Take a look at the force and potential-energy calculations:
|
|
Here you use both powers and an unnecessary sqrt evaluation: all of those powers are even powers of r. Provided you have already calculated RSQ=R*R then you will vastly reduce the computational time by something like:
|
|
This gets rid of all the expensive pow() calls.
You can remove the power calls in your calculation of kinetic energy as well - you only need velocity squared: do a simple multiply for each component, not a call to pow().
I've also commented before that you simply do not need the separate array F_ij (which you dynamically construct and then delete at each timestep). It is perfectly possible to rearrange your code to evaluate the force on the ith particle without having to separately calculate these interactions.
Last edited on
| I'd parallelise it! |
What does that mean =(?
| there are a bunch of 'tiny' tweaks you could apply to this code for example integer power issue (pow is significantly slower than X*X for a square and this is generally true for all small integer powers). These may add up given the run time, things that can normally be ignored might shave off 5-10 min! |
Thanks for the tipp, I did not know that.
The energies seem to be well behaved, except for larger timesteps. At dt=0.01 the total energy starts oscillating around its original value and at even larger timesteps (dt=0.05), it explodes to astronomical numbers. I hope this comes from the numerical behaviour of the algorithm, not from some code error.
Now, the next program is supposed to be the same fluid with a thermostat. Usually one could simply extend the previous code, but I will rewrite it and try to incorporate all your hints, especially:
- no arrays, use vectors
- work with structures
- no pow() if not necessary
Again, thank you all for your help, especially lastchance :)
Your plan for reformatting certainly seems sensible. Future development would be much easier if it was in more "modern" C++.
Removing pow() where it is not needed should make a big difference to the run time.
Parallelisation - running on multiple processors - depends on you having access to suitable computer facilities and being aware of a suitable paradigm: MPI or openMP. Useful where you are doing fundamentally the same operation on a large number of entities (as here), but rather a steep learning curve.
Good luck with it.
Removing pow() where it is not needed should make a big difference to the run time.
Parallelisation - running on multiple processors - depends on you having access to suitable computer facilities and being aware of a suitable paradigm: MPI or openMP. Useful where you are doing fundamentally the same operation on a large number of entities (as here), but rather a steep learning curve.
Good luck with it.
I have two questions to the following code proposed by lastchance:
Firstly, I don't know that kind of writing the while-condition. Does it mean "as long as inFile >> x >> y is possible, do..."?
Secondly: p is a vector of elements of type point.
When you say
|
|
Firstly, I don't know that kind of writing the while-condition. Does it mean "as long as inFile >> x >> y is possible, do..."?
Secondly: p is a vector of elements of type point.
When you say
p.push_back( { x, y } ) will C++ automatically see that x and y are the members of the structure? Is it necessary for them to be exactly the same variables as in the structure definition, or could I also write p.push_back( { a, b } )?
while ( inFile >> x >> y )Tries to read x and y successively from the stream. If it manages to (because there is enough data) the stream stays "good" and, when required, this can be converted to a boolean true.
When there isn't enough data the stream becomes "bad" (there's a more formal definition) but its status can be converted to a boolean false. 'true' and 'false' then determine the while loop operation.
It is a standard idiom for "read while there is still data to be read" - better than testing for EOF.
If p is a vector of "things" then
p.push_back( { x, y } ) will try to construct an object of type "thing" from x and y. {x,y} is just a standard, implicit, member-based constructor. Actually, you used it in your code when you wrote F_ij[i][j] = { fx, fy };
Last edited on
Alright, thanks :)
Now, I have compressed my code. When it took 400+s before, it takes only 20+s now.
However, taking a really close look at the plots of both programms, they seem to differ a little bit (they show qualitatively same behaviour and have the same order of magnitude).
Could this come from different data processing / calculating in the new code or is there something wrong?
Plot from new code: http://www.directupload.net/file/d/4927/pdzyh8d5_png.htm
Plot from old code: http://www.directupload.net/file/d/4927/ta7dy25h_png.htm
Now, I have compressed my code. When it took 400+s before, it takes only 20+s now.
|
|
However, taking a really close look at the plots of both programms, they seem to differ a little bit (they show qualitatively same behaviour and have the same order of magnitude).
Could this come from different data processing / calculating in the new code or is there something wrong?
Plot from new code: http://www.directupload.net/file/d/4927/pdzyh8d5_png.htm
Plot from old code: http://www.directupload.net/file/d/4927/ta7dy25h_png.htm
Could you put your plots somewhere where one doesn't have to click on external html files: paranoia and all that. I'm certainly not taking the risk.
There's something wrong with your velocity updates - they don't need a factor of 1/2 on the RHS:
change in velocity = acceleration x time .... you only needed a factor of 1/2 when you were including the "old" time-level forces as well (which you were doing in your previous version: not quite sure why you've changed from a second-order scheme to a first-order one). Your acceleration is just force / mass (with mass=1, presumably). The velocity error will throw out your kinetic energy, obviously. Position updates look OK ("ut + 0.5 a t2")
You have eliminated a lot of the pow() calls, which is good - but you can easily eliminate all of them.
Write 0.5 rather than 1./2. - even if an optimising compiler deals with it, it would still be easier to read.
There's something wrong with your velocity updates - they don't need a factor of 1/2 on the RHS:
change in velocity = acceleration x time .... you only needed a factor of 1/2 when you were including the "old" time-level forces as well (which you were doing in your previous version: not quite sure why you've changed from a second-order scheme to a first-order one). Your acceleration is just force / mass (with mass=1, presumably). The velocity error will throw out your kinetic energy, obviously. Position updates look OK ("ut + 0.5 a t2")
You have eliminated a lot of the pow() calls, which is good - but you can easily eliminate all of them.
Write 0.5 rather than 1./2. - even if an optimising compiler deals with it, it would still be easier to read.
Last edited on
You really need to work on your indentation, as presented your code is difficult to read, especially in your calculations.
I would also recommend that you consider working on removing code duplication. For example:
This code is not only horribly formatted but it has a lot of code duplication. With a little better formatting you should be able to better see the duplication.
You should be able to see that the math is duplicated for both x and y. You should be able to reduce the duplication by using an additional function that does the actual math.
Also I recommend you use parentheses when switching between addition/subtraction and multiplication/division for clarity.
Lastly since one of your complaints is execution speed you may want to try to eliminate the need to use the at() function when accessing those vectors. An alternative would be to check the sizes before the loop and throw() an exception if the three vectors are not the same size.
Something like the following:
I would also recommend that you consider working on removing code duplication. For example:
|
|
This code is not only horribly formatted but it has a lot of code duplication. With a little better formatting you should be able to better see the duplication.
|
|
You should be able to see that the math is duplicated for both x and y. You should be able to reduce the duplication by using an additional function that does the actual math.
Also I recommend you use parentheses when switching between addition/subtraction and multiplication/division for clarity.
Lastly since one of your complaints is execution speed you may want to try to eliminate the need to use the at() function when accessing those vectors. An alternative would be to check the sizes before the loop and throw() an exception if the three vectors are not the same size.
Something like the following:
|
|
| Could you put your plots somewhere where one doesn't have to click on external html files: paranoia and all that. I'm certainly not taking the risk. |
but where? Is there an upload function in the forum?
| There's something wrong with your velocity updates - they don't need a factor of 1/2 on the RHS: change in velocity = acceleration x time .... you only needed a factor of 1/2 when you were including the "old" time-level forces as well (which you were doing in your previous version: not quite sure why you've changed from a second-order scheme to a first-order one). Your acceleration is just force / mass (with mass=1, presumably). |
The formular is v=v+ 0.5*dt²*(F+F_old) [edit: square over dt is wrong]
So in order to do it without keeping F_old somewhere, I just used: v=+0.5*dt² *F
So in main it goes: update v --> update F --> update v, which should work just as earlier.
@jlb You're right, thanks for the hint!
When exactly can I use [ ] instead of at()? I have red somewhere that it was safer to use at( ), so not knowing what I'm doing, I just used it.
Last edited on
@PhysicIsFun,
In that case your velocity update is wrong.
Either use
v -> v + 0.5 * dt * ( F + Fold )
or
v -> v + dt * F
but DON'T just form the completely wrong hybrid
v -> v + 0.5 * dt * F
because all accelerations will be wrong by a factor of 2!
You would have been better (second-order time accuracy) keeping Fold - the extra resources are negligible.
You don't need .at() anywhere: your code will be faster and more readable with normal array notation [].
In that case your velocity update is wrong.
Either use
v -> v + 0.5 * dt * ( F + Fold )
or
v -> v + dt * F
but DON'T just form the completely wrong hybrid
v -> v + 0.5 * dt * F
because all accelerations will be wrong by a factor of 2!
You would have been better (second-order time accuracy) keeping Fold - the extra resources are negligible.
You don't need .at() anywhere: your code will be faster and more readable with normal array notation [].
I thought .at and [] were the same result.
vector foo;
foo[10];
foo.at(10); //same thing, I think? I honestly havent used .at, but that was my understanding of it, and I would hope the compiler does the same thing for both. [] is considered more readable.
at() may do an unnecessary bounds check, which may be 'safer' but safety costs time for stuff like that and its totally unnecessary if you wrote the code correctly. Not sure, you can google at() if you want, I am at work and busy atm.
parallelism and some tweaks can be at odds... for example if you decide to make some of the internals into function calls, making any temporary values static can be faster (it depends) but that isn't thread safe. Still, better to do 2 threads if you can; an obvious split would be to fill the X vector on one thread and the Y vector on another, with much of the code the same across the 2 cpus.
as you take a couple more passes over the thing ... common sluggishness in code problems come from a few oft-seen things...
- doing the same work more than once, as noted already
- using a slow algorithm or process (here pow was an example)
- memory allocations issues (push-back can incur a full copy of the whole vector every so often, this is huge over time) or things like new/delete pairs inside a tight loop, etc.
- calling functions in a tight loop where the function can't be inlined or has an issue (parameter data copy mistake, IE passed a large thing by value, allocation of memory, too many local variables created each call, etc)
- and the old 'you didnt even need to do that at all' stuff... unnecessary or redundant comparisons, unnecessary assignment to temporaries / intermediates, and so on.
also human knowledge tweaks can be huge. Lets take a basic 'did this hit that' 3-d object check. You want to know if the distance between this object and that one is less than some constant, like the rough radius from the center of a more or less round object, for example. You can do the distance formula, sqrt of stuff, but you can ALSO compare the squared values without it to get the same answer! You don't need the actual distance, you need to know if they HIT, and to know that, the sqrt isnt needed! Such understanding of the problem can greatly reduce the work done on things. This falls into the rubber duck debugging as someone put it :) If you can't explain why you needed to compute a particular value to be used or output, maybe you didnt need it.
vector foo;
foo[10];
foo.at(10); //same thing, I think? I honestly havent used .at, but that was my understanding of it, and I would hope the compiler does the same thing for both. [] is considered more readable.
at() may do an unnecessary bounds check, which may be 'safer' but safety costs time for stuff like that and its totally unnecessary if you wrote the code correctly. Not sure, you can google at() if you want, I am at work and busy atm.
parallelism and some tweaks can be at odds... for example if you decide to make some of the internals into function calls, making any temporary values static can be faster (it depends) but that isn't thread safe. Still, better to do 2 threads if you can; an obvious split would be to fill the X vector on one thread and the Y vector on another, with much of the code the same across the 2 cpus.
as you take a couple more passes over the thing ... common sluggishness in code problems come from a few oft-seen things...
- doing the same work more than once, as noted already
- using a slow algorithm or process (here pow was an example)
- memory allocations issues (push-back can incur a full copy of the whole vector every so often, this is huge over time) or things like new/delete pairs inside a tight loop, etc.
- calling functions in a tight loop where the function can't be inlined or has an issue (parameter data copy mistake, IE passed a large thing by value, allocation of memory, too many local variables created each call, etc)
- and the old 'you didnt even need to do that at all' stuff... unnecessary or redundant comparisons, unnecessary assignment to temporaries / intermediates, and so on.
also human knowledge tweaks can be huge. Lets take a basic 'did this hit that' 3-d object check. You want to know if the distance between this object and that one is less than some constant, like the rough radius from the center of a more or less round object, for example. You can do the distance formula, sqrt of stuff, but you can ALSO compare the squared values without it to get the same answer! You don't need the actual distance, you need to know if they HIT, and to know that, the sqrt isnt needed! Such understanding of the problem can greatly reduce the work done on things. This falls into the rubber duck debugging as someone put it :) If you can't explain why you needed to compute a particular value to be used or output, maybe you didnt need it.
Last edited on
| When exactly can I use [ ] instead of at()? |
Basically you can use the array notation any time you use the at() function. The at() function insures that you will never access the vector out of range, at the cost of slightly slower access. The array notation doesn't check for out of bounds indexes and is therefore slightly faster than the at() function. However in your loop, the pos vector will never be accessed out of bounds since you're using the size() function to terminate the loop. However before the loop you will need to insure that the other two vectors are of equal or larger sizes to prevent out of bounds access.
| lastchance wrote: |
|---|
| In that case your velocity update is wrong. Either use v -> v + 0.5 * dt * ( F + Fold ) or v -> v + dt * F but DON'T just form the completely wrong hybrid v -> v + 0.5 * dt * F because all accelerations will be wrong by a factor of 2! You would have been better (second-order time accuracy) keeping Fold - the extra resources are negligible. You don't need .at() anywhere: your code will be faster and more readable with normal array notation []. |
This time, I can't follow you...
v + 0.5 * dt * ( F + Fold ) = v + 0.5 * dt * F + 0.5 * dt * F_old
So lets say I have a current velocity of v, and a current force F.
update_vel --> new velocity u = v + 0.5 *dt * F
update_F --> the F above is basically F_old now
update_vel --> new velocity w = u + 0.5 * dt * F = v + 0.5 * dt * F_old + 0.5 * dt * F = v + 0.5 * dt *(F_old +F)
Right?
@jonnin, jlb Thanks for the insight. I will remove the at( ) stuff.
OK - I missed the fact that you were updating velocities twice per timestep: my error and my apologies.
Last edited on
A couple of other small things.
For your periodic boundary conditions, I think you are now dealing with this correctly in update_pos(). In calculate_forces() it is reasonable to adjust the rx, ry displacements for a force coming from a hypothetical particle in the adjacent box ... but I don't think you should adjust the particle positions at this point, because a later ith particle will find this jth particle at a completely different position, which is inconsistent. So
should just be
This will significantly affect your function call, because you can now pass the particle positions by constant reference rather than value:
At present, passing by value means the program has to copy the entire, large contents of that array. Provided you don't change the positions here - and I don't think you should - it can simply pass an implicit pointer to the start of the vector. This should give a big performance boost. (I'd simply call them pos(), not image_pos() as well.)
Finally, it's merely a suggestion, but you are simulating a free molecular gas by applying periodic boundary conditions (if something goes out of the right it comes back in at the left.) Just a thought, but you could, instead, just keep the same particles in the box by bouncing them elastically off the walls (reflect their position in the wall and reverse the corresponding velocity component if they go out). The really nice thing about this is that they will have a change of momentum -2mv, so giving an impulse to the wall of +2mv. If you total this up over the course of the calculation, then divide by time and the area of the wall you will get the pressure P in the box. The average kinetic energy per particle should be proportional to absolute temperature T and you have the volume of the box V. By playing around with these you should be able to simulate something akin to the ideal gas law PV=nRT (chemistry) or PV=NkT (physics). Just a thought.
For your periodic boundary conditions, I think you are now dealing with this correctly in update_pos(). In calculate_forces() it is reasonable to adjust the rx, ry displacements for a force coming from a hypothetical particle in the adjacent box ... but I don't think you should adjust the particle positions at this point, because a later ith particle will find this jth particle at a completely different position, which is inconsistent. So
|
|
if(rx>=L/2) rx-=L; // effect of a particle in a neighbouring box. This will significantly affect your function call, because you can now pass the particle positions by constant reference rather than value:
void calculate_forces(vector<particle> &F, const vector<particle> &pos, double R, double L, double& Epot)At present, passing by value means the program has to copy the entire, large contents of that array. Provided you don't change the positions here - and I don't think you should - it can simply pass an implicit pointer to the start of the vector. This should give a big performance boost. (I'd simply call them pos(), not image_pos() as well.)
Finally, it's merely a suggestion, but you are simulating a free molecular gas by applying periodic boundary conditions (if something goes out of the right it comes back in at the left.) Just a thought, but you could, instead, just keep the same particles in the box by bouncing them elastically off the walls (reflect their position in the wall and reverse the corresponding velocity component if they go out). The really nice thing about this is that they will have a change of momentum -2mv, so giving an impulse to the wall of +2mv. If you total this up over the course of the calculation, then divide by time and the area of the wall you will get the pressure P in the box. The average kinetic energy per particle should be proportional to absolute temperature T and you have the volume of the box V. By playing around with these you should be able to simulate something akin to the ideal gas law PV=nRT (chemistry) or PV=NkT (physics). Just a thought.
Last edited on
Yes, I realize that it is completely unnecessary to do change image_pos. I just need the positions to calculate rx and ry and for the if() conditions, so I will remove it.
But I learnt something again, that copying by value takes more effort than by reference. Sounds obvious, but without you saying it, I wouldn't have thought about it.
Yeah, this time it is all about periodic boundaries because we want to have a homogeneous system. And since it is so small with N=144, the surface interactions would disturb that.
What do you think of my earlier question with the slight change in plots?
Can it, in principle happen that ineffective code yields slightly different results than efficient code? I know you did not see the plots, but they show the same qualitative behaviour.
The energy for semi long dt fluctuates around the same value, the fluctuation is equally strong, only some peaks are a bit thinner or broader or something like that, so you can see a small difference with blanc eye.
But I learnt something again, that copying by value takes more effort than by reference. Sounds obvious, but without you saying it, I wouldn't have thought about it.
Yeah, this time it is all about periodic boundaries because we want to have a homogeneous system. And since it is so small with N=144, the surface interactions would disturb that.
What do you think of my earlier question with the slight change in plots?
Can it, in principle happen that ineffective code yields slightly different results than efficient code? I know you did not see the plots, but they show the same qualitative behaviour.
The energy for semi long dt fluctuates around the same value, the fluctuation is equally strong, only some peaks are a bit thinner or broader or something like that, so you can see a small difference with blanc eye.
Last edited on
I think I would reserve judgement until the position changes have been removed from the force calculations and the simulation rerun. There are some other checks you can do with your code, like no net drift (average velocity is zero), or you could deliberately put a drift velocity in your initial conditions and check that it is maintained (conservation of momentum).
I should make the advised changes to your code, including switching to normal array notation, rather than using .at(), and re-post.
There are a number of websites that you can use to paste images, including imgur.com
I should make the advised changes to your code, including switching to normal array notation, rather than using .at(), and re-post.
There are a number of websites that you can use to paste images, including imgur.com
Hi, sorry for being inactive the last week, was busy with other stuff.
Now we have talked about the property of systems being chaotic like this N particle system. Chaoticity means small changes in the input variables yield large deviations in the end results - independent of the numerical algorithm. This means that the same given initial conditions will lead to slightly different end results when the code runs on different engines, OR when one solves the problem with an inefficient algorithm compared to an efficient algorithm. So that I see small deviations was to be expected.
Now I need to expand my program to model a Nosé-Hoover-Thermostat. Sounds like fun.
Thank you all for the help up until here.
Now we have talked about the property of systems being chaotic like this N particle system. Chaoticity means small changes in the input variables yield large deviations in the end results - independent of the numerical algorithm. This means that the same given initial conditions will lead to slightly different end results when the code runs on different engines, OR when one solves the problem with an inefficient algorithm compared to an efficient algorithm. So that I see small deviations was to be expected.
Now I need to expand my program to model a Nosé-Hoover-Thermostat. Sounds like fun.
Thank you all for the help up until here.
FWIW the same ideas are used in hashing, encryption and security algorithms. You want John and Jon and Jan to all three give *completely* different results though the data changes are smallish. There are plenty of techniques to get this result. If they gave similar results, it would be easier to reverse engineer / hack the security stuff and for hashing, similar data would cluster into similar results causing collisions and destroying the whole point of the technique.