Add a bar potential to MWPotential2014 in AMUSE

[stevenalfonso]

Hello,
I am trying to add a bar to the MWPotential2014 as pot = MWPotential2014+SoftenedNeedleBarPotential(), and then convert it to AMUSE mwp_amuse=to_amuse(pot, reverse=True). However, when added to the AMUSE integration, there is an error due to potential shape, I guess.

Below is the code:

from amuse.lab import *
from amuse.couple import bridge
from amuse.datamodel import Particles
from galpy.potential import to_amuse, SoftenedNeedleBarPotential, MWPotential2014
from galpy.util import plot as galpy_plot
from amuse.community.bhtree.interface import BHTree
from amuse.units import quantities
from tqdm import tqdm

Convert galpy MWPotential2014 to AMUSE representation

pot = MWPotential2014 + SoftenedNeedleBarPotential()
mwp_amuse= to_amuse(pot, reverse=True)

Set initial cluster parameters

N= 1000
Mcluster= 1000. | units.MSun
Rcluster= 10. | units.parsec
Rinit = [x_init.value, y_init.value, z_init.value] | units.pc
Vinit = [-vx_init.value, -vy_init.value, -vz_init.value] | units.km/units.s # flip velocities

Setup star cluster simulation

def setup_cluster(N,Mcluster,Rcluster,Rinit,Vinit):
converter= nbody_system.nbody_to_si(Mcluster,Rcluster)
stars= new_plummer_sphere(N,converter)
stars.x+= Rinit[0]
stars.y+= Rinit[1]
stars.z+= Rinit[2]
stars.vx+= Vinit[0]
stars.vy+= Vinit[1]
stars.vz+= Vinit[2]
return stars, converter

tend = 120.0 | units.Myr
dt = 1.0 | units.Myr
times = quantities.arange(0.0 | units.Myr, tend, dt)

Setup cluster

stars,converter= setup_cluster(N,Mcluster,Rcluster,Rinit,Vinit)
cluster_code= BHTree(converter,number_of_workers=1) #Change number of workers depending no. of CPUs
cluster_code.parameters.epsilon_squared= (3. | units.parsec)**2
cluster_code.parameters.opening_angle= 0.6
cluster_code.parameters.timestep= dt
cluster_code.particles.add_particles(stars)

Setup channels between stars particle dataset and the cluster code

channel_from_stars_to_cluster_code= stars.new_channel_to(cluster_code.particles,
attributes=["mass", "x", "y", "z", "vx", "vy", "vz"])
channel_from_cluster_code_to_stars= cluster_code.particles.new_channel_to(stars,
attributes=["mass", "x", "y", "z", "vx", "vy", "vz"])

Setup gravity bridge

gravity= bridge.Bridge(use_threading=False)

Stars in cluster_code depend on gravity from external potential mwp_amuse (i.e., MWPotential2014)

gravity.add_system(cluster_code, (mwp_amuse,))

External potential mwp_amuse still needs to be added to system so it evolves with time

gravity.add_system(mwp_amuse,)

Set how often to update external potential

gravity.timestep= cluster_code.parameters.timestep/2.

Evolve

time= 0.0 | tend.unit
for i in tqdm(times):
gravity.evolve_model(i)

channel_from_cluster_code_to_stars.copy()
gravity.stop()

galpy_plot.plot(stars.x.value_in(units.kpc),stars.y.value_in(units.kpc),'.', xlabel=r'$X,(\mathrm{kpc})$',ylabel=r'$Y,(\mathrm{kpc})$')

The error comes when evolving the gravity code, this is the fuller error message:

---

ValueError Traceback (most recent call last)
time= 0.0 | tend.unit
for i in tqdm(times):
--> gravity.evolve_model(i)
channel_from_cluster_code_to_stars.copy()
gravity.stop()

File ~/.local/lib/python3.8/site-packages/amuse/couple/bridge.py:555, in Bridge.evolve_model(self, tend, timestep)
timestep = self.timestep
if self.method is None:
--> return self.evolve_joined_leapfrog(tend,timestep)
else:
return self.evolve_simple_steps(tend,timestep)

File ~/.local/lib/python3.8/site-packages/amuse/couple/bridge.py:570, in Bridge.evolve_joined_leapfrog(self, tend, timestep)
while self.time < (tend-timestep/2.):
if first:
--> self.kick_codes(timestep/2.)
first=False
else:

File ~/.local/lib/python3.8/site-packages/amuse/couple/bridge.py:711, in Bridge.kick_codes(self, dt)
for x in self.codes:
if hasattr(x,"kick"):
--> de += x.kick(dt)
self.kick_energy += de

