pysph icon indicating copy to clipboard operation
pysph copied to clipboard

Wrong coefficient for 1D-QuinticSpline

Open Jasmine969 opened this issue 6 months ago • 1 comments

Image

The normalization coefficient for 1D-QuinticSpline is $1/(60h)$ rather than $1/(120h)$.

Validation:

def quintic_spline(r, cut, dimension=3):
    h = cut / 3
    q = r / h
    w = np.where(q < 3, (3 - q) ** 5, 0)
    w -= np.where(q < 2, 6 * (2 - q) ** 5, 0)
    w += np.where(q < 1, 15 * (1 - q) ** 5, 0)
    if dimension == 1:
        alphaD = 1 / (60 * h)
    elif dimension == 2:
        alphaD = 7 / (478 * np.pi * h ** 2)
    else:  # dimension == 3
        alphaD = 1 / (120 * np.pi * h ** 3)
    return w * alphaD

kernels = {'Quintic-Spline': quintic_spline}

def test_kernel(name, cutoff):
    from scipy.integrate import quad, dblquad, tplquad
    res1d, _ = quad(lambda r: kernels[name](r, cut=cutoff, dimension=1),
                    0, cutoff)
    res2d, _ = dblquad(lambda r, theta: kernels[name](r, cut=cutoff, dimension=2) * r,
                       0, 2 * np.pi, 0, cutoff)
    res3d, _ = tplquad(lambda r, phi, theta: kernels[name](r, cut=cutoff, dimension=3) * r * r * np.sin(phi),
                       0, 2 * np.pi, 0, np.pi, 0, cutoff)
    flag_pass = True
    for i, res in enumerate([res1d, res2d, res3d]):
        print(f'Dimension {i+1}: result={res}.')
        if abs(res - 1) > 1e-6:
            print(f'Dimension {i + 1} does not satisfy the normalization condition!')
            flag_pass = False
    if flag_pass:
        print(f'Kernel {name} has passed the normalization condition test!')

if __name__ == '__main__':
    test_kernel('Quintic-Spline', 3)

Output:

Dimension 1: result=0.9999999999999204. Dimension 2: result=0.9999999999972096. Dimension 3: result=0.9999999999634721. Kernel Quintic-Spline has passed the normalization condition test!

Jasmine969 avatar Jul 30 '25 07:07 Jasmine969

@Jasmine969 -- Many thanks for this report and the code! Will fix this later this week in the main branch. It looks like none of our tests are checking this correctly so a proper test case will also need to be added.

prabhuramachandran avatar Aug 03 '25 15:08 prabhuramachandran