PCI-SPH porting thread

[giovanni.idili]

I am starting this thread to collect questions that may arise when porting the most recent version of the PCI-SPH algorithm to the simulation engine bundle.

First question (for Andrey I think): I noticed that the indexPostPass kernel [line 667] that was being executed via OpenCL in previous versions has been moved to the host code [line 396] in the version I am porting.

Comments say you have done this for performance, which is great, but my question is, are the old version in the kernel and the new version supposed to be doing the exact same thing? Looks like it but just wanted to be sure.

I am asking basically to understand if I can keep the one in OpenCL to speed up porting of the algorithm.

Thanks!

Giovanni 

[s.khayrulin]

Hi Giovanni,

Yes you are right both of this realizations do the same things, but one-thread version (new version) works significantly faster. So now we can remove  indexPostPass kernel from program or comment this.

вторник, 15 января 2013 г., 21:17:59 UTC+7 пользователь giovanni.idili написал:

[giovanni.idili]

Hi Sergey - thanks for this.

New question for you or Andrey: I noticed that there is a new block [line 1291] dealing with elastic matter. This code runs only for the first step.

Can you help me understand what happens in the kernel preElasticMatterPass and in the following function call prepareElasticMatterData() and why do we need to run this only for the first step?

I am asking because I need to understand what happen in the case I need to run the steps "one by one", meaning I run a step, get the results out (and possibly send them to display), then use the results as the initial conditions for the next step. Should we run preElasticMatterPass and prepareElasticMatterData every step in that case?

This is the way things are setup in the current Java version, we call solve and it does one step returning updated values. We will change this to do multiple steps at once in the future. 

Thanks!

[s.khayrulin]

Yes sure, 

Every particle in a system have it is liquid particle or elastic or boundary ("matter index")... Elastic particles form permanent connections between each other, we can define this in initialization, but it is not right. So the definition is producing in first iteration. First of all, we find all neighbour for all particles, so after this operation all matter indexes are disordered and kernel preElasticMatterPass fix this. Why this important? It is because after work of this kernel we run prepareElasticMatterData() that for all ellastic particles finds elastic particles among the founded neighbors and creates the connections between current particle and neighbor elastic particle that is constant:)(I hope you understand). So all seance of this two function it is creates connections between elastic particle. And it should work only once on first iteration of simulation.

Hope this helps.

If you want to realized it more detailed we can make a meeting.

среда, 16 января 2013 г., 23:10:01 UTC+7 пользователь giovanni.idili написал:

[giovanni.idili]

OK - if I understand this correctly then if my simulation is just one step and then I get the results out and wipe all the buffers, and after that I do another simulation from scractch using results of the previous step as the initial conditions for the new simulation, which is also one single step, I would have to run this elastic matter preparation routines every time :)

Will definitely change the sph.solver to run multiple steps at once (as the other solvers do or should do) but still someone could decide to set the number of steps to run to 1 (as in the example above) and in that case this elastic matter preparation stuff will have to run every time a "first step" is calculated, if I understand correctly.

Does this make sense?

[s.khayrulin]

Hi Giovanni,

I'm not quite sure that I understand you enough, I'll discuss this with Andrey. But I have a proposition for you we can remove this elastic matter preparation routines, we can store information about ellastic connections to the file, as we did with positions. And initialize it when we create the particles in simmulation. What do you think?

Thanks,

  Sergey

четверг, 17 января 2013 г., 1:20:59 UTC+7 пользователь giovanni.idili написал:

[Giovanni Idili]

Hi Sergey - if we can actually store that info in the file and treat it same as position I think that would be ideal! Great to hear it's the case :)

--
You received this message because you are subscribed to the Google Groups "OpenWorm-discuss" group.
To post to this group, send email to openworm...@googlegroups.com.
To unsubscribe from this group, send email to openworm-discu...@googlegroups.com.
Visit this group at http://groups.google.com/group/openworm-discuss?hl=en.
 
 

[giovanni.idili]

New question for you guys about this predictive loop [sph.cpp, line 1340]: why is it running exactly 3 times? Also I see that there is code commented out about running until the error becomes less than 2%, is that code obsolete and which version should I port? For now I took the non-commented version (just runs 3 times regardless of error).

