[JOSS Review] Sign of acceleration components?

[santisoler]

While playing with the minimal Python example shown in the README.md I defined a cube in the (+,+,+) octant, so one of its vertices is in [1, 1, 1] and its opposite one in [2, 2 ,2], with a density of 1.0. Then I defined a computation point in the origin ([0, 0, 0]) and used the evaluate() function to compute its fields.

import numpy as np
import polyhedral_gravity

cube_vertices = np.array(
    [
        [1, 1, 1],
        [2, 1, 1],
        [2, 2, 1],
        [1, 2, 1],
        [1, 1, 2],
        [2, 1, 2],
        [2, 2, 2],
        [1, 2, 2],
    ]
)
cube_faces = np.array(
    [
        [1, 3, 2],
        [0, 3, 1],
        [0, 1, 5],
        [0, 5, 4],
        [0, 7, 3],
        [0, 4, 7],
        [1, 2, 6],
        [1, 6, 5],
        [2, 3, 6],
        [3, 7, 6],
        [4, 5, 6],
        [4, 6, 7],
    ]
)
cube_density = 1.0
computation_point = np.array([0, 0, 0])

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)
print(acceleration)

After running this code I get:

[-5.715199089414289e-12, -5.715199089414289e-12, -5.7151990894142846e-12]

The three components of the acceleration are negative, which would imply that any particle located in the computation point would feel an acceleration towards the (-,-,-) octant, meaning it would feel a repulsive force from the gravitational interaction with the cube. Given that the cube has a positive density, I'd assume that the three components should be positive, not negative.

Is this the expected result? Or are there any sign conventions that I missed?

---

This issue is part of the JOSS review: openjournals/joss-reviews#6384

[JoachimSchwabe]

After checking again, I think you are right. In context of my issue #29 it appears that both the signs of all quantities are flipped in general (at least for my test body and according to the sign convention for the potential in my field), and additionally the Vz/Vy/Vz are flipped again with respect to the potential and the tensor.

According to the definition in my field (physical geodesy), the potential outside a body with positive density is also defined positive. The gravitational acceleration is then the NEGATIVE derivative of the potential, e.g. dg_z = -Vz. Since the potential becomes smaller at further distance of the body, the sign of Vx/Vy/Vz should then be negative, whereas the acceleration is positive (but becomes also smaller in magnitude with increasing distance).

UPDATE 2024-03-11:

1. I checked once more with the example cube from the description. My computation points were located at (0, 0, 3) and (0, 0, 3.001) above the test cube, and the signs were as follows:
   V : positive (decreasing upwards) -> Vz > 0
   acc[2] : positive acceleration (decreasing upwards) --> negative of Vz !
   tensor[2] : positive --> correct sign of Vzz = d/dz(Vz), but flipped w.r.t. to output of acc[2]

2. From this I conclude that the mesh of my own test body is indeed "inverted", so that all quantities are additionally flipped.

[schuhmaj]

The problem

So, the primary concern is that the acceleration's sign is inverted compared to the potential and the second derivative tensor.
Hence, the potential fix is the simple multiplication with -1.
I totally agree with that, especially given that, as @JoachimSchwabe reported, the numerical difference of the potential yields a value with an inverted sign to the returned acceleration.

Handling until now

However, I want to put the following things for reference in the issue (since this sign question still is a bit unclear to me since I am not a domain expert):

• Tsoulis et al.'s reference implementation does not mention how to handle the sign. He returns the accelerations in his Fortran and Matlab Implementations as absolute values.
• The analytical cube gravity, which we used as a reference, defines the derivative of the cube as the negative derivative of the potential. However, the potential is defined as negative in the first place. (So, this is aligned to our acceleration but inverted to our potential - One of my advisors meant at that time that this is just a "definition thing.")
• We utilized this polyhedral gravity model as the ground truth for neural density fields. Here, potential and acceleration were both multiplied with -1. This yielded the same results as the Mascon Model utilized in a previous study iteration.

Towards a solution

To summarize, I am highly convinced by the fact that the numerical derivatives sign is flipped when compared to the derivative as returned by the implementation.

So, I would definitely add a minus multiplication to the acceleration output.
@santisoler @JoachimSchwabe I would still ask you for your confirmation that this is the change you desire?

