• FEATURED
• FEATURED
• FEATURED
• FEATURED
• FEATURED

View more

View more

View more

### Image of the Day Submit

IOTD | Top Screenshots

### The latest, straight to your Inbox.

Subscribe to GameDev.net Direct to receive the latest updates and exclusive content.

# Simple SPH questions. Kernel interpolation - *(its starting to flow)

Old topic!

Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

13 replies to this topic

### #1Peacekeeper  Members

Posted 30 March 2014 - 02:01 PM

Hi!

Iam trying to develop a very basic 2d SPH program in python and it seems that in this forum there are actually some people who know more about this topic^^

I have read a ton of literature but as English is not my main language im running into a few problems here. Right now I created a basic SPH program that calculates the movement based purely on the pressure gradients.  I think my main problem is the understanding of the basic kernel interpolation.

I see the same formulas everywhere but sometimes I wish someone supplied a basic example with values to compare with.

to calculate the density of a specific particle i need to use this:

1. 2D Kernels:

I don’t really understand what they normalization factors (15/pih^6) are for and no kernel that is supplied online does actually equal  to “1”. Some are even called the same but the volume created by the kurve (for 2D) is not even close to one.  Also i found different formulas for the same kernel. the Spiky Kernel can be written as “factor * (h-r)^3” or  “factor * (1-(r/h))^3”.

“factor * (h-r)^3” with the factor being “10/(PI()*I3^5)” does actually provide the same volume even if “h” is changed. So I assumed it should be correct. The value is actually ~ 4.6 (not  1!)

I tried it in excel for 2d! i hope i did it the right way. maybe my whole understanding of how it works is wrong...

As I use different Kernels for calculating the density (poly6) and pressure (spiky) I think the totally different values do fu** up my calculations. Using poly6 my density estimation is around “455” while using the spiky kernel leads to densitys in the range of millions! This brings me to my second question

2. People talk about using the rest density ( for  “p = (p – p0)” ) like a real density! They use 1000 kg/m³ for it. But isn’t the density basically a pseudo density? How can I compute realistic densities and values out of a summation of weighted mass points? Like I said – my rest density is highly influenced by the kernel and is between 100 and 10000000

attachments:

- simple excel sheet in rar. maybe someone can edit it to be correct?

it goes all boom:

no viscosity so kinetic energy is preserved?

hm?

Edited by Peacekeeper, 31 March 2014 - 07:02 PM.

### #2Javier Meseguer de Paz  Members

Posted 31 March 2014 - 05:07 AM

I don’t really understand what they are for and no kernel that is supplied online does actually equal  to “1”.

How are you computing the volume? Because the volume of the kernel is its integral in the range [-h,h], and it has to be 1.