Thanks,

Giovanni

[Andrey Palyanov]

"We found a minimum of 3 iterations generally sufficient to achieve
a low level of pressure fluctuations." // B. Solenthaler's dissertation, page 51
I've also observed the same number of iterations mentioned by other
authors of SPH implementations.
I've chosen to use 3 iterations as well to reduce computational
performance and because sometimes for some small amount of particles
the process of iterative corrections seems to be non-converging. I'll
prefer to keep commented code which calculated and checks error value
for further experiments with this.

Best regards,
Andrey

2013/1/23 giovanni.idili <giovann...@gmail.com>:

[Giovanni Idili]

Thanks Andrey for the quick response. I'll keep the 3 iterations too so in the java version - if we ever switch to the error check I'll pick it up.

--
Giovanni Idili - http://www.linkedin.com/in/giovanniidili

[Mike Vella]

From my memory of the code there should be no need to hard code 3 iterations? I would generalize it.

--

[Giovanni Idili]

Yeah will make it configurable as soon as the porting is done. 

[giovanni.idili]

New question for you guys (Andrey, Sergey).

As I am finishing the porting I am integrating new buffers in the host code (C++ in the original version / Java in the simulation engine version), in particular I am looking at the particleIndexBack buffer. This buffer appears many times in the C++ code but is never "used" (filled out with data) other than being passed down to kernel functions. I see it is first used on the sortPostPass kernel function where it seems to be initialized with IDs, and then passed to all the pcisph new kernels that use it retrieving IDs from it (pcisph_computeDensity, pcisph_computeForcesAndInitPressure, pcisph_computeElasticForces, pcisph_predictPositions, pcisph_predictDensity, pcisph_correctPressure, pcisph_computePressureForceAcceleration, pcisph_integrate).

My questions is, is the above correct (the buffer is filled with data in sortPossPass and then that data is used by all the pcisph kernel functions) and also what is the use of the buffer in this new version of the algorithm?

Thanks!

- Giovanni

[giovanni.idili]

From Andrey on my last question:

Hi Giovanni,