File ~/.local/lib/python3.8/site-packages/amuse/couple/bridge.py:454, in GravityCodeInField.kick(self, dt)
if(self.verbose):
print(self.code.class.name,"receives kick from",field_code.class.name, end=' ')
--> self.kick_with_field_code(
particles,
field_code,
dt
if(self.verbose):
print(".. done")

File ~/.local/lib/python3.8/site-packages/amuse/couple/bridge.py:490, in GravityCodeInField.kick_with_field_code(self, particles, field_code, dt)
def kick_with_field_code(self, particles, field_code, dt):
--> ax,ay,az=field_code.get_gravity_at_point(
self._softening_lengths(particles),
particles.x,
particles.y,
particles.z
self.update_velocities(particles, dt, ax, ay, az)

File ~/.local/lib/python3.8/site-packages/galpy/potential/amuse.py:159, in galpy_profile.get_gravity_at_point(self, eps, x, y, z)
phi = numpy.arctan2(y.value_in(units.kpc), x.value_in(units.kpc))

Cylindrical force

--> Rforce = potential.evaluateRforces(
self.pot,
R / self.ro,
zed / self.ro,
phi=phi,
t=self.tgalpy,
use_physical=False,
)
phitorque = potential.evaluatephitorques(
self.pot,
R / self.ro,
(...)
use_physical=False,
) / (R / self.ro)
zforce = potential.evaluatezforces(
self.pot,
R / self.ro,
(...)
use_physical=False,
)

File ~/.local/lib/python3.8/site-packages/galpy/potential/Potential.py:82, in potential_positional_arg..wrapper(Pot, *args, **kwargs)
@wraps(func)
def wrapper(Pot, /, *args, **kwargs):
--> return func(Pot, *args, **kwargs)

File ~/.local/lib/python3.8/site-packages/galpy/util/conversion.py:1113, in potential_physical_input..wrapper(*args, **kwargs)
if (
"zmax" in kwargs
and _APY_LOADED
and isinstance(kwargs["zmax"], units.Quantity)
kwargs["zmax"] = kwargs["zmax"].to(units.kpc).value / ro
--> return method(*args, **kwargs)

File ~/.local/lib/python3.8/site-packages/galpy/util/conversion.py:1008, in physical_conversion..wrapper..wrapped(*args, **kwargs)
if use_physical_explicitly_set:
warnings.warn(
"Returning output(s) in internal units even though use_physical=True, because ro and/or vo not set"
)
return method(*args, **kwargs)

File ~/.local/lib/python3.8/site-packages/galpy/potential/Potential.py:2168, in evaluateRforces(Pot, R, z, phi, t, v)
@potential_positional_arg
@potential_physical_input
@physical_conversion("force", pop=True)
def evaluateRforces(Pot, R, z, phi=None, t=0.0, v=None):
"""
Evaluate the radial force F_R(R,z,phi,t) of a potential, force or a list of potentials/forces.

(...)
"""
--> return _evaluateRforces(Pot, R, z, phi=phi, t=t, v=v)

File ~/.local/lib/python3.8/site-packages/galpy/potential/Potential.py:2190, in _evaluateRforces(Pot, R, z, phi, t, v)
out += pot._Rforce_nodecorator(R, z, phi=phi, t=t, v=v)
else:
--> out += pot._Rforce_nodecorator(R, z, phi=phi, t=t)
return out
elif isinstance(Pot, Potential):

File ~/.local/lib/python3.8/site-packages/galpy/potential/Potential.py:204, in Potential._Rforce_nodecorator(self, R, z, phi, t)
def _Rforce_nodecorator(self, R, z, phi=0.0, t=0.0):
# Separate, so it can be used during orbit integration
try:
--> return self._amp * self._Rforce(R, z, phi=phi, t=t)
except AttributeError: # pragma: no cover
raise PotentialError(
"'_Rforce' function not implemented for this potential"
)

File ~/.local/lib/python3.8/site-packages/galpy/potential/SoftenedNeedleBarPotential.py:98, in SoftenedNeedleBarPotential._Rforce(self, R, z, phi, t)
def _Rforce(self, R, z, phi=0.0, t=0.0):
--> self._compute_xyzforces(R, z, phi, t)
return numpy.cos(phi) * self._cached_Fx + numpy.sin(phi) * self._cached_Fy

File ~/.local/lib/python3.8/site-packages/galpy/potential/SoftenedNeedleBarPotential.py:126, in SoftenedNeedleBarPotential._compute_xyzforces(self, R, z, phi, t)
def _compute_xyzforces(self, R, z, phi, t):
# Compute all rectangular forces
--> new_hash = hashlib.md5(numpy.array([R, phi, z, t])).hexdigest()
if new_hash != self._force_hash:
x, y, z = self._compute_xyz(R, phi, z, t)

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (4,) + inhomogeneous part.

How can I solve this? Or is something that I am doing wrong?

Thanks!

[jobovy]

Thanks for the detailed example! I think the issue is that AMUSE expects the potential to be able to be evaluated at multiple points simultaneously, while SoftenedNeedleBarPotential doesn't support that. @webbjj wrote the AMUSE interface and has used it, so he might be able to comment. For how, you could try adding a DehnenBarPotential instead, which should work (we used it in https://ui.adsabs.harvard.edu/abs/2020MNRAS.496.5101P/abstract), but note that that one is just a perturbation, so can just be added to MWPotential2014 (no need to remove bulge mass to make room for the bar).

[webbjj]

@jobovy is correct. The Galactic potential is evaluated at each star's position, so arrays of each star's coordinates is sent to the galpy potential for evaluation. If SoftenedNeedleBarPotential doesn't support that, then it won't work.

If you have strong reasons for preferring the SoftenedNeedleBarPotential over the DehnenBarPotential that Jo recommends, the only work around I can think of would be to integrate the cluster's orbit and the cluster's internal dynamical evolution separately. However in this case tidal stripping would not be modelled accurately and you would have to provide your own criteria for a star being stripped from the cluster. So it would depend on your science goals.