About 2: what you can do (at least that's what I did) was to set the mass of your particles so that their SPH rest density is your desired rest density. it's fairly easy, basically you do:

• create 1 particle P with a full neighborhood filled of particles around it
• set mass of all particles to one.
• compute the density of P, call it invMass.
• Since the density is proportional to the mass of your particles, set the mass of your particles is then 1/invMass * rest_density.

That way, your rest density is going to be the desired one.

PS: About going boom. SPH is highly unstable, really. Well, at least it is more unstable the more you want to conserve volume. In my (little) experience you'll need a well-tuned viscous force to compensate for your pressure forces to avoid explosions. I used Monaghan's.

Edited by Javier Meseguer de Paz, 31 March 2014 - 07:14 AM.

“We should forget about small efficiencies, say about 97% of the time; premature optimization is the root of all evil” -  Donald E. Knuth, Structured Programming with go to Statements

"First you learn the value of abstraction, then you learn the cost of abstraction, then you're ready to engineer" - Ken Beck, Twitter

### #3Peacekeeper  Members

Posted 31 March 2014 - 08:53 AM

thanks for the input! the idea with the "invMass" makes sense to me and i will try it. but it still doesnt answer how people can archieve physical values ?

regarding the volume of the kernel. i hacked it together in excel like this:

the idea was to multiply the surface below one side of the function with "2* PI * r" - (Rotational Symmetry)

1. i discretized my values from min (e.g. "0") to max (e.g. "1") into 50 points.

2. calculated W(r,h) with the formula at all points

3. simple little rectangles for the integral ( r0 * W(r0,h) + (r1-r0) * W(r1,h) + (r2-r1) * W(r1,h) + (r3-r2) .....)

4. multiplied the sum with the 2 * PI * r with r beeing equal to h.

i agree that the excel sheet is a bit unclear ^^ i attached a slightly improved version and here is a gif so you dont have to download it.

when you look at the red box for the volume it stays the same around ~4.6

#### Attached Files

Edited by Peacekeeper, 31 March 2014 - 08:56 AM.

### #4Javier Meseguer de Paz  Members

Posted 31 March 2014 - 03:16 PM

Concerning the integral, I have attached another excel page so that you can verify the Poly6 kernel in 2d.

As for how to use physical quantities in SPH well, I don't understand your question...

I told you how to adjust the mass of the particles so that SPH tries to conserve the correct density of the fluid via the state equation (and the forces that arise from there...) In that regard, the SPH is going to be correct.

The mass of the particles indeed is going to depend on the kernel. However, keep in mind that your particles are ultimately samples, not "physical balls". So for example, if you want to know the mass of a box of fluid of 1x1x1 meters a SPH, you don't actually add up the masses of the particles that are inside1. You get the particles that are inside, you get their densities, you average them (weighted by how much of their smoothing circle/sphere is inside of the box, if you want accuracy), and you multiply by your volume (1m^3). Since you are approximately conserving the density, your box will weight approximately 1000Kg.

Does that explanation help?

As for viscosity, that's a trickier one. The viscosity I used in the past was Monaghan's artificial viscosity, and it is called artificial for a reason (it has no physical meaning whatsoever). I found the proper values by tuning...

1 I had seen in some lecture notes that you actually add up the masses of the particles of SPH to find the mass inside a volume. That is correct (if not really accurate unless you take the smoothing volume into account) with some methods for initializing the particle's mass. However, with those methods what you don't get is a physical rest density value. Keep in mind that there is no problem with that. The simulation is going to be the same as long as the ratio between your rest density and mass is the same. You are actually free to choose how you want to proceed.

#### Attached Files

Edited by Javier Meseguer de Paz, 31 March 2014 - 03:20 PM.

“We should forget about small efficiencies, say about 97% of the time; premature optimization is the root of all evil” -  Donald E. Knuth, Structured Programming with go to Statements

"First you learn the value of abstraction, then you learn the cost of abstraction, then you're ready to engineer" - Ken Beck, Twitter

### #5Peacekeeper  Members

Posted 31 March 2014 - 07:01 PM

i really appreciate your detailed help. thanks alot for fixing the ecxel. i hope i dont annoy you with my frequent questions ;) i found my mistake in the excel sheet. obviously the functions where alright but my calcualation of the integrals wrong.

now its starting to flow! iam quite exited ^^

note:

- no viscosity yet - instead iam simply dampening the velocity by -0.1 in every timestep

- walls just move the particle back to their surface if crossed

- very low stiffnes of the fluid. so its compressible

ok i think the next step is a closer look at what you wrote about calculating the rest density tomorrow. right now i choosed working random values. sorry that iam a bit slow here. can you give an example for a 0.1m x 0.1m space with 100 particles and water ? your method seems to need an additional first loop for initialising the simulation?

edit: i think thats indeed what lets my simulation require such a tiny timestep (t=0.0003) or it explodes. as initially the particles are too close or to far away from each other  HUGE accelerations are calculated. after it slowed down i quess i could change the timestep to be quicker.

