Interpolation returns computed density rather than the interpolated density

[efaulhaber]

As discussed on Slack, this represents the computed density rather than the interpolated density. I would suggest using the interpolation approach described in the attached formula instead.

Originally posted by @LasNikas in #763 (comment)

[efaulhaber]

So, according to Bonet & Lok (1999), this is the default way to interpolate a quantity:

The corrected kernel interpolation (Shepard filter) is this:

This is exactly what @svchb implemented for our interpolation routines.

What is the interpolation strategy called that you mentioned above?

[LasNikas]

Note that they calculate the number density (SummationDensity).
So this is basically the same as in #808 when you replace rho_b with sum(m_b * W_b).
The equation I mentioned is the inverse distance weighting by Shepard.

[efaulhaber]

That doesn't make sense. Even if you calculate the summation density, this is not the same. Also, where does it say in the paper that this is the summation density?

[LasNikas]

In eq. (1) ?

[LasNikas]

I created a small example that demonstrates the effect of mass weighting on the result. In this case, only the mass is used for weighting, so the outcome is independent of the particle positions. You can see this by setting, e.g., perturbation_factor=0.4. This means that here, the mass acts purely as a weighting factor.
The setup below shows two domains. The top domain is with pressure=200 and the bottom domain with pressure=100. The interpolated pressure should be 150 exactly at the interface.
pressure_1 is at the red point and pressure_2 at the green point.

without perturbation:

┌ Info: Interpolated pressure_1 without mass weighting:
└   p_1 = 103.84664450727665
┌ Info: Interpolated pressure_1 with mass weighting:
└   p_2 = 150.0033166660959
┌ Info: Interpolated pressure_2 without mass weighting:
└   p_1 = 200.00000000000009
┌ Info: Interpolated pressure_2 with mass weighting:
└   p_2 = 199.99999999999991

with perturbation:

┌ Info: Interpolated pressure_1 without mass weighting:
└   p_1 = 104.2901137294761
┌ Info: Interpolated pressure_1 with mass weighting:
└   p_2 = 152.8436149322566
┌ Info: Interpolated pressure_2 without mass weighting:
└   p_1 = 199.99999999999994
┌ Info: Interpolated pressure_2 with mass weighting:
└   p_2 = 199.99999999999997

Note that the pressure deviates from 150 because the interface is no longer sharp due to the perturbation.

using TrixiParticles
using Plots
using LinearAlgebra

function interpolate_pressure(point, pressure, smoothing_length, domain_1, domain_2;
                              mass_weighting=true)
    sum_W_a_1 = 0.0
    m_1 = mass_weighting ? domain_1.mass[1] : 1

    for particle in 1:nparticles(domain_1)
        distance = norm(point - SVector{2, Float64}(domain_1.coordinates[:, particle]))

        if distance <=
           TrixiParticles.compact_support(SchoenbergQuinticSplineKernel{2}(),
                                          smoothing_length)
            W_a = TrixiParticles.kernel(SchoenbergQuinticSplineKernel{2}(), distance,
                                        smoothing_length)
            pressure += domain_1.pressure[particle] * W_a * m_1
            sum_W_a_1 += W_a * m_1
        end
    end

    sum_W_a_2 = 0.0
    m_2 = mass_weighting ? domain_2.mass[1] : 1

    for particle in 1:nparticles(domain_2)
        distance = norm(point - SVector{2, Float64}(domain_2.coordinates[:, particle]))
        if distance <=
           TrixiParticles.compact_support(SchoenbergQuinticSplineKernel{2}(),
                                          smoothing_length)
            W_a = TrixiParticles.kernel(SchoenbergQuinticSplineKernel{2}(), distance,
                                        smoothing_length)
            pressure += domain_2.pressure[particle] * W_a * m_2
            sum_W_a_2 += W_a * m_2
        end
    end

    pressure /= (sum_W_a_1 + sum_W_a_2)

    return pressure
end

fluid_density = 1000.0

particle_spacing_1 = 0.01
perturbation_factor = 0.0
domain_1 = RectangularShape(particle_spacing_1, (50, 50), (0.0, 0.0);
                            coordinates_perturbation=perturbation_factor,
                            density=fluid_density, pressure=100.0)

particle_spacing_2 = 0.05
domain_2 = RectangularShape(particle_spacing_2, (10, 10), (0.0, 0.5), density=fluid_density,
                            coordinates_perturbation=perturbation_factor,
                            pressure=200.0)

smoothing_length = 1.2 * particle_spacing_2

point_1 = SVector(0.25, 0.5)
point_2 = SVector(0.25, 0.7)

pressure_1 = 0.0
pressure_2 = 0.0

p_1 = interpolate_pressure(point_1, pressure_1, smoothing_length, domain_1, domain_2;
                           mass_weighting=false)
@info "Interpolated pressure_1 without mass weighting:" p_1

p_2 = interpolate_pressure(point_1, pressure_1, smoothing_length, domain_1, domain_2;
                           mass_weighting=true)
@info "Interpolated pressure_1 with mass weighting:" p_2

p_1 = interpolate_pressure(point_2, pressure_2, smoothing_length, domain_1, domain_2;
                           mass_weighting=false)
@info "Interpolated pressure_2 without mass weighting:" p_1

p_2 = interpolate_pressure(point_2, pressure_2, smoothing_length, domain_1, domain_2;
                           mass_weighting=true)
@info "Interpolated pressure_2 with mass weighting:" p_2

p1 = plot(domain_1, domain_2)
scatter!(p1, [point_1[1]], [point_1[2]], markersize=10, color=:red,
         label="point 1")
scatter!(p1, [point_2[1]], [point_2[2]], markersize=10, color=:green,
         label="point 2")

[svchb]

fixed by #808