Fix tafuni extrapolation

[LasNikas]

This PR fixes the extrapolation for BoundaryModelTafuni. Unfortunately, when I implemented it, I swapped the indices, because I didn't notices the kernel ( $W_{kl}$ ) and the position difference ( $\mathbf{x}_l - \mathbf{x}_k$ ) in the equation below are swapped.

It has also become clear to me now why Negi et al. (2020) modified the equation by introducing a sign change. However, I do not yet fully understand the reasoning behind the explanation for this modification.

Here are the screenshots from the code below.

On main

This PR

using TrixiParticles
particle_spacing = 0.05
domain_length = 1.0

function pressure_function(pos)
    t = 0
    U = 1.0
    b = -8pi^2 / 10
    x = pos[1]
    y = pos[2]

    return -U^2 * exp(2 * b * t) * (cos(4pi * x) + cos(4pi * y)) / 4
end

function velocity_function(pos)
    t = 0
    U = 1.0
    b = -8pi^2 / 10
    x = pos[1]
    y = pos[2]

    vel = U * exp(b * t) *
          [-cos(2pi * x) * sin(2pi * y), sin(2pi * x) * cos(2pi * y)]

    return SVector{2}(vel)
end

n_particles_xy = round(Int, domain_length / particle_spacing)

domain_fluid = RectangularShape(particle_spacing, (2, 1) .* n_particles_xy,
                                (0.0, 0.0), density=1000.0,
                                pressure=pressure_function,
                                velocity=velocity_function)

smoothing_length = 1.5 * particle_spacing
smoothing_kernel = WendlandC2Kernel{ndims(domain_fluid)}()
fluid_system = EntropicallyDampedSPHSystem(domain_fluid, smoothing_kernel,
                                           smoothing_length, 1.0)

fluid_system.cache.density .= domain_fluid.density

plane_out = ([2.0, 0.0], [2.0, domain_length])

outflow = BoundaryZone(; plane=plane_out, boundary_type=OutFlow(),
                       plane_normal=[-1.0, 0.0],
                       open_boundary_layers=10, density=1000.0, particle_spacing)
open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
                                          boundary_model=BoundaryModelTafuni(),
                                          buffer_size=0)

semi = Semidiscretization(fluid_system, open_boundary_out)
TrixiParticles.initialize_neighborhood_searches!(semi)

v_open_boundary = zero(outflow.initial_condition.velocity)
v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure')

TrixiParticles.set_zero!(open_boundary_out.pressure)

TrixiParticles.extrapolate_values!(open_boundary_out, v_open_boundary, v_fluid,
                                   outflow.initial_condition.coordinates,
                                   domain_fluid.coordinates, semi, 0.0;
                                   prescribed_pressure=false,
                                   prescribed_velocity=false)

trixi2vtk(fluid_system.initial_condition, filename="fluid",
          v=domain_fluid.velocity, p=domain_fluid.pressure)

trixi2vtk(open_boundary_out.initial_condition, filename="open_boundary_out",
          v=v_open_boundary, p=open_boundary_out.pressure)

plane_in = ([0.0, 0.0], [0.0, domain_length])

inflow = BoundaryZone(; plane=plane_in, boundary_type=InFlow(),
                       plane_normal=[1.0, 0.0],
                       open_boundary_layers=10, density=1000.0, particle_spacing)
open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
                                          boundary_model=BoundaryModelTafuni(),
                                          buffer_size=0)

semi = Semidiscretization(fluid_system, open_boundary_in)
TrixiParticles.initialize_neighborhood_searches!(semi)

v_open_boundary = zero(inflow.initial_condition.velocity)
v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure')

TrixiParticles.set_zero!(open_boundary_in.pressure)

TrixiParticles.extrapolate_values!(open_boundary_in, v_open_boundary, v_fluid,
                                   inflow.initial_condition.coordinates,
                                   domain_fluid.coordinates, semi, 0.0;
                                   prescribed_pressure=false,
                                   prescribed_velocity=false)

trixi2vtk(fluid_system.initial_condition, filename="fluid",
          v=domain_fluid.velocity, p=domain_fluid.pressure)

trixi2vtk(open_boundary_in.initial_condition, filename="open_boundary_in",
          v=v_open_boundary, p=open_boundary_in.pressure)

[codecov]

Codecov Report

Attention: Patch coverage is 78.57143% with 3 lines in your changes missing coverage. Please review.

Project coverage is 70.43%. Comparing base (3ed23c5) to head (0e594c0).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/io/write_vtk.jl	0.00%	2 Missing ⚠️
src/schemes/boundary/open_boundary/mirroring.jl	91.66%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #829      +/-   ##
==========================================
- Coverage   70.51%   70.43%   -0.08%     
==========================================
  Files         106      106              
  Lines        6865     6867       +2     
==========================================
- Hits         4841     4837       -4     
- Misses       2024     2030       +6     

Flag	Coverage Δ
unit	70.43% <78.57%> (-0.08%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[efaulhaber]

I think your example fails to highlight how incorrect the extrapolation on main is. Here is a simple velocity gradient on main:

This is the same with this PR:

[LasNikas]

/run-gpu-tests