Yes, that's correct -- the buffer is filled with data in sortPossPass
and then that data is used by all the pcisph kernel functions. While
sorting particleIndex paricle numbers are shuffled (i-th particle's
coordinates x_i, y_i, z_i are no longer stored in array at index 'i'.
So when after end of the step we need to write new coordinates (which
are stored in 'sortedPositions') to 'positions' array we don't know
what indexes to use to access correct array cells. particleIndexBack
serves to keep this information and use it when necessary. I've
introduced this feature when we started to use elastic matter which
needs constant connections(springs) between pairs of particles (for
liquid this was not necessary).

Best regards,
Andrey

[giovanni.idili]

Hey Sergey,

I am having a look at the refactoring you checked in on the master branch of PCI-SPH, very nice work!

One of the issues I was having was trying to understand how to generate the new model with particles types and elastic matter connection data. You seem to be now doing this (filling buffers) in the helper class, on the loadConfiguration function (Line 72), which is much clearer to understand compared to the previous version we had in the "Integrated Version" branch.

Questions: can you confirm everything to do with particle model generation happens in that method and all the data is coming from config files now (position.txt, velocity.txt, elasticConnections.txt)?

This is exactly what we (Matteo and I) wanted to talk about in the meeting btw I tried to reschedule it for tuesday if you can make it, but I'll keep working on this and maybe there's no need for the meeting now that the code is clearer. I'll let you know.

Thanks,

Giovanni 

[giovanni.idili]

Sergey - after looking at it a bit more I think we are going to need some help with this last piece.

Just to give an idea of what's left to do for the porting to be completed, the createModel method generates this file at the moment (representing a scene with random particle of the same default type as in the old version of the SPH demo), so we need to update createModel to generate the same scene you are generating (worm elastic body splashing into a liquid, I am guessing). 

Any help you can give us with this would greatly accelerate completion of this task - because at the moment the only ways I have to replicate the scene are trying to understand how the scene is generated from the old code (quite hard) or parse the configurations files in the new refactored code (easier but time consuming). Maybe you can share the code you used to generate those files (position.txt, velocity.txt, elasticConnections.txt)? That would probably help a lot.

Thanks!

- Giovanni

[s.khayrulin]

Hi Giovanni!

Sorry for late answer. Honestly I didn't create any universal method for genereting sceen) I use old version of PCI SPH generate start position and put them to the files (position velocity etc.). Yes sure I help you. I can change format of output file and generate for you many configuration in XML format.

Thanks,

 Sergey
пятница, 8 февраля 2013 г., 22:59:17 UTC+7 пользователь giovanni.idili написал:

Sergey - after looking at it a bit more I think we are going to need some help with this last piece.

Just to give an idea of what's left to do for the porting to be completed, the createModel method generates this file at the moment (representing a scene with random particle of the same default type as in the old version of the SPH demo), so we need to update createModel to generate the same scene you are generating (worm elastic body splashing into a liquid, I am guessing). 

Any help you can give us with this would greatly accelerate completion of this task - because at the moment the only ways I have to replicate the scene are trying to understand how the scene is generated from the old code (quite hard) or parse the configurations files in the new refactored code (easier but time consuming). Maybe you can share the code you used to generate those files (position.txt, velocity.txt, elasticConnections.txt)? That would probably help a lot.

Thanks!

- Giovanni

On Friday, February 8, 2013 2:52:58 PM UTC, giovanni.idili wrote:

[giovanni.idili]

Hi Sergey - let's talk about it during a meeting this week then, I have a few questions about the elastic matter connections parameters too. 

Thanks for the offer to help generating xml files with different scenes! I think we should do it from the Java version so we don't write it twice and also so we can regenerate the files if we need to change stuff without asking you every time. We need to extend the xml schema to have particle types and elastic connections (hoping to do it before our meeting), once that's done you can help us generate the xml directly from the Java code (we can have different methods that generate different scenes!). Even just helping us to set up the loops correctly (with different particle types and elastic connections) to generate different scenes would be a great help, once that's in place generating the xml will be very easy :)

Thanks,

- Giovanni

[giovanni.idili]

Hey Andrey, Sergey,

after our last meeting, as discussed, I spent time trying to get a liquid only simulation (including boundary particles) on the ported version of PCI-SPH (while waiting for scripts/code to generate the scene including elastic connections). After some issues to get the kernels to compile on mac-os, I had some success getting something to run (see screenshot with some boundaries and only 2 particles) but there is still something wrong (crashes unexpectedly after a few steps if I add liquid particles or more boundaries) and I suspect it has to do with the initial conditions and/or parameters.

I am creating the boundary particles and a layer of water on the bottom same as Andrey does in this code (IntegratedVersion branch), this is the file that I generate and this is the code I am using to do it (ported from C++ to Java).

My 1st question: is there an assumed relationship between MIN/MAX X/Y/Z (defining the grid dimensions) and the number of particles that needs to be respected ( and could cause the simulation to do weird stuff if not respected)?

We picked: 

X: 0-100

Y: 0-40

Z: 0-40

PARTICLE_COUNT: 1024*42 / 4 = 10752. 

But I noticed that Andrey has the following  in the Integrated branch (sph.h): 

X: 0-120.24

Y: 0-80.16

Z: 0-180.36

PARTICLE_COUNT: 1024*42 = 43008

My 2nd question: why not integers for MAX values as it used to be? I am guessing there is a very specific reason for those decimal values! :) 

Any help appreciated!

- Giovanni

[s.khayrulin]

Hi Giovanny,

I'll try to answer on your questions. Between any neighbour boundary particles should be distance equal to r0 = h * 0.5 (where h is smooth radius h = 3.34) so we need arrange integer number of particle on each plane of sceen (XMAZ / r0 = integer number ) as you can see r0 is float so XMAX can not be integer. Hope it help.

Thanks,

 Sergey

пятница, 15 февраля 2013 г., 8:30:36 UTC+7 пользователь giovanni.idili написал:

Hey Andrey, Sergey,

