Fix OBC extrapolation (again)

[LasNikas]

This PR fixes the bug with the swapped function call (check_domain! after update_boundary_quantities!).

In addition, this PR contains:

• Three different mirroring methods (see images below).
• A fallback to zeroth-order interpolation when the determinant of the correction matrix drops below a critical threshold.

in the following we'll see some snapshots at t= 1.5.

On main

Tafuni inflow and outflow

scenario 0

bug: check_domain! after update_boundary_quantities!

# ...
density[particle] = f_d[1] + dot(pos_diff, df_d) # extrapolated
# ...
pressure[particle] = state_equation(density[particle])
# ...

This PR

Tafuni inflow and outflow

Zero-th order mirroring

First order mirroring

Simple mirroring

[codecov]

Codecov Report

❌ Patch coverage is 73.17073% with 55 lines in your changes missing coverage. Please review.
✅ Project coverage is 70.49%. Comparing base (9972df0) to head (fcfd477).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/schemes/boundary/open_boundary/mirroring.jl	73.79%	38 Missing ⚠️
...oundary/open_boundary/method_of_characteristics.jl	11.76%	15 Missing ⚠️
src/schemes/boundary/open_boundary/system.jl	0.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #855      +/-   ##
==========================================
+ Coverage   70.45%   70.49%   +0.03%     
==========================================
  Files         106      106              
  Lines        6895     7033     +138     
==========================================
+ Hits         4858     4958     +100     
- Misses       2037     2075      +38     

Flag	Coverage Δ
unit	70.49% <73.17%> (+0.03%)	⬆️

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.

[LasNikas]

Update

I added 3 different mirror methods and plotted the pressure along the center line of the entire domain

p(x) = x

p(x) = x^2

p(x) = cos(2pi * x)

Reproduce

using TrixiParticles
using Plots

function mirror(pressure_function, mirror_method;
                particle_spacing=0.05, domain_size=(2.0, 1.0))
    domain_fluid = RectangularShape(particle_spacing,
                                    round.(Int, domain_size ./ particle_spacing),
                                    (0.0, 0.0), density=1000.0,
                                    pressure=pressure_function)

    smoothing_length = 1.2 * particle_spacing
    smoothing_kernel = WendlandC2Kernel{2}()
    fluid_system = EntropicallyDampedSPHSystem(domain_fluid, smoothing_kernel,
                                               smoothing_length, 1.0)

    fluid_system.cache.density .= domain_fluid.density

    plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]])

    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, mirror_method,
                                       v_open_boundary, v_fluid,
                                       outflow.initial_condition.coordinates,
                                       domain_fluid.coordinates, semi, 0.0;
                                       prescribed_pressure=false)

    plane_in = ([0.0, 0.0], [0.0, domain_size[2]])

    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)

    TrixiParticles.set_zero!(open_boundary_in.pressure)

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

    return fluid_system, open_boundary_in, open_boundary_out, v_fluid
end

function interpolate_pressure(mirror_method, pressure_func; particle_spacing=0.05)
    fluid_system, open_boundary_in, open_boundary_out,
    v_fluid = mirror(pressure_func, mirror_method)

    p_fluid = [TrixiParticles.current_pressure(v_fluid, fluid_system, particle)
               for particle in TrixiParticles.active_particles(fluid_system)]

    fluid_system.initial_condition.pressure .= p_fluid
    open_boundary_in.initial_condition.pressure .= open_boundary_in.pressure
    open_boundary_out.initial_condition.pressure .= open_boundary_out.pressure

    entire_domain = union(fluid_system.initial_condition,
                          open_boundary_in.initial_condition,
                          open_boundary_out.initial_condition)

    smoothing_length = 1.2 * particle_spacing
    smoothing_kernel = WendlandC2Kernel{2}()

    # Use a fluid system to interpolate the pressure
    interpolation_system = WeaklyCompressibleSPHSystem(entire_domain,
                                                       ContinuityDensity(),
                                                       nothing, smoothing_kernel,
                                                       smoothing_length)
    interpolation_system.pressure .= entire_domain.pressure

    semi = Semidiscretization(interpolation_system)
    ode = semidiscretize(semi, (0, 0))
    v_ode, u_ode = ode.u0.x

    result = interpolate_line([-0.5, 0.5], [2.5, 0.5], 100, semi,
                              interpolation_system, v_ode, u_ode)

    return result.pressure
end

pressure_func_1(pos) = pos[1]
pressure_func_2(pos) = pos[1]^2
pressure_func_3(pos) = cos(2pi * pos[1])

pressures = interpolate_pressure.([
                                      SimpleMirroring(),
                                      FirstOrderMirroring(),
                                      ZerothOrderMirroring()
                                  ],
                                  pressure_func_3)

p = plot(Shape([0.0, 2.0, 2.0, 0.0],
               [
                   minimum(stack(pressures)),
                   minimum(stack(pressures)),
                   maximum(stack(pressures)),
                   maximum(stack(pressures))
               ]), color=nothing, label="fluid domain")

label_ = ["Simple Mirroring" "First Order Mirroring" "Zeroth Order Mirroring"]
linestyles = [:solid :dash :dot]

x = collect(range(-0.5, stop=2.5, length=length(first(pressures))))
plot!(p, x, pressures, label=label_, xlabel="x", ylabel="p", linewidth=2,
      linestyle=linestyles, legend=:outerright)

[LasNikas]

/run-gpu-tests

[efaulhaber]

I'll have to look at the tests again.

[LasNikas]

/run-gpu-tests

[LasNikas]

/run-gpu-tests

[LasNikas]

/run-gpu-tests