Edited by Peacekeeper, 31 March 2014 - 07:18 PM.

### #6Javier Meseguer de Paz  Members

Posted 01 April 2014 - 04:01 AM

Glad to be of help!

About the rest density. Ok, what your SPH solver will do is to enforce your rest density via pressure forces. Indeed, whatever your state equation, it will give you pressures that depend on your current density and your rest density.

Now, as said, SPH is unstable, and its unstable the higher the forces it generates (thus, the higher the timestep or the higher your current situation deviates from your target situation, in which the fluid is not compressed). So what you want, is to avoid compression as much as possible at all times. This begins with your initial state.

You want to initialize your simulation in a correct state. So you have two options: for a given mass and rest density, put the particles in the right positions (which can be hell to do) or given some particle positions and rest density, compute which has to be the mass of the particles. How to do this?

Before the simulation (we will discard the particles we are going to create):

- Create a particle (let's call it P) and set its mass to 1 (or whatever value, for that matter).

- Create its neighboring particles. Use whatever algorithm you are going to use in your real simulation. If for example you create particles every N meters in X and Y uniformly, just do that. The closer P resembles any particle inside your final fluid volume the better.

- Now, using the same algorithms and kernels you are going to use in your simulation, compute the density of P.

Now you know, for your setup, whats the ratio between your density and your particle mass. Now, do you want to have a particular rest density (eg. water's)? Easy! The mass you have to use is currentMass * watersDensity / currentDensity.

As for your timestep, there are formulas to compute it automatically. Basically you want a timesteps so that your particles don't move in one timestep more than a fraction of the smoothing length. Although it is a little more involved than that. I am sorry I don't recall the formula right now. I'll look for it if you are interested.

PS: You're not going slow, don't worry. SPH is a simple idea but it is filled with little details that make life... interesting.

“We should forget about small efficiencies, say about 97% of the time; premature optimization is the root of all evil” -  Donald E. Knuth, Structured Programming with go to Statements

"First you learn the value of abstraction, then you learn the cost of abstraction, then you're ready to engineer" - Ken Beck, Twitter

### #7Peacekeeper  Members

Posted 04 April 2014 - 11:07 AM

Hi!

i was rather busy the last days so my progress isnt that impressive, but i got your method to work today! works like charm.

its still rather slow due to python. in the future when everything is setup i might transfer everything to cython to speed it up. but thats not a problem right now. next comes the artificial viscosity.

i have a few general questions left here:

- all papers simulate water. but as sph is actually compressible what about an ideal gas? is that possible? actually ideal gas doesnt even have viscosity but probably that´d make the whole simulation quite unstable.

- i noticed that i use the spiky kernel for the pressure calculation AND the density calculation. as you can see above it works quite good. using the poly kernel for the density like its done in the literature results in clumping of the particles. yes, they clump together altough the spiky(-pressure) kernel should avoid that. i have no idea how this can happen.

Edited by Peacekeeper, 04 April 2014 - 11:08 AM.

### #8h4tt3n  Members

Posted 04 April 2014 - 11:18 AM

Hi there Peacekeeper. To avoid clumping, it might be a good idea to apply a repulsive-only kernel between each interacting particle pair, where force is proportional to either distance or a local pressure value ony dependent on those two particles. As far as I know (I'm no expert but have some experience with sph) it's the only way to avoid clumping. In my experience it works best if the particle pair kernel is spikier than the pressure kernel, which means that repulsive force is usually weaker than pressure force, except when particles are very near each other.

Cheers Mike

Edited by h4tt3n, 04 April 2014 - 11:19 AM.

### #9Peacekeeper  Members

Posted 04 April 2014 - 12:00 PM

yep and that should be the case here. as i can use the poly 6 kernel for calculating the density and pressure at each particle but for the pressure force (gradient) i use the spiky kernel. so if particles come close to each other there should still be a repulsive force

### #10h4tt3n  Members

Posted 04 April 2014 - 02:45 PM

Well okay, might just be a matter of fine-tuning then. What happens if you turn off pressure kernel, so the only forces influencing the particles are from the repulsive-only kernel?

### #11Peacekeeper  Members

Posted 04 April 2014 - 04:17 PM

oh sry maybe i misunderstood you. right now the particles move only based on the pressure gradient that is calculated using the spiky kernel. there is no extra kernel for repulsive forces when the particles are close because - as you can see - it works well if the spiky kernel is used on both (density/pressure and pressure forces). i just wonder why it works pretty well that way, as in many papers people tell you to use the poly 6 kernel for the density as it should result in more physically realistic behavior.

so i cant figure out why particles that get too close to each other do not repulse each other as the spiky kernel should provide repulsive forces. the question is why they even start to get so close in the first place. must have something todo with the kernel i do not understand. as if the particles have to get really close to reach the desired density.

edit: what about the ideal gas btw ?

Edited by Peacekeeper, 04 April 2014 - 04:18 PM.

### #12Javier Meseguer de Paz  Members

Posted 05 April 2014 - 05:45 AM

You are pretty much simulating a gas already (surprise!) but with the density of water and in the vacuum, so it ultimately behaves like water.

Actually very dense gases are not that different from liquids:

EDIT: Well, the "in the vacuum" part is pretty bad explained. In the end if we have a gas we are not in the vacuum, sorry about that. What I meant here is that we only have one phase (just our simulated fluid, without air surrounding it). What really matters here is that for dense gases the action of gravity is stronger than the repulsion forces of pressure until the particles are quite close together.

Edited by Javier Meseguer de Paz, 05 April 2014 - 02:59 PM.

“We should forget about small efficiencies, say about 97% of the time; premature optimization is the root of all evil” -  Donald E. Knuth, Structured Programming with go to Statements

"First you learn the value of abstraction, then you learn the cost of abstraction, then you're ready to engineer" - Ken Beck, Twitter

### #13h4tt3n  Members

Posted 06 April 2014 - 05:28 AM

oh sry maybe i misunderstood you. right now the particles move only based on the pressure gradient that is calculated using the spiky kernel. there is no extra kernel for repulsive forces when the particles are close because - as you can see - it works well if the spiky kernel is used on both (density/pressure and pressure forces). i just wonder why it works pretty well that way, as in many papers people tell you to use the poly 6 kernel for the density as it should result in more physically realistic behavior.

so i cant figure out why particles that get too close to each other do not repulse each other as the spiky kernel should provide repulsive forces. the question is why they even start to get so close in the first place. must have something todo with the kernel i do not understand. as if the particles have to get really close to reach the desired density.

edit: what about the ideal gas btw ?

Ah okay, slight misunderstanding. I compute a local density for each interacting particle pair and use that value for a repulsive-only "local" pressure force. Then I sum up all local densities to a "global" density value, which is used in an attractive-repulsive pressure kernel. This is a fail-safe way to 100% avoid clumping. There are many, many different kernels for computing density, pressure, and viscosity, and you are completely free to invent your own. My kernels are self-invented for working without knowing the distance between particles. This way I don't need to normalize vectors and hence not use square-root, which makes the simulation much faster.

In my experience, clumping happens when two particles are pushed very close together, and the kernels aren't capable to separate them; they make each others density kernels reach a value higher than rest-density, resulting in a high pressure. Then the pressure kernel pushes all other particles away from the two, but it doesn't manage to separate the particles causing the high density.

Edited by h4tt3n, 06 April 2014 - 07:40 AM.

### #14Peacekeeper  Members

Posted 06 April 2014 - 10:58 AM

ok thanks for the input i´ll stay with my working kernels for now.

do you have any further advice on vicsosity or rigid body collision? there are tony of possibilites aswell. which methods do you use? i´d prefer simple solutions. ;)

Old topic!

Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.