Fix accel func

[The9Cat]

Summary

This is a draft fix for the expansion center behavior described by @chervs in IntegrateOrbits()

Details

I put in all the infrastructure in pyEXP for propagating the expansion center to the AccelFunc basis evaluation but then didn't use it for providing translated coordinates for the force evaluation.

This PR is a fix for that.

@M1ssing-N0 : Would you compile this branch and try your failing example? Thanks.

[M1ssing-N0]

I wanted to test this code yesterday, but I had to fix a compatibility issue of my code with the latest version of pynbody...
This build works great in the test. I put identical particles in 3 potential settings, a stationary one, one moving with constant velocity, and one moving with constant acceleration. The trajectories are pretty consistent across these settings and the particle velocity matches the potential velocity.

I'm running a final test with a larger acceleration tomorrow, but I think the result will be similarly consistent.

[chervs]

Thanks @The9Cat!

Teng @M1ssing-N0, those tests worked already with the previous version, try to do the test where you put the Sgr potential at some fixed position (-X_sgr, 0,0) from the origin and integrate the orbit of a particle at rest at (0,0,0) see if it now moves towards Sgr in the negative x-direction...

[M1ssing-N0]

Also checked that one. The particle is pulled towards Sgr correctly.

I noticed that the particle may still fly away freely if the potential has a large acceleration. The acceleration is very large though, 1e5 km/s/Gyr. Is that a reasonable behavior for the code?

[The9Cat]

The code is simply using time-centered leap frog to solve the equations of motion given the force. So orbits should not be erratic unless the gravitational field itself is very rough Since you have a fixed potential, I recommend checking energy conservation for your orbit using the evaluation operator. E.g. if your basis instance is basis, then basis(x, y, z) or basis.getFields(x, y, z) will return a list of density, potential and force.

[michael-petersen]

I've exposed the center in pyEXP, such that along with the time of a coefficient structure, we can also inspect (and set!) the center. I was getting frustrated building an example that I couldn't quickly get the center.

On being readwrite: I could imagine the center being readonly as well: it's possible to cause some real problems with this. But by the same token, the time is also readwrite, and I've used that specifically, so maybe this is right.

I have a draft version that also duplicates the getCoefStruct syntax to be getCoefCenter, but in the end I wasn't sure that was needed, particularly if center is readwrite.

[The9Cat]

The initial design intentionally makes the center immutable. That is, the center is set when the coefficients are created and remains a fixed attribute. So explicitly readonly.

That may be too restrictive and I'm fine with a more flexible API. After all, one can modify coefficient values for various purposes.

Perhaps having setCoefCenter and getCoefCenter member function manipulators is a good compromise? As in, not allowing the center to be set as member data but allowing it to be set with an intentional call.

[michael-petersen]

I think that's a fine compromise -- mainly I want some way to get and set the centre in pyEXP.

1. Should .center be available at all? I've changed it to readonly, but the attribute could just disappear. Probably easier if it does disappear, to avoid confusion.
2. Should the name be more specific than getCenter and setCenter? I sort of feel like no since it is an already a member of a specific CoefStruct object, but I could be persuaded.

[The9Cat]

Looks good to me. On your two points:

1. I favor eliminating .center to avoid confusion. It's harmless enough as const member data but not sure we need or want both.
2. We could be more descriptive, as in getExpansionCenter and setExpansionCenter to identify the role of that quantity more explicitly. Those are long names, though. I think I suggested getCoefCenter etc. Shorter, but maybe that's too vague?

[M1ssing-N0]

The code is simply using time-centered leap frog to solve the equations of motion given the force. So orbits should not be erratic unless the gravitational field itself is very rough Since you have a fixed potential, I recommend checking energy conservation for your orbit using the evaluation operator. E.g. if your basis instance is basis, then basis(x, y, z) or basis.getFields(x, y, z) will return a list of density, potential and force.

The particle escaped the potential soon after the simulation starts (positive total energy) and the gravity from the potential become negligible, so I think the behavior is indeed reasonable.

I'm excited about the new member functions. The getter function is important for testing and the setter function can be particularly useful when working on a model with multiple components like the Milky Way.

Plus getCoefCenter sounds good to me, it's consistent with other functions so it's actually pretty easy to understand?

[michael-petersen]

Thanks, @M1ssing-N0 and @The9Cat for good naming suggestions. I've now

• Created getCoefTime and setCoefTime
• Changed the center names to getCoefCenter and setCoefCenter
• Removed the .center attribute

Last question: should .time still be accessible, and should it still be read/write? It's been that way since introduction as far as I can tell, so we might break some workflows if we change now. The path of least resistance is to not touch it, and then simply not recommend it to users (in favour of getCoefTime).

[The9Cat]

For convenience, I think it's fine to have .time (and .center for that matter) as readonly in the wrapper. They are redundant with the get/set members. But if they are readonly the user can't mistakenly change the values as a side effect of some workflow.

[michael-petersen]

Good idea. Okay, they are exposed readonly in the most recent commit!

I think if we're all happy we can go ahead and merge this.

[The9Cat]

These changes look fine to me. Do you have a MWE for pyEXP-examples? I haven't had a chance to work on this yet.

[michael-petersen]

I'm working on a MWE -- will have it set up soon and link here. The last commit sets an explicit zero center if the user doesn't pass anything (rather than being implicitly 0,0,0). This isn't ideal, I think, as the different structures could have a different number of dimensions, but having three zeros seems safest.

The center needs to be set for .setCoefCenter to work, as it checks the dimensionality. Another tactic could be to remove the guards from .setCoefCenter, and then no default initialisation is needed. Thoughts @The9Cat?

[michael-petersen]

There is now a draft MWE here, with the relevant notebook here.

[The9Cat]

Good catch. Thought we were doing that somewhere already, but no...

Yes, should be a 3d center. I think the design assumes a three-vector even if the system is 2d.

[Copilot]

Pull Request Overview

This PR fixes the acceleration function by propagating the expansion center to the force evaluation. Key changes include updating the Python bindings in pyEXP/CoefWrappers.cc, initializing and managing the center data in expui/CoefStruct.H, and applying the center offset in the field evaluations in expui/BiorthBasis.cc.

Reviewed Changes

Copilot reviewed 4 out of 5 changed files in this pull request and generated 3 comments.

File	Description
pyEXP/CoefWrappers.cc	Updated bindings: changed time to readonly and added center access and mutators with accompanying docstrings.
expui/CoefStruct.H	Initialized the center vector and added getter and setter functions for center and time.
expui/BiorthBasis.cc	Adjusted acceleration evaluation to subtract the expansion center from field coordinates.
expui/BasisFactory.H	Exposed the basis expansion center through a new getter.

Files not reviewed (1)

• CMakeLists.txt: Language not supported