NBody Simulation questions
Moderators: gmalivuk, Moderators General, Prelates
NBody Simulation questions
I'm trying to model an Nbody system, and in addition to gravity, I want to have a force which models the repulsive forces of stellar wind, star formation, and other phenomenon. Hopefully something of the form F(m1, m2, distance).
I tried using PV = nRT to solve for T, then using StefanBoltzmann, but that gave me power, not force, and I started to think using the ideal gas law wasn't exactly a good idea, given how stars deviate from ideal gasses.
Can anyone come up with a simple approximation for this force?
I tried using PV = nRT to solve for T, then using StefanBoltzmann, but that gave me power, not force, and I started to think using the ideal gas law wasn't exactly a good idea, given how stars deviate from ideal gasses.
Can anyone come up with a simple approximation for this force?
Last edited by >) on Sat Jul 18, 2015 8:36 pm UTC, edited 1 time in total.
Re: Equation for repulsive forces in Nbody simulation?
Drag forces (which would include stellar wind forces) tend to go like the velocity at low speeds (Reynolds number much less than 1) and the square at high speeds (Reynolds number much higher than 1). You can probably look up the radiusspeed relation online or, if not, I can try looking it up in the morning (you get it by modelling the flow as a pipe of crosssection 4pi*r^2).
Using Stefanboltzman you can get the rate if massloss and from all of this get a repulsive force of the form F(m_{1},distance)
Edit: I looked up an actually usable version of the formula
In the case of an isothermal interstellar medium stationary at infinity:
u^{2} = 2 c_{s}^{2} ln(rho_{infinity}/rho) + 2GM/r
Where c_{s} is the speed of sound in the surrounding medium.
This version was derived for spherical accretion (hence the assumption of it being stationary at infinity), but accretion and stellar winds are just timereversals of each other and everything applies the same. There are refiniements you can do but, seeing as your forces are just proportion u or u^{2}, most of the interesting differences can be bumped off into the coefficient of proportionality.
More general forms are a pain, have unphysical solutions (due to the approximations necessary in their derivation) and are complicated to relate back to physical quantities but, if you want you, you can use (normalising velocities by c_{s} and distances by r_{s} = GM/2c_{s}^{2})
u^{2}  2ln(u) = 4 ln(r) + 4/r + C
For some constant C related to the mass flux. C=3 is the only value with sonic transitions.
Using Stefanboltzman you can get the rate if massloss and from all of this get a repulsive force of the form F(m_{1},distance)
Edit: I looked up an actually usable version of the formula
In the case of an isothermal interstellar medium stationary at infinity:
u^{2} = 2 c_{s}^{2} ln(rho_{infinity}/rho) + 2GM/r
Where c_{s} is the speed of sound in the surrounding medium.
This version was derived for spherical accretion (hence the assumption of it being stationary at infinity), but accretion and stellar winds are just timereversals of each other and everything applies the same. There are refiniements you can do but, seeing as your forces are just proportion u or u^{2}, most of the interesting differences can be bumped off into the coefficient of proportionality.
More general forms are a pain, have unphysical solutions (due to the approximations necessary in their derivation) and are complicated to relate back to physical quantities but, if you want you, you can use (normalising velocities by c_{s} and distances by r_{s} = GM/2c_{s}^{2})
u^{2}  2ln(u) = 4 ln(r) + 4/r + C
For some constant C related to the mass flux. C=3 is the only value with sonic transitions.
my pronouns are they
Magnanimous wrote:(fuck the macrons)
Re: Equation for repulsive forces in Nbody simulation?
Thank you. Just to be clear, u is the speed of the stellar wind?
Re: Equation for repulsive forces in Nbody simulation?
Yup.
Strictly speaking the forces will be a bit more complicated because the viscosity of the interstellar medium will lead to interactions between the various winds but hopefully it should be a reasonable assumption to treat the vairous winds as independent.
You'll also need to pick the coefficient of proportionality for the drag force, and that choice should be informed by a physical example if possible.
Strictly speaking the forces will be a bit more complicated because the viscosity of the interstellar medium will lead to interactions between the various winds but hopefully it should be a reasonable assumption to treat the vairous winds as independent.
You'll also need to pick the coefficient of proportionality for the drag force, and that choice should be informed by a physical example if possible.
my pronouns are they
Magnanimous wrote:(fuck the macrons)

 Posts: 6
 Joined: Wed Dec 29, 2010 11:24 pm UTC