---

Anyways, many thanks the two of you for bringing up this issue in detail and stating examples - this is highly appreciated.

[santisoler]

Hi @schuhmaj. Thanks for the detailed reply. I think the solution of adding a -1 to the acceleration will solve the issue I get in that particular example, although I think it would be more important to be sure why the implementation needs that -1 sign, so the fix applies to every possible case and not just this one.

Regarding your comments on why this is needed, I'm not very familiar with Tsoulis et al.'s implementation, and I'm not sure why they return absolute values. But I could say that your advisor is right when they stated that the issue the sign of the potential is a "definition thing".

Lets consider the gravity acceleration vector $\mathbf{g}(\mathbf{p})$ on an observation point $\mathbf{p}$ due to a particle located on $\mathbf{q}$ with mass $m$. If $m>0$, then this vector will always be attractive, therefore pointing in the $\mathbf{q} - \mathbf{p}$ direction. Since gravitational forces are conservative, we can define a potential $V$ such as:

$$ \mathbf{g} = - \nabla V $$

The minus sign used in the definition of $V$ is a convention. Mostly in Physics, potentials are defined with the minus sign, so negative values of the potential mean attractive fields, while positive values mean repulsive fields. This is why in most physical sciences (including the analytical cube gravity article you shared), the gravitational potential is always negative for $m>0$.

Nonetheless, in geodesy and geophysics, the positive sign convention is usually used, so $V$ is defined as (Hofmann-Wellenhof and Moritz, 2005, eq. 1-9):

$$ \mathbf{g} = + \nabla V $$

Meaning that gravitational potentials have positive values for $m>0$ bodies. This convention doesn't impact the sign of the acceleration vector components (which are physical observable quantities that don't depend on sign conventions) nor on the signs of the tensor components.

In conclusion, the sign of the potential could be either positive or negative, as long as it is clear what convention the implementation is applying. Regarding the acceleration vector, its sign must match the attractive nature of gravitational interactions for positive masses.

---

Hofmann-Wellenhof and Moritz (2005). Physical Geodesy. doi: 10.1007/b139113

[JoachimSchwabe]

I think I have to make a general comment here:

Firstly, as I now realized, with the correct mesh orientation of the test cube the package delivers a positive potential (outside the test body), decreasing when facing away from the mass. As outlined below by @santisoler, this is in line with the convention in physical geodesy.

Now, regarding this geodetic convention:

1. Indeed in physical geodesy, the potential of a mass with positive density is always considered positive. It is understood as the work to be done against the gravitational force to bring a test mass away from the body to infinity.
2. Yes, the actual physical (measurable) acceleration (i.e., length of the potential gradient vector | grad V |) is positive and facing towards the mass.
3. However, for the partial derivatives Vx/Vy/Vz of the potential ("acceleration"), we have to consider the location of the computation point with respect to the mass in the coordinate system. In gravimetry, the observation/computation point is usually above the mass (z_P > 0 and z_mass < z_P, z is upwards). This means that Vz of the positive potential is negative in the computation point (V decreasing upwards). On the other hand, what is measured with gravimeters is the actual acceleration acting downwards. This is why in gravimetry, per definition acc_z = -Vz! The tensor element Vzz = d/dz(Vz) = -d/dz(acc_z) is then again positive, since the negative Vz becomes smaller upwards (away from the mass). If we now assume the test cube hovering in the air and the computation point is below the cube (z_P < 0, z_mass > z_P), consequently Vz must become positive along the z axis (towards the cube). But because of the convention acc_z = -Vz, the code will output a negative acc_z in this case (although the actual physical attraction is also positive towards the z axis in this case, because according to the definition above "acceleration" is defined as "downwards"). Vzz is again positive, because the positive Vz becomes larger along the z axis (towards the mass).
4. The same also holds for all the other coordinate axis. The "acceleration" is defined as the negative of Vx/Vy/Vz, which both can be positive or negative depending whether the respective axis is facing towards or away from the mass! This makes sense because if we consider any body in a free xyz coordinate system, in fact there is no "up", "left" or "right". Only that the implementation of the package (right-handed coordinate system) happens to match with geodetic coordinates such as e.g. UTM coordinates: "East" (x), "North" (y), Height (z).