after our last meeting, as discussed, I spent time trying to get a liquid only simulation (including boundary particles) on the ported version of PCI-SPH (while waiting for scripts/code to generate the scene including elastic connections). After some issues to get the kernels to compile on mac-os, I had some success getting something to run (see screenshot with some boundaries and only 2 particles) but there is still something wrong (crashes unexpectedly after a few steps if I add liquid particles or more boundaries) and I suspect it has to do with the initial conditions and/or parameters.

I am creating the boundary particles and a layer of water on the bottom same as Andrey does in this code (IntegratedVersion branch), this is the file that I generate and this is the code I am using to do it (ported from C++ to Java).

My 1st question: is there an assumed relationship between MIN/MAX X/Y/Z (defining the grid dimensions) and the number of particles that needs to be respected ( and could cause the simulation to do weird stuff if not respected)?

We picked: 

X: 0-100

Y: 0-40

Z: 0-40

PARTICLE_COUNT: 1024*42 / 4 = 10752. 

But I noticed that Andrey has the following  in the Integrated branch (sph.h): 

X: 0-120.24

Y: 0-80.16

Z: 0-180.36

PARTICLE_COUNT: 1024*42 / 4 = 43008

[giovanni.idili]

Hi Sergey - thanks for this, as it would explain why things are broken in my test.

A practical question just to make sure I understand: so I can pick any XMAX, YMAX and ZMAX as long as they give me an integer when they are divided by r0?

Also: Is there any constraint on the total number of particles or can I pick any integer and it should work? I am thinking maybe cannot be > N if the boundary box is too small or maybe it needs to > M so that we have at least enough particles for the boundaries (M, N being any integers).

Thanks for you help!

- Giovanni

[s.khayrulin]

Yes Giovanni, you understand right.

About partcile count: number of boundary particle includes in number of all particle PARTICLE_COUNT = LIQUID_PARTICLE_COUNT + ELASTIC_PARTICLE_COUNT + BOUNDARY_PARTICLE_COUNT, so you can't pick any number of boundary particle. Now I'm working on configurations generator and I take it into account.

Thanks,

 Sergey

суббота, 16 февраля 2013 г., 2:46:09 UTC+7 пользователь giovanni.idili написал:

[Giovanni Idili]

Nice - thanks for confirming about MAX values and for the explanation about PARTICLE_COUNT.

One thing I am still not sure I understand though is if there is a minimum number of boundary particles for a given grid (with the grid being determined by min - max values) I need to have and if so how do I calculate what that minimum is (basically how many boundary particles I need, so that I know PARTICLE_COUNT needs to be more than that). This is probably very simple so I am sorry if I am asking a stupid question :)

And thanks for working on the configurations generator, that will help a lot!

--

You received this message because you are subscribed to the Google Groups "OpenWorm-discuss" group.

To unsubscribe from this group and stop receiving emails from it, send an email to openworm-discu...@googlegroups.com.

To post to this group, send email to openworm...@googlegroups.com.

Visit this group at http://groups.google.com/group/openworm-discuss?hl=en.

For more options, visit https://groups.google.com/groups/opt_out.
 
 

[s.khayrulin]

Hi Giovanni,

Sorry for late response. So you can calculate num of boundary particles from sizes of box(grid - xnin, xmax) as you can see here. But you should remember about particles that locate on  edge because this particle are common for more that one faces.

So total number of boundary particle is : 

      BOUNDARY_PARTICLE_COUNT = 2 *(k*n + n + k +  k*(m-2)+ k + m - 2 + ( n - 2 ) * ( m-2 )+ n + m - 4)

where n = (int)( ( XMAX - XMIN ) / r0 ) //X

m = (int)( ( YMAX - YMIN ) / r0 ) //Y

k = (int)( ( ZMAX - ZMIN ) / r0 ); //Z

Thanks, hope it help

Sergey

суббота, 16 февраля 2013 г., 15:24:34 UTC+7 пользователь giovanni.idili написал:

[giovanni.idili]

Hi Sergey - this helps a lot. I saw that bit of code but the boundary particle count but I didn't understand it was a valid comment, I thought it was just some commented out code! :)

[giovanni.idili]

Hi Sergey - I am trying to work with the new info you gave me but I am slightly confused, let's try with an example:

