Question in the introduction example

[BBBBBbruce]

Hi, I am trying to use pmesh to build a nbody simulation, but I cannot get the correct force out of the computation.

I try to use the example in the introduction. Since that example is not directly runnable, I made a few changes.

I think P['Accel'] is the gravitational force of each particle, so I compare this force with the force computed by brute force which loop through each particle and sumup the all the gravitational force. And the two forces I got are very different.

One problem I notice is when I increase Nmesh, the force from pmesh would get smaller in magnitude.

I know that particle mesh should be an approxiamtion, so I shouldnt expect two forces to be equal. but I think that the error level should be reasonable. and the level should be decreasing when I increase Nmesh.

Sorry if my question is fundamental. But I couldnt find any tutorial except the exmple in the introduction. Can you see what I am missing?

my code is:

G = 1 
mass = 1
BoxSize = 10 
Nparticles = 5
Nmesh = 200
np.random.seed(42)  # For reproducibility

P = {'Position': np.random.uniform(0, BoxSize, (Nparticles, 3)),
     'Accel': np.zeros((Nparticles, 3))}

def compute_gravitational_acceleration(Pin,Nmesh):
    P = copy.deepcopy(Pin)
    pm = ParticleMesh([Nmesh, Nmesh, Nmesh],BoxSize=BoxSize, dtype='f8', comm= MPI.COMM_WORLD)
    def potential_transfer_function(k, v):
        k2 = sum(ki ** 2 for ki in k)
        return -4 * np.pi * G * v / (k2 + (k2 == 0))

    layout = pm.decompose(P['Position'])
    density = pm.create(type='real')
    density.pm.paint(P['Position'], mass = mass, layout=layout,out = density)
    pot_k = density.r2c().apply(potential_transfer_function)

    for d in range(3):
        def force_transfer_function(k, v, d=d):
            ki = k[d]
            return 1j * ki * v

        force_d = pot_k.apply(force_transfer_function).c2r()
        P['Accel'][:, d] = force_d.readout(P['Position'], layout=layout)
    return P['Accel']

for brute force:

def compute_gravitational_bruteforce(Pin):
    P = copy.deepcopy(Pin)
    for i in range(Nparticles):
        accel = np.zeros(3)
        for j in range(Nparticles):
            if i != j:
                diff = P['Position'][j] - P['Position'][i]
                r = np.sqrt(np.sum(diff**2) )
                accel += G * mass * diff / r**3
        P['Accel'][i] = accel
    return P['Accel']