Problem spawning many test particles at a location

[Spiruel]

I am attempting to launch test particles isotropically from the Earth. However sometimes when I add up to 1000 particles in a for loop, they spawn at the Sun's location and not the Earth. This is unusual behavior and I am unsure how best to solve this. My example code is below: (issue at for i in range(100): line)

import numpy as np

sim = rebound.Simulation()
sim.integrator = "janus"
sim.units = ('AU', 'yr', 'Msun')
solar_system_objects = ["Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"]
sim.add(solar_system_objects)
sim.N_active = 9
sim.testparticle_type = 1
sim.move_to_com()
sim.save("checkpoint.bin")
sim.status()

# collision and boundary options
sim.collision = "direct"
sim.collision_resolve = "merge"
sim.collision_resolve_keep_sorted = 1
#sim.boundary = "open"
#boxsize = 200.
#sim.configure_box(boxsize)
sim.track_energy_offset = 1

# set timestep
sim.dt = -0.001

# add then remove test particle to fix bug
sim.add(x=0, y=0, z=0, vx=0, vy=0.5, vz=0)
sim.remove(sim.N-1)

earth = sim.particles[3]
# velocities in AU/yr
AU_yr = 4740.57172
R = 0.00267383485
AU = 6.68459e-12

for i in range(100):
    phi = np.random.random() * np.pi
    theta = np.random.random() * 2 * np.pi
    x = R * np.sin(theta) * np.cos(phi)
    y = R * np.sin(theta) * np.sin(phi)
    z = R * np.cos(theta)

    dir_vec = np.array([x,y,z])
    rand = np.random.randint(20E3,40E3) / AU_yr
    vx, vy, vz = rand * (dir_vec/np.sqrt(dir_vec.dot(dir_vec)))
    sim.add(x=earth.x+x, y=earth.y+y, z=earth.z+z, vx=vx, vy=vy, vz=vz)

Noutputs = 1000
year = 2.*np.pi # One year in units where G=1
times = np.linspace(0.,-1E4*year, Noutputs)
part_num = sim.N
x = np.zeros((part_num,Noutputs))
y = np.zeros((part_num,Noutputs))

E0 = sim.calculate_energy()

sim.integrator = "janus" # IAS15 is the default integrator, so we actually don't need this line
sim.move_to_com()        # We always move to the center of momentum frame before an integration
ps = sim.particles       # ps is now an array of pointers and will change as the simulation runs

for i,time in enumerate(times):
    sim.integrate(time)
    for j in range(part_num):
        if ps[j].vx < 0.2:
            x[j][i] = ps[j].x   # This stores the data which allows us to plot it later
            y[j][i] = ps[j].y

dE = abs((sim.calculate_energy() - E0)/E0)
print dE

sim.save("rewound.bin")

[hannorein]

A couple of comments/questions, not sure if they solve your problem:

• Is there a specific reason that you are using JANUS as an integrator? When using JANUS, you need to either use the default units of G=1, or carefully set the scales in the JANUS struct. Because JANUS works on an integer grid, the scale of the problem actually matters. I would suggest you start out with another integrator unless you really need the time reversibility of JANUS.
• You have turned on collision detection, but none of your particles has a radius. Thus, there will never be any collision in your simulation.
• Why do you add and then remove a particle?

[Spiruel]

-The aim of this example code is to integrate backwards in time, and generate a population of objects that definitely impacts the Earth at a future point. If using Janus for this reason is the wrong decision please let me know and I'll switch the integrator. -I think collision detection was left in by mistake because I was worried that particles were falling into the Sun and being accelerated to impossible speeds without colliding/merging. Feel free to omit this. -When testing the environment in a Jupyter Notebook I found that sometime adding a test particle, then removing it, meant that the test particles correctly spawned on the Earth and not incorrectly spawning on the Sun. It was weird behaviour that I was struggling to make much sense of.

Many thanks for your help.

[hannorein]

Got it. So JANUS makes sense. The easiest way to fix the scale issue is to use units where G=1, unit of length = 1AU, unit of time = year/(2pi). Otherwise you need to adjust the scales in the JANUS struct accordingly. All kind of things could go wrong it you use scales that are off, so see if this makes a difference. I don't think that adding and removing a test particle does anything to resolve your problem. Also consider setting the mass of the test particles exactly equal to 0.

[Spiruel]

Thank you. I have ensured the default units G=1 are in place, and explicitly set the mass of test particles to 0. I have removed the references to collisions (as the particles have no radius).