I picked X, Y, Z MAX all the same for simplicity as 5.01 (with MIN being 0), which means n, m, k are all 3.

According to your formula I expect boundary particle count to be 50:

BOUNDARY_PARTICLE_COUNT = 2 *(k*n + n + k + k*(m-2)+ k + m - 2 + ( n - 2 ) * ( m-2 )+ n + m - 4)

According to the formula on line 1576, boundary particle count should be 54:

BOUNDARY_PARTICLE_COUNT = 2 * n * m + 2 * k * n + 2 * k * m

But, using the code I ported to Java from the C++ implementation to generate boundary particles (I went over it many times and I cannot see any difference from the original), 56 boundary particles are generated!

Question: which of the 3 is the correct boundary particle count in this example? :)

I am pretty sure I need to sort this out before I can hope to see the code running properly.

[s.khayrulin]

I thing second result it's not quite right it's error in formula on line 1576 sory it's my fault. 

BOUNDARY_PARTICLE_COUNT = 2 *(k*n + n + k + k*(m-2)+ k + m - 2 + ( n - 2 ) * ( m-2 )+ n + m - 4) =  2 * n * m + 2 * k * n + 2 * k * m - 4. Can you giv me a link to java file whih boundary particle generation, maybe I'll see a problem and wrong place.

Thanks,

 Sergey

понедельник, 18 февраля 2013 г., 22:40:12 UTC+7 пользователь giovanni.idili написал:

[giovanni.idili]

Thanks Sergey - so it sounds like it should be 50 boundary particles in this example.

Here's the Java function that generates 56 boundary particles (X, Y, Z MAX = 5.01, n, m, k = 3):

https://github.com/openworm/org.openworm.simulationengine.model.sph/blob/pcisphUpdate/src/main/java/org/openworm/simulationengine/model/sph/utils/SPHXMLCreator.java#L53

Would be a interesting test to check if your C++ version generates 50 or 56 :)

[giovanni.idili]

A few updates and a question on this, we have made a lot of progress on this last week in London:

Updates: 

• Integrated the non-deterministic bug-fix from Andrey in the Java version (and also removed the rhoInv buffer no longer used).
• Found out with help from Mike that when loading only boundary particles the crash was due to my error in porting the scene generation code, causing duplicate boundary particles (exact same position) that were crashing it
• Found out that when building scenes with Liquid only the algorithm crashes because liquid particles sometimes end up having the exact same coordinates as some boundary particles (so the duplicate particle issue again). This always happens on the same step for the same scene, so at least it's reproducible, which is a lot, as documented here (last post from me showing particle coordinates).
• We were able to run the "liquid only" scene that Sergey is using (Matteo wrote some quick code to covert txt files to the xml format we are using) to test the C++ version and it does not crash the Java version, meaning we did something right. I will soon write code to convert xml files to txt so that we can run the same scenes :)

Question: 

• there seems to be a bug in the algorithm that causes (liquid) particles to sometime assume the exact same position as the boundary particles, crashing the algorithm at the following step. Are you guys aware of this, have you observed this same behaviour, and / or are investigating?

Let me know what I can do to help narrow it down - unfortunately it's difficult to me to understand what's wrong as I don't understand the algorithm as well as you guys do :)

[giovanni.idili]

I am getting closer to nailing this down - I created an issue for the crash. I have a test that crashes immediately (first step) when running this scene.

As you can see the last particle has coordinates (0, 0, 0) same as one of the boundary particles. If I change that to other coordinates (1, 1, 1) it doesn't crash. The problem is that eventually after a number of steps some other particle will overlap with a boundary particle and it will crash. So at the moment I don't know if it this is allowed and it just shouldn't crash when two particles have the same coordinates... or if they shouldn't get to that point in the first place.

Things that would be helpful for me to know in roder to get to the bottom of this: 

• Is (0, 0, 0) or any other position overlapping exactly with a boundary particle a valid location for a particle or should it never happen?
• Does the C++ version crash if you present a similar scene with particles overlapping boundary particles (just give one of the particles (0, 0, 0) and see if it crashes) or is it the case that the simulation never gets to that state because of different handling of floats / rounding etc.?
• What is the function that determines new position of the particles? Maybe I can look at that code and impose that the overlapping never happens.

