Open boundaries don't work if `UpdateCallback` is not the last callback in the list

[efaulhaber]

For some reason, we are getting 1e16 coordinates (which are used by inactive particles) into the neighborhood search update (which is only done for active particles):

┌ Info: 
│   point_coords =
│    2-element SVector{2, Float64} with indices SOneTo(2):
│     1.0e16
└     1.0e16
┌ Info: 
└   point in eachindex_y = true
ERROR: LoadError: particle coordinates are NaN or outside the domain bounds of the cell list
cell (3333333333333333505, 3333333333333333505) is out of bounds for cell grid of size [2:72, 2:37]

[efaulhaber]

This is happening when either a SolutionSavingCallback or a PostprocessCallback comes after the UpdateCallback. The error is thrown inside the get_du(integrator) that is called by the saving/pp callback.

[LasNikas]

I just tested. That works:

julia> trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
                     callbacks=CallbackSet(UpdateCallback(), InfoCallback(interval=100), SolutionSavingCallback(dt=0.02)))

[efaulhaber]

That doesn't:

using TrixiParticles
using OrdinaryDiffEqLowStorageRK

# ==========================================================================================
# ==== Resolution
n_particles_y = 9

# ==========================================================================================
# ==== Experiment Setup
tspan = (0.0, 2.0)

fin_length = 0.05
fin_thickness = 8e-3
poisson_ratio = 0.3
modulus = 1.5625e10

fiber_volume_fraction = 0.6
fiber_density = 1800.0
epoxy_density = 1250.0
density = fiber_volume_fraction * fiber_density +
          (1 - fiber_volume_fraction) * epoxy_density

clamp_radius = 0.05

tank_size = (0.2, 0.1)
center = (0.04, 0.05)
initial_fluid_size = tank_size
initial_velocity = (1.0, 0.0)

# The structure starts at the position of the first particle and ends
# at the position of the last particle.
particle_spacing = fin_thickness / (n_particles_y - 1)
fluid_particle_spacing = particle_spacing

smoothing_length_structure = sqrt(2) * particle_spacing
smoothing_length_fluid = 1.5 * fluid_particle_spacing
smoothing_kernel = WendlandC2Kernel{2}()

# Beam and clamped particles
length_clamp = round(Int, 0.005 / particle_spacing) * particle_spacing # m
n_particles_per_dimension = (round(Int, (length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x,
                             n_particles_y)
shape_sampled = RectangularShape(particle_spacing, n_particles_per_dimension,
                                 (-length_clamp, 0.0), density=density, place_on_shell=true)
length_clamp = 0.0
n_particles_per_dimension = (round(Int, (fin_length + length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x,
                             n_particles_y)

# Note that the `RectangularShape` puts the first particle half a particle spacing away
# from the boundary, which is correct for fluids, but not for structures.
# We therefore need to pass `place_on_shell=true`.
beam = RectangularShape(particle_spacing, n_particles_per_dimension,
                        (-length_clamp, 0.0), density=density, place_on_shell=true)

fixed_particles = setdiff(shape_sampled, beam)

# structure = union(beam, fixed_particles)

# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
boundary_layers = 3

# Make sure that the kernel support of fluid particles at an open boundary is always
# fully sampled.
# Note: Due to the dynamics at the inlets and outlets of open boundaries,
# it is recommended to use `open_boundary_layers > boundary_layers`
open_boundary_layers = 6

fluid_density = 1000.0
tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
                       n_layers=boundary_layers,
                       faces=(false, false, true, true), velocity=initial_velocity)
# fluid = setdiff(tank.fluid, structure)

open_boundary_size = (fluid_particle_spacing * open_boundary_layers, tank_size[2])

min_coords_inlet = (-open_boundary_layers * fluid_particle_spacing, 0.0)
inlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size,
                        fluid_density, n_layers=boundary_layers,
                        min_coordinates=min_coords_inlet,
                        faces=(false, false, true, true))

min_coords_outlet = (tank.fluid_size[1], 0.0)
outlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size,
                         fluid_density, n_layers=boundary_layers,
                         min_coordinates=min_coords_outlet,
                         faces=(false, false, true, true))


NDIMS = ndims(tank.fluid)
n_buffer_particles = 10 * tank.n_particles_per_dimension[2]^(NDIMS - 1)


structure = union(beam, fixed_particles)
# Move the fin to the center of the tank
structure.coordinates .+= center .+ (fluid_particle_spacing / 2, fluid_particle_spacing / 2)

fluid = setdiff(tank.fluid, structure)

n_clamped_particles = nparticles(structure) - nparticles(beam)

sound_speed = 40.0
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
                                   exponent=1, background_pressure=0.0)

# ==========================================================================================
# ==== Structure
boundary_density_calculator = AdamiPressureExtrapolation()
viscosity_fluid = ViscosityAdami(nu=1e-4)
viscosity_fin = ViscosityAdami(nu=1e-4)