Unfortunately, when testing this code in a Jupyter Notebook I am still getting particles spawned inside the Sun, unless I first make a bogus particle and then delete it. Are you also getting this behaviour?

[hannorein]

Could you just copy and paste the new code here. I want to make sure I run the same code as you...

[Spiruel]

Here you go:

import numpy as np

sim = rebound.Simulation()
sim.integrator = "janus"
solar_system_objects = ["Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"]
sim.add(solar_system_objects)
sim.N_active = 9
sim.testparticle_type = 1
sim.move_to_com()
sim.save("checkpoint.bin")
sim.status()

# set timestep
sim.dt = -0.001

# add then remove test particle to fix bug
sim.add(x=0, y=0, z=0, vx=0, vy=0.5, vz=0)
sim.remove(sim.N-1)

earth = sim.particles[3]
# velocities in AU/yr/2pi
AU_yr = 4740.57172
R = 0.00267383485
AU = 6.68459e-12

for i in range(1000):
    phi = np.random.random() * np.pi
    theta = np.random.random() * 2 * np.pi
    x = R * np.sin(theta) * np.cos(phi)
    y = R * np.sin(theta) * np.sin(phi)
    z = R * np.cos(theta)

    dir_vec = np.array([x,y,z])
    rand = np.random.randint(20E3,40E3) / AU_yr
    vx, vy, vz = rand * (dir_vec/np.sqrt(dir_vec.dot(dir_vec)))
    sim.add(m=0, x=earth.x+x, y=earth.y+y, z=earth.z+z, vx=vx, vy=vy, vz=vz)

Noutputs = 1000
year = 2.*np.pi # One year in units where G=1
times = np.linspace(0.,-1E4*year, Noutputs)
part_num = sim.N
x = np.zeros((part_num,Noutputs))
y = np.zeros((part_num,Noutputs))

E0 = sim.calculate_energy()

sim.move_to_com()        # We always move to the center of momentum frame before an integration
ps = sim.particles       # ps is now an array of pointers and will change as the simulation runs

for i,time in enumerate(times):
    sim.integrate(time)
    for j in range(part_num):
        if ps[j].vx < 0.2:
            x[j][i] = ps[j].x   # This stores the data which allows us to plot it later
            y[j][i] = ps[j].y

dE = abs((sim.calculate_energy() - E0)/E0)
print dE

sim.save("rewound.bin")

[hannorein]

Note that the units do not match up:

sim.units = ('AU', 'yr', 'Msun')

[hannorein]

I think a lot of the confusion is coming from the units. This code works for me as expected:

import numpy as np
import rebound

sim = rebound.Simulation()
sim.integrator = "janus"
solar_system_objects = ["Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"]
sim.add(solar_system_objects)
sim.N_active = 9
sim.testparticle_type = 1
sim.move_to_com()

sim.dt = -0.001
earth = sim.particles[3]
R = 0.00267383485
AU_yr = 4740.57172

for i in range(10):
    phi = np.random.random() * np.pi
    theta = np.random.random() * 2 * np.pi
    x = R * np.sin(theta) * np.cos(phi)
    y = R * np.sin(theta) * np.sin(phi)
    z = R * np.cos(theta)

    dir_vec = np.array([x,y,z])
    rand = np.random.randint(20E3,40E3) / AU_yr
    vx, vy, vz = rand * (dir_vec/np.sqrt(dir_vec.dot(dir_vec)))
    sim.add(m=0,x=earth.x+x, y=earth.y+y, z=earth.z+z, vx=vx, vy=vy, vz=vz)

sim.getWidget()

// in new cell:

sim.integrate(-20, exact_finish_time=0)

// backwards:

sim.dt *=-1
sim.integrate(0, exact_finish_time=0)

[Spiruel]

I think there was a lot of confusion for me with running this in a Jupyter notebook, where me changing the units was not being persistent. Your code is working well on my end. Thank you very much for your help Hanno.

[hannorein]

Great to heat it's working now. I'll close this issue but feel free to open a new one if a new problem comes up.

[dsspiegel]

FYI, for some issues along these lines, you might want to try nodebook http://multithreaded.stitchfix.com/blog/2017/07/26/nodebook/.

On Tue, Nov 14, 2017 at 9:06 AM, Hanno Rein notifications@github.com wrote:

Closed #311 https://github.com/hannorein/rebound/issues/311.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/hannorein/rebound/issues/311#event-1341320191, or mute the thread https://github.com/notifications/unsubscribe-auth/ABGUvhdgC2VWHuCqBcRjPcf1sInBkws8ks5s2cgxgaJpZM4Qbc3S .

-- Dave Spiegel, Ph.D.