mikegrudic / pytreegrav

Fast N-body gravitational force and potential in Python
MIT License
63 stars 9 forks source link

Strange divergence in potential estimate #18

Open ttepperg opened 3 weeks ago

ttepperg commented 3 weeks ago

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 commented 3 weeks ago

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 commented 3 weeks ago

Hi Mike,

Thanks so much for your prompt reply; that was really unexpected =)

I’m digging a bit deeper to find which particles exactly are causing this; but no problem, I’ll be happy to share the data with you asap!

I’ll get back to you soon,

Thor


Dr. Thor Tepper García | Astrophysicist Sydney Institute for Astronomy (SIfA) The University of Sydney www.thorsten.mxhttp://www.thorsten.mx

The University of Sydney A globally ranked universityhttps://www.sydney.edu.au/about-us/our-world-rankings.html Top 20 university in the world | 1st for sustainability in Australia *QS World University Rankings 2024

From: Mike Grudić @.> Date: Monday, 26 August 2024 at 15:25 To: mikegrudic/pytreegrav @.> Cc: Thorsten Tepper Garcia @.>, Author @.> Subject: Re: [mikegrudic/pytreegrav] Strange divergence in potential estimate (Issue #18)

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?

— Reply to this email directly, view it on GitHubhttps://github.com/mikegrudic/pytreegrav/issues/18#issuecomment-2309345134, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AK3O47VHOL4WIFUS6WRAWF3ZTK347AVCNFSM6AAAAABNDIRVZGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBZGM2DKMJTGQ. You are receiving this because you authored the thread.Message ID: @.***>

ttepperg commented 3 weeks ago

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);
    1. 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 commented 3 weeks ago

Thanks, this is very helpful. Looking into it.

mikegrudic commented 3 weeks ago

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 commented 3 weeks ago

I have added warnings to inform the user of this issue where present: https://github.com/mikegrudic/pytreegrav/commit/b4b0663d8667e0167178222f95f6be93c07bff48

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 commented 3 weeks ago

Mike,

Thank you for looking into the issue. and for the very clear explanation of what is going on! Also thanks for the code update.

I’m afraid I’m not savvy enough about the tree algo to modify the code, and contribute a PR, but I’ll certainly keep using / testing it, and report bugs if I happen to come across some.

Cheers

Thor


Dr. Thor Tepper García | Astrophysicist Sydney Institute for Astronomy (SIfA) The University of Sydney www.thorsten.mxhttp://www.thorsten.mx

The University of Sydney A globally ranked universityhttps://www.sydney.edu.au/about-us/our-world-rankings.html Top 20 university in the world | 1st for sustainability in Australia *QS World University Rankings 2024

From: Mike Grudić @.> Date: Tuesday, 27 August 2024 at 04:36 To: mikegrudic/pytreegrav @.> Cc: Thorsten Tepper Garcia @.>, Author @.> Subject: Re: [mikegrudic/pytreegrav] Strange divergence in potential estimate (Issue #18)

I have added warnings to inform the user of this issue where present: b4b0663https://github.com/mikegrudic/pytreegrav/commit/b4b0663d8667e0167178222f95f6be93c07bff48

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.

— Reply to this email directly, view it on GitHubhttps://github.com/mikegrudic/pytreegrav/issues/18#issuecomment-2310827443, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AK3O47R7X3D2QK7QDBUW5EDZTNYTTAVCNFSM6AAAAABNDIRVZGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJQHAZDONBUGM. You are receiving this because you authored the thread.Message ID: @.***>