# For the FSI we need the hydrodynamic masses and densities in the structure boundary model
hydrodynamic_densites = fluid_density * ones(size(structure.density))
hydrodynamic_masses = hydrodynamic_densites * particle_spacing^2

boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densites,
                                                   hydrodynamic_masses,
                                                   state_equation=state_equation,
                                                   boundary_density_calculator,
                                                   smoothing_kernel, smoothing_length_fluid,
                                                   viscosity=viscosity_fin)

structure_system = WallBoundarySystem(structure, boundary_model_structure)

# ==========================================================================================
# ==== Fluid
fluid_density_calculator = ContinuityDensity()
density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
# density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1)

fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
                                           state_equation, smoothing_kernel,
                                           smoothing_length_fluid, viscosity=viscosity_fluid,
                                           density_diffusion=density_diffusion,
                                           shifting_technique=ParticleShiftingTechnique(sound_speed_factor=0.2, v_max_factor=0.0),
                                        #    pressure_acceleration=tensile_instability_control,
                                           buffer_size=n_buffer_particles)
# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
#                                            sound_speed, viscosity=ViscosityAdami(; nu),
#                                            transport_velocity=TransportVelocityAdami(10 * sound_speed^2 * fluid_density))

# ==========================================================================================
# ==== Open Boundaries
periodic = false
if periodic
    min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2
    max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2
    min_corner = convert.(typeof(fluid_particle_spacing), min_corner)
    max_corner = convert.(typeof(fluid_particle_spacing), max_corner)
    periodic_box = PeriodicBox(; min_corner, max_corner)
    open_boundary_system = nothing
    wall = tank.boundary
else
    periodic_box = nothing

    open_boundary_model = BoundaryModelDynamicalPressureZhang()
    # open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring())
    reference_velocity_in = [1.0, 0.0]
    reference_pressure_in = nothing
    reference_density_in = nothing
    boundary_type_in = InFlow()
    face_in = ([0.0, 0.0], [0.0, tank_size[2]])
    flow_direction = [1.0, 0.0]
    inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction,
                        open_boundary_layers, density=fluid_density, particle_spacing,
                        reference_density=reference_density_in,
                        reference_pressure=reference_pressure_in,
                        reference_velocity=reference_velocity_in,
                        initial_condition=inlet.fluid, boundary_type=boundary_type_in)

    reference_velocity_out = nothing
    reference_pressure_out = 0.0
    reference_density_out = nothing
    boundary_type_out = OutFlow()
    face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], tank_size[2]])
    outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction),
                        open_boundary_layers, density=fluid_density, particle_spacing,
                        reference_density=reference_density_out,
                        reference_pressure=reference_pressure_out,
                        reference_velocity=reference_velocity_out,
                        initial_condition=outlet.fluid, boundary_type=boundary_type_out)

    open_boundary_system = OpenBoundarySystem(inflow, outflow; fluid_system,
                                    boundary_model=open_boundary_model,
                                    buffer_size=n_buffer_particles,
                                    shifting_technique=nothing)

    wall = union(tank.boundary, inlet.boundary, outlet.boundary)
end

# ==========================================================================================
# ==== Boundary
boundary_density_calculator = AdamiPressureExtrapolation()
boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass,
                                             state_equation=state_equation,
                                             boundary_density_calculator,
                                             smoothing_kernel, smoothing_length_fluid)

boundary_system = WallBoundarySystem(wall, boundary_model)

# ==========================================================================================
# ==== Simulation
min_corner = minimum(wall.coordinates, dims=2) .- fluid_particle_spacing / 2
max_corner = maximum(wall.coordinates, dims=2) .+ fluid_particle_spacing / 2
cell_list = FullGridCellList(; min_corner, max_corner)
neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list,
                                                update_strategy=ParallelUpdate())

semi = Semidiscretization(fluid_system, boundary_system, open_boundary_system, structure_system; neighborhood_search,
                          parallelization_backend=PolyesterBackend())
ode = semidiscretize(semi, tspan)

info_callback = InfoCallback(interval=100)
saving_callback = SolutionSavingCallback(dt=0.01, prefix="")

function total_volume(system::WeaklyCompressibleSPHSystem, data, t)
    return sum(data.mass ./ data.density)
end
function total_volume(system, data, t)
    return nothing
end
pp_cb = PostprocessCallback(; total_volume, interval=100,
                            filename="total_volume", write_file_interval=50)

callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), pp_cb)

dt_fluid = 2e-5
sol = solve(ode,
            CarpenterKennedy2N54(williamson_condition=false),
            dt=dt_fluid, # This is overwritten by the stepsize callback
            save_everystep=false, callback=callbacks, maxiters=10^8);