Re: Equation for repulsive forces in Nbody simulation?
Writing simulations like this is basically what I'm doing my phd on. If you want something realistic, a pairwise force depending only on position isn't really going to work. I see two relatively easy paths you can take.
The simplest follows directly from the track you were on with P = nkT. If you're willing to take your stars and other temperature sources as static, you can bake a static pressure field onto the simulation space (either analytically or using a grid). The gradient of the pressure field will give you a force, which you can apply to every particle. So long as the effect of the particles you are simulating on the background isn't significant, this is the way I'd go.
If you actually want to simulate the dynamics of a significant mass of the fluid, I'd go with SPH (Smoothed Particle Hydrodynamics). There are plenty of other fluid simulation methods, but SPH is pretty easy to codeespecially since you're already using a particle method to simulate gravity already.
What context are you planning on using this in?
The simplest follows directly from the track you were on with P = nkT. If you're willing to take your stars and other temperature sources as static, you can bake a static pressure field onto the simulation space (either analytically or using a grid). The gradient of the pressure field will give you a force, which you can apply to every particle. So long as the effect of the particles you are simulating on the background isn't significant, this is the way I'd go.
If you actually want to simulate the dynamics of a significant mass of the fluid, I'd go with SPH (Smoothed Particle Hydrodynamics). There are plenty of other fluid simulation methods, but SPH is pretty easy to codeespecially since you're already using a particle method to simulate gravity already.
What context are you planning on using this in?
Re: Equation for repulsive forces in Nbody simulation?
Just for fun, I'm trying to simulate the large scale structure of the universe by placing a bunch of particles and hoping they collapse into filaments and clusters and such. I wanted a repulsive force both to prevent everything from clumping together and to prevent force from going to infinity as radius goes to zero, so stray particles don't go flying off at ridiculous speeds. I was actually hoping that the force would be proportional to distance^k, k < 2, so that the force would dominate at small distances. I'll see about SPH

 Posts: 6
 Joined: Wed Dec 29, 2010 11:24 pm UTC
Re: Equation for repulsive forces in Nbody simulation?
If you just want to avoid the singularity in the 1/r^{2} force, there's actually a very easy answer! Most NBody sims use some form of softened gravity to avoid this. The simplest is Plummer softening, where you replace 1/r^{2} with 1/(r^{2} + \epsilon^{2}), where \epsilon (the "forcesoftening length") is set by the scale of the structures you want to resolve. For cosmology, you can probably get away with around ~200 kpc if this is just for fun. This will resolve massive halos and larger structures nicely.
Re: Equation for repulsive forces in Nbody simulation?
thanks!
By the way, a bit offtopic, but would you happen to have any idea on what would be a good initial starting condition? I've read that a lot of large scale simulations use actual cosmological data, but if I just wanted to place the particles, how should I? I was thinking the mass of the particles should either be uniformly distributed or follow some sort of power law. I was going to start the position of the particles off almost uniformly distributed, but I'm not sure if velocities should be zero too.
By the way, a bit offtopic, but would you happen to have any idea on what would be a good initial starting condition? I've read that a lot of large scale simulations use actual cosmological data, but if I just wanted to place the particles, how should I? I was thinking the mass of the particles should either be uniformly distributed or follow some sort of power law. I was going to start the position of the particles off almost uniformly distributed, but I'm not sure if velocities should be zero too.

 Posts: 6
 Joined: Wed Dec 29, 2010 11:24 pm UTC
