Discrepancy in SPH Density Summation Results Using OpenFPM

Question

I'm working on a Smoothed Particle Hydrodynamics (SPH) simulation using the OpenFPM framework. I've set up a 2D non-dimensionalized Poiseuille flow simulation and am observing a minor but consistent discrepancy in the density summation results.

Problem Description:

Particle spacing is dp = 4, leading to a particle volume of V = dp^2 = 16. Ideally, the SPH summation should sum up to 1/16 = 0.06250. However, my results are slightly off:

kernel_summation : 0.062503951537163940366248482405354

Kernel Implementation: I'm using the quintic spline kernel:

inline double densitySummation(const double r_)
    const double ALPHA_2D = 7.0 / (478.0 * M_PI);
    double q = r_ / H;
    double result;
    if (q <= 1.0)
        result = std::pow((3.0 - q),5) - 6.0 * std::pow((2.0 - q),5) + 15.0 * std::pow((1.0 - q),5);
    else if (q <= 2.0)
        result = std::pow((3.0 - q),5) - 6.0 * std::pow((2.0 - q),5);
    else if (q <= 3.0)
        result = std::pow((3.0 - q),5);
    else
        result = 0.0;
    return ALPHA_2D * result / (H * H);

This is the function which I am using to validate the density summation. I am setting up the domain and iterating over it, and doing this summation, without any shifting, nevertheless, the summation keeps having the error.

inline void densitySummationNoCorrection(particles &vd, myCelllist &NN)
{
    // particle iterator
    auto part = vd.getDomainIterator();

    // For each particle ...
    while (part.isNext())
    {
        auto a = part.get();
        if (vd.template getProp<type>(a) == FLUID_PARTICLES)
        {
            Point<2, double> xa = vd.getPos(a);
            double massa = (vd.getProp<type>(a) == FLUID_PARTICLES) ? MassFluid : MassBound;
            auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(a)));

            double kernel
            double summationInLoop = 0.0;
            while (Np.isNext())
            {
                auto b = Np.get();
                double massb = (vd.getProp<type>(a) == FLUID_PARTICLES) ? MassFluid : MassBound;
                double rhob = vd.template getProp<rho>(a);

                // Get particle properties b
                Point<2, double> xb = vd.getPos(b);
                Point<2, double> dr = xa - xb;
                Point<2, double> DW;
                double r_ = sqrt(norm2(dr));
               

                double wab = Wab(r_);
                kernel += wab
                

                ++Np;
            }

            vd.template getProp<Volume>(a) = kernel;
            vd.template getProp<rho_densitySummation>(a) = massa*kernel;
           
        }
        ++part;
    }
}

Constants:

constexpr double a2_ = 7.0 / 478.0 * oneOverPi * one_over_H * one_over_H;
constexpr double a2_grad = -5.0 * a2_ * one_over_H;

Validation: I've validated the simulation parameters using another simulation suite where they worked as expected. The domain setup and particle spacing appear correct, with a precise spacing of 4.000.

I'm puzzled by this minor discrepancy. While it's small, it's consistent and I'm keen to understand its origin. Any insights or suggestions would be greatly appreciated.

Thank you in advance!

Best, Ivan

as described above