Best,

- Giovanni

[giovanni.idili]

As promised here's a video of the latest PCI-SPH Java porting showing weird behavior of elastic matter exploding (blue). Also I am not too sure liquid is behaving as it should either, chances are it's not!

Here's the video: http://www.youtube.com/watch?v=v8_Lw1P1e6Q

P.S> also showcasing the latest advancements on the front-end - good job Matteo!

[Mike Vella]

The exact same thing happened to me when I messed with the timesteps in the c++ implementation.

--

[Giovanni Idili]

Mmm now that you say that Mike, I think I might have an idea of what's going on - our implementation is basically stateless, runs one step and returns the results. I am guessing the current implementation relies on some buffer internally populated that is being wiped clean after every cycle in our implementation  - maybe similar as you were doing when messing with the C++ version?

[Mike Vella]

That sounds plausible. The current implementation only returns every 3rd timestep to the gui afaik, the other two modify state.

On 24 Apr 2013 18:09, "Giovanni Idili" <giovann...@gmail.com> wrote:

Mmm now that you say that Mike, I think I might have an idea of what's going on - our implementation is basically stateless, runs one step and returns the results. I am guessing the current implementation relies on some buffer internally populated that i being wiped clean after every cycle in our implementation  - maybe similar as you were doing when messing with the C++ version?

On Wed, Apr 24, 2013 at 5:53 PM, Mike Vella <vell...@gmail.com> wrote:

--
Giovanni Idili - http://www.linkedin.com/in/giovanniidili

--

[Giovanni Idili]

It gets sampled every 3 steps but some buffers never get wiped clean from what I can see, right? Our internal buffers are completely reallocated every cycle. It's on my TODO list to allow for a given amount of steps to run (instead of one by one) but still whenever it returns it would wipe everything so there still be some version of this problem.

[Mike Vella]

I don't really understand what you mean Giovanni? What information do these internal buffers contain?

--

[Giovanni Idili]

I wish I knew! I am just guessing based on structural differences of the code base. But I hope it's something like that - because it's either that or a needle in a haystack :)

[Mike Vella]

Well things like velocity of each particle need to be preserved for the integrator to function correctly. If you are not maintaining this state between time steps it can't work.

--

[Giovanni Idili]

We are maintaining position and velocity for each particle so that should be fine. Another guess is that all the elastic connections are screwed up so the particles are free to fly away as single particles. I am generating the scene same as the guys did on the C++ version - need time to double check all that's done right too. 

[muffinma...@gmail.com]

hello,

I am searching for a sample implementation of PCI-SPH and Google brought me here as one of the top search results. I looked through some of the comments and postings on Openworm but I couldn't find the source code that you guys mentioned.

in the wiki, the links for the source code is gone
http://code.google.com/p/openworm/wiki/GettingStartedWithPhysics
--> http://code.google.com/p/openworm/downloads/list?q=label:SoftBodySolver

and the descriptions in this page show it only has basic SPH implementations
https://github.com/openworm/Smoothed-Particle-Hydrodynamics/tree/IntegratedVersion

I tried reading through the PCI-SPH paper to understand it myself, but I was stuck on how the author corrected the pressure based. The prediction step of the algorithm was confusing too, so I am searching for some sample code and trying to understand it from there.

Thanks

[Giovanni Idili]

Hi there - this is the latest and greatest PCI-SPH implementation in C++, you can refer to that to understand the algorithm: 

https://github.com/openworm/Smoothed-Particle-Hydrodynamics

We are also working on a Java porting but keep in mind it's still work in progress and it doesn't work quite yet:

https://github.com/openworm/org.openworm.simulationengine.solver.sph

Anything you find on the google code page is obsolete by the way - we are in the process of decommissioning that stuff.

If you have questions on the algorithm itself feel free to start another thread here - you might get help here if you're lucky ;)

--
You received this message because you are subscribed to the Google Groups "OpenWorm-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openworm-discu...@googlegroups.com.
To post to this group, send email to openworm...@googlegroups.com.
Visit this group at http://groups.google.com/group/openworm-discuss?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.