diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 3e2870628..9912b9172 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -423,12 +423,6 @@ function compute_pressure!(boundary_model, boundary_pressure_extrapolation!(Val(parallelize), boundary_model, system, neighbor_system, system_coords, neighbor_coords, v, v_neighbor_system, semi) - - @threaded semi for particle in eachparticle(system) - # Limit pressure to be non-negative to avoid attractive forces between fluid and - # boundary particles at free surfaces (sticking artifacts). - pressure[particle] = max(pressure[particle], 0) - end end @trixi_timeit timer() "inverse state equation" @threaded semi for particle in @@ -447,13 +441,17 @@ function compute_adami_density!(boundary_model, system, system_coords, particle) # The summation is only over fluid particles, thus the volume stays zero when a boundary # particle isn't surrounded by fluid particles. # Check the volume to avoid NaNs in pressure and velocity. - if volume[particle] > eps() - pressure[particle] /= volume[particle] + if @inbounds volume[particle] > eps() + @inbounds pressure[particle] /= volume[particle] # To impose no-slip condition compute_wall_velocity!(viscosity, system, system_coords, particle) end + # Limit pressure to be non-negative to avoid attractive forces between fluid and + # boundary particles at free surfaces (sticking artifacts). + @inbounds pressure[particle] = max(pressure[particle], 0) + # Apply inverse state equation to compute density (not used with EDAC) inverse_state_equation!(density, state_equation, pressure, particle) end @@ -519,26 +517,30 @@ end v_neighbor_system, particle, neighbor, pos_diff, distance, viscosity, cache, pressure, pressure_offset) - density_neighbor = current_density(v_neighbor_system, neighbor_system, neighbor) + density_neighbor = @inbounds current_density(v_neighbor_system, neighbor_system, + neighbor) + # Fluid pressure term + fluid_pressure = @inbounds current_pressure(v_neighbor_system, neighbor_system, + neighbor) + + # Hydrostatic pressure term from fluid and boundary acceleration resulting_acceleration = neighbor_system.acceleration - - current_acceleration(system, particle) + @inbounds current_acceleration(system, particle) + hydrostatic_pressure = dot(resulting_acceleration, density_neighbor * pos_diff) - kernel_weight = smoothing_kernel(boundary_model, distance, particle) + # Additional dynamic pressure term (only with `BernoulliPressureExtrapolation`) + dynamic_pressure_ = dynamic_pressure(boundary_density_calculator, density_neighbor, + v, v_neighbor_system, pos_diff, distance, + particle, neighbor, system, neighbor_system) - pressure[particle] += (pressure_offset - + - current_pressure(v_neighbor_system, neighbor_system, - neighbor) - + - dynamic_pressure(boundary_density_calculator, density_neighbor, - v, v_neighbor_system, pos_diff, distance, - particle, neighbor, system, neighbor_system) - + - dot(resulting_acceleration, density_neighbor * pos_diff)) * - kernel_weight + sum_pressures = pressure_offset + fluid_pressure + dynamic_pressure_ + + hydrostatic_pressure + + kernel_weight = smoothing_kernel(boundary_model, distance, particle) - cache.volume[particle] += kernel_weight + @inbounds pressure[particle] += sum_pressures * kernel_weight + @inbounds cache.volume[particle] += kernel_weight compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, kernel_weight, particle, neighbor) @@ -586,8 +588,8 @@ function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, neighbor_s v_neighbor_system, kernel_weight, particle, neighbor) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - for dim in 1:ndims(neighbor_system) - cache.wall_velocity[dim, particle] += kernel_weight * v_b[dim] + for dim in eachindex(v_b) + @inbounds cache.wall_velocity[dim, particle] += kernel_weight * v_b[dim] end return cache @@ -606,17 +608,17 @@ end # This velocity is zero when not using moving boundaries. v_boundary = current_velocity(system_coords, system, particle) - for dim in 1:ndims(system) - # The second term is the precalculated smoothed velocity field of the fluid. - wall_velocity[dim, - particle] = 2 * v_boundary[dim] - - wall_velocity[dim, particle] / volume[particle] + for dim in eachindex(v_boundary) + # The second term is the precalculated smoothed velocity field of the fluid + new_velocity = @inbounds 2 * v_boundary[dim] - + wall_velocity[dim, particle] / volume[particle] + @inbounds wall_velocity[dim, particle] = new_velocity end return viscosity end @inline function inverse_state_equation!(density, state_equation, pressure, particle) - density[particle] = inverse_state_equation(state_equation, pressure[particle]) + @inbounds density[particle] = inverse_state_equation(state_equation, pressure[particle]) return density end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 4027cf756..acf269d7b 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -245,6 +245,8 @@ end function update_positions!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t) (; current_coordinates) = system + # `current_coordinates` stores the coordinates of both integrated and fixed particles. + # Copy the coordinates of the integrated particles from `u`. @threaded semi for particle in each_moving_particle(system) for i in 1:ndims(system) current_coordinates[i, particle] = u[i, particle]