Re: Equation for repulsive forces in Nbody simulation?
The standard starting point in the field is to use the Zel'dovich approximation with the observationally determined density fluctuation powerspectrum. This is much less intimidating than it might sound.
You're going to want to start with a uniform 3d grid of uniform mass particles. Each particle on the grid is then going to be perturbed off of it's initial position by a small displacement. The velocity for the particle will be proportional to the displacement. These displacements will form a Gaussian random fielda distribution that takes a powerspectrum as a parameter. If you look around online, there are a lot of resources that will go through generating this in a variety of programming languages.
A really simple starting point would be to make the displacement vector a small random number. This effectively using a constant powerspectrum (white noise), which is physically pretty unrealistic, but doesn't involve any additional math. Set the initial velocities to zero. This isn't consistent, but the error it introduces is relatively small, and decays with time.
You're going to want to start with a uniform 3d grid of uniform mass particles. Each particle on the grid is then going to be perturbed off of it's initial position by a small displacement. The velocity for the particle will be proportional to the displacement. These displacements will form a Gaussian random fielda distribution that takes a powerspectrum as a parameter. If you look around online, there are a lot of resources that will go through generating this in a variety of programming languages.
A really simple starting point would be to make the displacement vector a small random number. This effectively using a constant powerspectrum (white noise), which is physically pretty unrealistic, but doesn't involve any additional math. Set the initial velocities to zero. This isn't consistent, but the error it introduces is relatively small, and decays with time.
Re: NBody Simulation questions
I tried it out, and it seems to work ok. But I'm still not seeing filaments. The uniform grid seems to collapse, forming some geometric patterns (probably because of the barneshut quadtree). Then, large areas of high and low density temporarily appear before the simulation reverts to a more chaotic, uniform state. This gif should illustrate what I mean  I simulated 4096 bodies over 31.7 billion years on a toroidal grid with sidelength ~32 Mpc.
I'm not sure when I should expect to see filaments or structures  during the collapse or after it? I suspect that my timestep might be too high, but I'm not sure when I should slow it down if I want to see the detail. Or perhaps I need more particles?
If I'm supposed to be looking before the grid collapses completely, the geometric patterns are slightly offputting. But I have a feeling I'm not supposed to be looking after the grid collapses, because it takes around 25 billion years just for everything to "stabilize".
I was also thinking I might need to take into account cosmological inflation by modifying all the distances by an increasing multiplier. Equivalently, I could just weaken gravity over time. I'm unsure if this is significant though.
I'm not sure when I should expect to see filaments or structures  during the collapse or after it? I suspect that my timestep might be too high, but I'm not sure when I should slow it down if I want to see the detail. Or perhaps I need more particles?
If I'm supposed to be looking before the grid collapses completely, the geometric patterns are slightly offputting. But I have a feeling I'm not supposed to be looking after the grid collapses, because it takes around 25 billion years just for everything to "stabilize".
I was also thinking I might need to take into account cosmological inflation by modifying all the distances by an increasing multiplier. Equivalently, I could just weaken gravity over time. I'm unsure if this is significant though.
 Xanthir
 My HERO!!!
 Posts: 5426
 Joined: Tue Feb 20, 2007 12:49 am UTC
 Location: The Googleplex
 Contact:
