Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Strange divergence in potential estimate #18

Open
ttepperg opened this issue Aug 26, 2024 · 7 comments
Open

Strange divergence in potential estimate #18

ttepperg opened this issue Aug 26, 2024 · 7 comments

Comments

@ttepperg
Copy link

Hi Mike!

First of all, so many thanks for making pytreegrav publicly available; it's simply AWESOME!

At the moment, I'm using it to estimate the potential induced by the particles in an N-body simulation,
and in general I find very reasonable results. However, I've noticed that there are some points at which
the estimated tree potential differs by orders of magnitude (!) from the 'true' potential. The latter is
obtained directly from the simulation output.
Attached below is a plot showing the behaviour of the 'true' and pytreegrav potential with distance for this
particular simulation snapshot. Ideally, the orange and blue curves should be on top of one another.
That they diverge at large distance is expected I guess, because the tree code’s accuracy diminishes with distance.
Bu the 'spikes' close to 0 are what is worrying me. I've no idea where they come from. These spikes are caused by
a small number of locations (particles) with very large potential values.

I'm estimating the potential using:

    pytreegrav_pot, pytreegrav_octtree = \
        Potential(pos, mass, G=Grav,
        theta=0.7,
        method='tree',
        return_tree=True, parallel=True)

where pos and mass are the position and mass arrays.

ANY help to resolve this issue would be highly appreciated!

Cheers

Thor

pytreegrav

@mikegrudic
Copy link
Owner

oh yikes. Have not seen something like this. Is it possible for you to share the positions and masses so i can try to reproduce it?

@ttepperg
Copy link
Author

ttepperg commented Aug 26, 2024 via email

@ttepperg
Copy link
Author

Hi Mike,

While digging deeper I came across another weird problem: it seems as is there is some stochasticity (!) in the calculation of the potential.

I’m attaching a script, a data file (with positions and masses), and a plot. These serve two purposes:

  1. reproduce a case where there are ‘spikes’ in the otherwise smooth potential for no real apparent reason (the particles masss and positions that cause these look OK to me);
  2. demonstrate the somewhat stochastic nature of the calculation for the same data set, which is REALLy weird to me.

The plot shows the resulting potential ‘profile’ for each iteration, each of which merely consists in loading the same
data set and calculate the potential, over and over again.

Please note that you’ll need to change the script extensions (our firewall does not allow to send anything that may look
like an executable of any kind).

Maybe I’m doing something very wrong, but I truly cannot understand the reason for any of these
two (perhaps related?) issues; I hope you can!

I’m running this on a Mac Catalina, Python 3.9, Pytreegrav 1.1.2

Please let me know if I can assist in any way to find out what is going on.

Cheers

Thor
pytreegrav_iter
test_pytreegrav.txt
pos_mass.zip

@mikegrudic
Copy link
Owner

Thanks, this is very helpful. Looking into it.

@mikegrudic
Copy link
Owner

OK so the standout feature of your point set is that your coordinates are not unique. The source of the stochasticity is that the tree-build procedure deals with this case by perturbing positions randomly on the order of machine precision.

The actual solution for the potential is undefined at the positions of the overlapping points, because the potential diverges at 0. What pytreegrav does right now certainly doesn't satisfy the principle of least astonishment - the most reasonable behaviour would be to return -inf for the undefined points. At the very least, I think a check for uniqueness and a warning is warranted.

For your problem, you will need either unique positions or a finite softening length to get a well-defined answer.

@mikegrudic
Copy link
Owner

I have added warnings to inform the user of this issue where present: b4b0663

I will leave this issue open because there is a less-weird way to handle this: don't perturb the positions and just assign to negative infinity. Just requires some rejiggering of the tree-build and -walks. I might get around to this, but anyone else should feel free to contribute a PR.

@ttepperg
Copy link
Author

ttepperg commented Aug 26, 2024 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants
@mikegrudic @ttepperg and others