GalacticDynamics-Oxford / Agama

Action-based galaxy modeling framework
Other
72 stars 36 forks source link

Wrong Potential for Gas and DarkMatter particles #11

Closed MartDonn closed 2 years ago

MartDonn commented 3 years ago

Hi! I'm using Agama with python 2 (anaconda/2/2019.03) to evaluate the vertical action of MW-like galaxies in the IllustrisTNG simulation. I'm having problems in evaluating the potentials of stars and gas particles.

The code runs but when I plot the circular velocity profile of gas and stars, it goes like crazy, and compared within the ones evaluated with the mass enclosed they are completely wrong.

Do you have any idea of what can be the problem?

eugvas commented 3 years ago

would be helpful to see the plot! or you can send me the problematic snapshot and your script. There are several places where things can go wrong, for instance, CylSpline is rather sensitive to the orientation of the coordinate system (should be aligned with the disk plane - if at all possible), and both Multipole and CylSpline can be sensitive to grid parameters (i.e. increasing the size of the grid doesn't help and can actually make things worse, amplifying discreteness noise). OTOH if you use Multipole with lmax=0 (i.e. a spherical system), I'd expect that the circular velocity curve should be pretty reliable. But the vertical action may be challenging because of fluctuations in the disk plane.

MartDonn commented 3 years ago

Here you find the script https://www.dropbox.com/s/h8lwe9ww9p52zx8/Action_AGAMA.py?dl=0

and here the plot: https://www.dropbox.com/s/xl5dxwmfwb71h3e/Vcirc.pdf?dl=0

I'm using CylSpline for stars and gas and Multipole for dark matter (which is the only curve that makes sense)

Thanks a lot!

MartDonn commented 3 years ago

Sorry, I forgot to mention that all particles are centered in the centre of mass of the galaxy and I have rotated X, Y, Z such that XY is face on and XZ is edge-on

eugvas commented 3 years ago

yeah, this certainly shouldn't be that bad! Nothing strikes me as terribly wrong in the script file, so I guess I need to look at the snapshots - can you also upload it somewhere? probably can keep only the essential data (positions and masses) by saving them into a text or nemo file:
agama.writeSnapshot('gas', (gas_coord, gas_mass), 'nemo')

MartDonn commented 3 years ago

Hi, any chance you can read the HDF5 file?

https://www.dropbox.com/s/as4y9hz6nrgfe1q/MW_342447_rotate.hdf5?dl=0

MartDonn commented 3 years ago

Sorry, this is the right folder: https://www.dropbox.com/sh/plhxu1obmtw0bnq/AADUEYPGNcadnn2QmoC3H3rAa?dl=0

eugvas commented 3 years ago

not sure what happens with your machine: when I run your script, I get a nice rotation curve: Vcirc of course, it's not without defects - the sudden jump at 30 kpc is due to the grid being too small: while stellar particles largely end at 30 kpc, gas particles are found up to 150-200 kpc, so you should increase Rmax and Zmax roughly to virial radius; this eliminates the discontinuity. Also I'd suggest to increase Rmin and Zmin a bit, probably to 0.5 kpc - this will avoid problems in case of slight centering defects (as happens in this snapshot - the stellar bar has a pronounced clump at ~0.2 kpc from the centre), will make the potential smoother and hence reduce the chance of action finder failures (it is quite sensitive to the noise). Don't worry too much about losing resolution: setting zmin comparable to the disc scale height is still fine, since the potential is interpolated cubically so the actual resolution is better than grid spacing. Also the disc plane varies by a comparable amount, so no point of hunting for smaller details. Could it be that your particle masses are somehow messed up? Can you check pot_stars.totalMass() vs np.sum(stars_mass)? Can you also try replacing CylSpline by Multipole for both stars and gas to see if that removes the glitch? I'm not suggesting this as a permanent fix - the potential of these highly flattened (in case of gas disc) components is much better represented by CylSpline; with Multipole you would need to go to high order in lmax, but this also amplifies high-frequency noise.

MartDonn commented 3 years ago

So indeed if I print: pot_stars.totalMass() = 5.05074441997e+12 np.sum(stars_mass) = 9.4613987328e+10 pot_gas.totalMass() = 1.04300736786e+12 np.sum(gas_mass) = 9.2946579456e+10 pot_dm.totalMass() = 9.28344600029e+11 np.sum(dm_mass) = 92.83446e+10

So something happens with mstars and mgas. BUT, I have also tried to use Multipole for all the three potentials and it works perfectly, with a reasonable rotation curve.

Do you think could be a problem with Python 2?

eugvas commented 3 years ago

no it shouldn't be different between Py2 and Py3. But you can try, just in case. Are you using a recent version? (there haven't been any significant changes in CylSpline for a while, though)

eugvas commented 3 years ago

did you manage to check if the problem persists with another version of the code or Python?