To conclude: Sign of the accelerations here are defined always "downwards" (irrespective of where the mass, the computation point and the coordinate system actually are) for z, or towards "west"/negative x and "north"/negative y, accordingly. Therefore, I can only strongly advise not to just take the absolute value of any of the vector components. This should be clearly stated in the description of the package. In particular, the table in paragraph "Output" is wrong and should be corrected to -Vx, -Vy and -Vz.

[schuhmaj]

So, if I get this correctly, sorry, if I am asking so stupidly, but I am quite a lot confused with the sign:

---

At the moment:

• My acceleration is NOT facing towards the mass, but inverted (100% sure about that given the geodesynet study in which we inverted the acceleration so that it actually points towards the mass)
• Hence, I am returning: V, (-Vx, -Vy, -Vz), (Vxx, Vyy, Vzz, Vxy, Vxz, Vyz)?

---

So, to comply with the geodesy and geophysics conventions, I must:

• TODO 1a) Change the documentation, saying I am returning -Vx, -Vy, -Vz and not as it is currently: Vx, Vy, Vz

OR

• TODO 1b) I return the accelerations multiplied with -1.0 --> So that the returned acclereation points towards the mass

AND ANYWAYS

• TODO 2) Specify in the docs that the implementation follows the geodesy and geophysics conventions

---

Is this correct?

[santisoler]

So, if I get this correctly, sorry, if I am asking so stupidly, but I am quite a lot confused with the sign:

There are not stupid questions 🙂. I think is super valuable that you are asking for clarification.

I would suggest to take the 1b path. I think it's easier to understand and document the definition of the potential as positive and that the acceleration vector is returned as:

$$ \mathbf{g} = \nabla V = \left( \frac{\partial V}{\partial x}, \frac{\partial V}{\partial y}, \frac{\partial V}{\partial z} \right) $$

If the 1a path is chosen, then it would be wrong to call the (-Vx, -Vy, -Vz) vector the acceleration, because it's not any acceleration vector that could be observed in such scenario.

---

Minor comment regarding orientation of axis. I appreciate the point that @JoachimSchwabe makes regarding the gravity acceleration measurements in geophysics. Nonetheless, I think the main confusion that comes from (and this are historical reasons) the fact that in most geophysical literature, authors refer to $g_z$ as the downward component of the acceleration vector. By doing that mapping between $g_z$ and the downward direction, they were implicitly setting the $z$ axis pointing downwards. This imposed a problem when we need to generalize the forward models on systems where we don't have a natural up or down direction (asteroids for example), or even when our masses are above the observation point.

In conclusion, given the fact that polyhedral-gravity-model is mainly aimed to model bodies like asteroids, where there might not be a natural choice for the coordinate system, I suggest to follow a general formulation and keep definitions simpler.

From personal preference, I would use the geodetic convention of defining the potential $V$ as always positive for positive masses, and return the acceleration as the gradient of $V$ with respect to whatever x, y and z coordinates the user is defining (as long as the observation point and the mesh are defined in the same coordinate system). In summary, it would be the 1b path (as I stated before).

[JoachimSchwabe]

I pretty much agree with santosoler. Again, whether Vx etc. or -Vx etc. points towards the mass depends on the coordinate system and the location of mass and computation point. And as I had written before and santosoler pointed out again, the convention with the negative sign for the acceleration comes from the case of surface gravimetry where the acceleration is considered as the vertical component only and points downwards, whereas the height coordinate points upwards. This concept does not really make sense for general 3d coordinates, or also for the "horizontal" component of the gravity vector, if x and y are associated with east and north.
Thus, I would suggest to output positive V and Vx, Vy, Vz and tensor accordingly to prevent confusion and make a clear statement about it. It is the task of the user to e.g. take the absolute value of either component or the length of the gradient vector for the positive acceleration towards the mass, if that is needed, depending on the mass body and the coordinate system.

[schuhmaj]

Thanks a lot to both of you for taking the time to detail your explanation a second time. Now, I get it 🙂

#31 fixes the sign issues by adding a -1 multiplication to the acceleration while also improving the documentation with regard to the convention being utilized.