search

openworm / sibernetic

This is a C++/OpenCL implementation of the PCISPH algorithm supplemented with a set of biomechanics related features applied to C. elegans locomotion
Other
353 stars 106 forks source link

Choose the best method of integration of equations of movement #19

Closed a-palyanov closed 9 years ago

a-palyanov commented 10 years ago

Choose the best method of integration of equations of movement taking into account recently discovered problem of finite precision of float datatype (more information here: https://github.com/openworm/OpenWorm/issues/158 ) There are many methods available. I've implemented and compared a set of them: Explicit Euler, Implicit Euler, Improved Euler, Velocity Verlet, Runge-Kutta 2nd order and LeapFrog along with analytical solution for a very simple system - single mass point attached to a vertical spring (its another end is attached to the ceiling). LeapFrog should be the best among them, since it has 2nd order of precision and it is symplectic. If I use double values for a, v, and x, it shows excellent results, but for float the problem of finite precision leads to errors which make it completely useless (see image). 17

a-palyanov commented 10 years ago

Fortunately, I've also tried another one symplectic method (1st order) - semi-implicit Euler method. Surprisingly, it does not suffer from finite precision problem - the solution is very close to analytical one at timeStep values up to 1e-04 s (previously we were able to use timeStep up to 5e-6 s). 18 In this simple test it really works at this timeStep, but in a simulation including liquid I see stable behavior of the system only at timeStep values up to 1e-05 s, but even this is twice better than 5e-06. Our recent worm body model and liquid around it works fine with this 1st order semi-implicit Euler method. So, I consider it as a good solution which singificanly improved adequacy of Sibernetic and I plan to update the master branch of the Sibernetic code soon with this new method. Also there is one more variant - we should keep it in mind for the case if even better precision will be necessary. It is 2nd order symplectic LeapFrog algorithm, which works only if we use double values for accelerations, velocities and positions, but right now, in my opinion, it is not necessary.

Neurophile commented 10 years ago

What is the computation cost per step of semi-implicit Euler vs the existing method? The total computation time to advance the simulation by time delta-t is the important metric.

Double precision calculations do take longer than 32 bit, but if we are able to use a simpler numerical method, or to significantly increase the time step, then the trade-off may be worth it.

slarson commented 10 years ago

@a-palyanov believes we have the best option for floats but we may need a different option for doubles

slarson commented 10 years ago

@a-palyanov thanks for the update on this today in the hangout -- feel free to write a bit more to illustrate where this is so other folks can get the update too.

a-palyanov commented 10 years ago

Float_max = 3.402823e+038 Float_min = 1.175494e-038 Double_max = 1.797693e+308 Double_min = 2.225074e-308

Double vs float: double_vs_float

Float: e0074287_48ce64c0b5f24 As we can see at this example, -118.625 is composed of [-1] * [118625] * [1e-03]. So, for fraction part we can work with values from 0 to some MAX_VALUE.

a-palyanov commented 10 years ago

And here is an illustration of how this influences on precision which is available for us in case of calculations (numerical integration). limited precision explanation

a-palyanov commented 10 years ago

My calculations showed that for 32bit float coordinates (in all three dimensions) inside the interval ~ (-100000_r0...100000_r0) (r0 = 'radius' of the SPH particle) the maximal integral error connected with this effect will be limited by ~0.002_r0 per 100 integration steps (or ~0.00002_r0 per 1 integration step), which will not be noticeable. This numbers are for integration time step = 3.00e-05 s. We currently use in simulation 5e-06 s. I hope I will soon recalculate for this value too, but anyway results shouldn't be worse.

Neurophile commented 10 years ago

Remember to ensure that every stage of calculation must stay within the bounds 1E38 1E-38. There are some places in the original kernel that violate this restriction.

Also when summing you need to ensure that the operands are of a similar magnitude. On Jun 19, 2014 12:54 AM, "a-palyanov" notifications@github.com wrote:

My calculations showed that for 32bit float coordinates (in all three dimensions) inside the interval ~ -100000_r0...100000_r0 http://r0%20=%20'radius'%20of%20the%20SPH%20particle the maximal integral error connected with this effect will be limited by ~0.002_r0 per 100 integration steps (or ~0.00002_r0 per 1 integration step), which will not be noticeable.

— Reply to this email directly or view it on GitHub https://github.com/openworm/Smoothed-Particle-Hydrodynamics/issues/19#issuecomment-46526041 .

a-palyanov commented 10 years ago

More precisely, we must stay within the bounds (3.40e+38...1.17e-38, 0.00e+00 , -1.17e-38...-3.40e+38) And I hope in this new code (https://github.com/openworm/Smoothed-Particle-Hydrodynamics/tree/development) which I committed today this restrictions are satisfied.