Re: NBody Simulation questions
Those clean rectangular collapse regions at the beginning show that something's clearly off, for one.
But also, iiuc the largescale structure is a result of quantum fluctuations causing random density changes immediately before hyperinflation, which turned them from "quantum scale" to "cosmological scale". Gravity could then work on the differences and collapse the slightlydenser regions into filaments.
But also, iiuc the largescale structure is a result of quantum fluctuations causing random density changes immediately before hyperinflation, which turned them from "quantum scale" to "cosmological scale". Gravity could then work on the differences and collapse the slightlydenser regions into filaments.
(defun fibs (n &optional (a 1) (b 1)) (take n (unfold '+ a b)))
Re: NBody Simulation questions
The grid looks like an aliasing effect from the spacing of the bodies and the size of the torus causing some fourier mode to be excited. Changing the size of the torus by just a few pixels should drastically alter the number of grid squares (if I'm right).
The fact everything starts off stationary might also be a reason you're not seeing the sorts of things you're expecting (and getting regular structure instead) if you give each dot a small random velocity it'd hopefully help disrupt the grid further and might get stuff more like what you want.
The fact everything starts off stationary might also be a reason you're not seeing the sorts of things you're expecting (and getting regular structure instead) if you give each dot a small random velocity it'd hopefully help disrupt the grid further and might get stuff more like what you want.
my pronouns are they
Magnanimous wrote:(fuck the macrons)
Re: NBody Simulation questions
Actually it turned out to be a pretty severe bug in my Barnes Hut implementation which I've fixed.
Here's the corrected simulation, covering just 15 billion years
Here's frame four hundred. There's definitely clusters and voids, but I'm not sure about the filaments. With a bit of creative license, I can point out a few possible filaments, but I'm not sure if they're dense enough to qualify as such.
This image covers about a quarter the size as my simulation, but the filaments look much more welldefined. In addition to the specks of light, the image also seems to contain long, stretched out strings of light, and I can't tell if that's making the filaments look more welldefined than they actually are.
Basically I'm not sure if my filaments are as prominent as they should be.
Here's the corrected simulation, covering just 15 billion years
Here's frame four hundred. There's definitely clusters and voids, but I'm not sure about the filaments. With a bit of creative license, I can point out a few possible filaments, but I'm not sure if they're dense enough to qualify as such.
This image covers about a quarter the size as my simulation, but the filaments look much more welldefined. In addition to the specks of light, the image also seems to contain long, stretched out strings of light, and I can't tell if that's making the filaments look more welldefined than they actually are.
Basically I'm not sure if my filaments are as prominent as they should be.
Re: NBody Simulation questions
It might be that your attractive and repulsive forces are balanced the same as in real life or of the right form. Alternatively, I could see it being due to your grid being too coarse/your universe too small (I am just guessing on both of these by the way).
my pronouns are they
Magnanimous wrote:(fuck the macrons)
Re: NBody Simulation questions
Originally, I wanted to comment on probable fact that you aren't accounting for relativity (both gravity delay and curvature), but the wiki page states it isn't significant anyway.
Do you have to account for "particles" exerting different amounts of radiation? I can imagine some particles represent mostly dark matter and some represent mostly normal matter. It's unlikely that this leads to more filaments, but it's worth a try I guess.
Do you have to account for "particles" exerting different amounts of radiation? I can imagine some particles represent mostly dark matter and some represent mostly normal matter. It's unlikely that this leads to more filaments, but it's worth a try I guess.
Re: NBody Simulation questions
I was actually trying to find the equation for the repulsive/radiation force, but ended up only implementing plummer softening. Since radiation drops off proportional to distance squared, I reasoned modeling it would be equivalent to weakening the gravitational constant a little. Plus, as icicletasty mentioned, I'd need some SPH to get that working. Essentially, I have everything acting as dark matter right now. I'm actually not sure if the Millenium Simulation simulated any baryonic matter, or whether it has significant effect though.
I think my universe should be somewhere near the right size, as the wikipedia image I linked before is a quarter the size of my universe and it has filaments too. I do plan on trying with more particles, but even 4096 is a bit taxing for my computer.
I'm also not actually sure why the collapse of the original nearuniform state of the matter should form filaments at all, rather than just clusters. What force is holding the filaments together? Gravity should collapse the filaments into clusters
I think my universe should be somewhere near the right size, as the wikipedia image I linked before is a quarter the size of my universe and it has filaments too. I do plan on trying with more particles, but even 4096 is a bit taxing for my computer.
I'm also not actually sure why the collapse of the original nearuniform state of the matter should form filaments at all, rather than just clusters. What force is holding the filaments together? Gravity should collapse the filaments into clusters

 Posts: 6
 Joined: Wed Dec 29, 2010 11:24 pm UTC
Re: NBody Simulation questions
4096 particles is probably on the lower end of where you're going to see any real large scale structure. You can also try a more realistic power spectrum for the initial conditions. As Xanthir pointed out, in real cosmology this comes from the inital quantum fluctuations that are blown up to large scales by inflation/the big bang (with an additional feature on top of them of cosmological origin, the Baryon Acoustic Oscillations).
As for why filaments form, it's because objects collapse exponentially faster allong their shorter dimensions, forming first pancakes, then filaments. The filaments will then gradually feed into clusters, so eventually they will disappear.
Looking at your gifs, I worry that your time steps may be too large. What integration method are you using and how are you choosing your timesteps?
As for why filaments form, it's because objects collapse exponentially faster allong their shorter dimensions, forming first pancakes, then filaments. The filaments will then gradually feed into clusters, so eventually they will disappear.
Looking at your gifs, I worry that your time steps may be too large. What integration method are you using and how are you choosing your timesteps?
Re: NBody Simulation questions
If it's necessary I can try a run with 16384 or 65536 particles.
I looked up the zel'dovich approximation and the math seems to be a bit beyond me. By power spectrum do you mean the distribution of the magnitudes of the displacements? Can it be described as a normal/poisson/gamma distribution? Currently I'm just displacing each point <= 0.2 times the lattice spacing on each of the axes. I've also set random initial velocity to be <= 1E5 m/s for every object, which is significantly less than the Local Group's velocity of ~600 km/s, just to be safe.
I'm using RK4 to integrate, with timestep = 1E15 seconds because that allows the simulation to run through the lifetime of the universe in one night's worth of sleep . Is there any canonical method to pick timesteps?
edit: This paper suggests setting dt so that average displacement < epsilon/2, where epsilon is the Plummer softening length.
Also, considering that Hubble's Law (70 km/s/Mpc * 32 Mpc * 5E17s = 1.12E24 m) means my simulation more than doubles in "size" over the time period, do you think scaling up the dimensions of the universe to simulate metric expansion is worthwhile?
I looked up the zel'dovich approximation and the math seems to be a bit beyond me. By power spectrum do you mean the distribution of the magnitudes of the displacements? Can it be described as a normal/poisson/gamma distribution? Currently I'm just displacing each point <= 0.2 times the lattice spacing on each of the axes. I've also set random initial velocity to be <= 1E5 m/s for every object, which is significantly less than the Local Group's velocity of ~600 km/s, just to be safe.
I'm using RK4 to integrate, with timestep = 1E15 seconds because that allows the simulation to run through the lifetime of the universe in one night's worth of sleep . Is there any canonical method to pick timesteps?
edit: This paper suggests setting dt so that average displacement < epsilon/2, where epsilon is the Plummer softening length.
Also, considering that Hubble's Law (70 km/s/Mpc * 32 Mpc * 5E17s = 1.12E24 m) means my simulation more than doubles in "size" over the time period, do you think scaling up the dimensions of the universe to simulate metric expansion is worthwhile?
Re: NBody Simulation questions
One problem with Runge–Kutta methods is that they don't conserve energy, and that could be killing your filaments. Try using a symplectic integrator like Verlet, or Leapfrog. I generally use the synchronized version of Leapfrog (for doing simple orbit sims) but the standard version's probably adequate for your purposes.
As the Wiki quote mentions, the Hamiltonian that Leapfrog conserves isn't exactly equal to the true Hamiltonian; you can reduce the perturbation by making the time step smaller. But you might not need to worry about that too much  I suspect that a perturbed system with energy conservation will be very similar in structure to the true system, since it won't suffer the same "meltdown" problem that you're getting with RK4.
Wikipedia wrote:One use of this equation is in gravity simulations, since in that case the acceleration depends only on the positions of the gravitating masses (and not on their velocities), although higher order integrators (such as Runge–Kutta methods) are more frequently used.
There are two primary strengths to Leapfrog integration when applied to mechanics problems. The first is the timereversibility of the Leapfrog method. One can integrate forward n steps, and then reverse the direction of integration and integrate backwards n steps to arrive at the same starting position.
The second strength of Leapfrog integration is its symplectic nature, which implies that it conserves the (slightly modified) energy of dynamical systems. This is especially useful when computing orbital dynamics, as many other integration schemes, such as the (order 4) RungeKutta method, do not conserve energy and allow the system to drift substantially over time.
As the Wiki quote mentions, the Hamiltonian that Leapfrog conserves isn't exactly equal to the true Hamiltonian; you can reduce the perturbation by making the time step smaller. But you might not need to worry about that too much  I suspect that a perturbed system with energy conservation will be very similar in structure to the true system, since it won't suffer the same "meltdown" problem that you're getting with RK4.
Re: NBody Simulation questions
I implemented metric expansion and also leapfrog integration. For some reason, the simulations with leapfrog integration seem to collapse much slowly.When I add metric expansion (linear expansion from 5E23 m to 1E24 m), barely any clusters have time to form. However, when I computed the error in the change in positions between leapfrog and rk4, the majority of the particles were off by less than 10%, with no bias either way, so it's not obviously an implementation error. I didn't think it was possibly to have a bug in something as simple as leapfrog, but I'm not sure this is how physics is supposed to work either.
RK4 with metric expansion: gif, last frame
RK4 with no metric expansion: gif, last frame
timesteps are 1E15 seconds
Leapfrog with metric expansion: gif, last frame
Leapfrog with no metric expansion: gif, last frame
timesteps are 2.5E14 seconds. note that the gifs run four times as slow because there are four times as many frames.
Anyway at this point I think I'll see what 16384 bodies gets me with a lower timestep
RK4 with metric expansion: gif, last frame
RK4 with no metric expansion: gif, last frame
timesteps are 1E15 seconds
Leapfrog with metric expansion: gif, last frame
Leapfrog with no metric expansion: gif, last frame
timesteps are 2.5E14 seconds. note that the gifs run four times as slow because there are four times as many frames.
Anyway at this point I think I'll see what 16384 bodies gets me with a lower timestep
Who is online
Users browsing this forum: No registered users and 12 guests