diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c316cd801..8193ec9cf 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -45,6 +45,7 @@ include("callbacks/callbacks.jl") include("general/general.jl") include("setups/setups.jl") include("schemes/schemes.jl") +include("refinement/refinement.jl") # Note that `semidiscretization.jl` depends on the system types and has to be # included separately. `gpu.jl` in turn depends on the semidiscretization type. @@ -72,7 +73,7 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BernoulliPressureExtrapolation - +export ParticleRefinement, SpatialRefinementCriterion export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 1aefdff72..3d8f5701a 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -31,3 +31,4 @@ include("density_reinit.jl") include("post_process.jl") include("stepsize.jl") include("update.jl") +include("refinement.jl") diff --git a/src/callbacks/refinement.jl b/src/callbacks/refinement.jl new file mode 100644 index 000000000..130de8756 --- /dev/null +++ b/src/callbacks/refinement.jl @@ -0,0 +1,129 @@ +struct ParticleRefinementCallback{I} + interval::I +end + +function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) + if dt > 0 && interval !== -1 + throw(ArgumentError("setting both `interval` and `dt` is not supported")) + end + + # Update in intervals in terms of simulation time + if dt > 0 + interval = Float64(dt) + + # Update every time step (default) + elseif interval == -1 + interval = 1 + end + + refinement_callback = ParticleRefinementCallback(interval) + + if dt > 0 + # Add a `tstop` every `dt`, and save the final solution. + return PeriodicCallback(refinement_callback, dt, + initialize=initial_update!, + save_positions=(false, false)) + else + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(refinement_callback, refinement_callback, + initialize=initial_update!, + save_positions=(false, false)) + end +end + +# initialize +function initial_refinement!(cb, u, t, integrator) + # The `ParticleRefinementCallback` is either `cb.affect!` (with `DiscreteCallback`) + # or `cb.affect!.affect!` (with `PeriodicCallback`). + # Let recursive dispatch handle this. + + initial_refinement!(cb.affect!, u, t, integrator) +end + +function initial_refinement!(cb::ParticleRefinementCallback, u, t, integrator) + semi = integrator.p + + foreach_system(semi) do system + resize_refinement!(system) + end + + cb(integrator) +end + +# condition +function (refinement_callback::ParticleRefinementCallback)(u, t, integrator) + (; interval) = refinement_callback + + return condition_integrator_interval(integrator, interval) +end + +# affect +function (refinement_callback::ParticleRefinementCallback)(integrator) + t = integrator.t + semi = integrator.p + v_ode, u_ode = integrator.u.x + + # Update NHS + @trixi_timeit timer() "update nhs" update_nhs(u_ode, semi) + + # Basically `get_tmp_cache(integrator)` to write into in order to be non-allocating + # https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Caches + v_tmp, u_tmp = integrator.cache.tmp.x + + v_tmp .= v_ode + u_tmp .= u_ode + + refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + + resize!(integrator, (length(v_ode), length(u_ode))) + + # Tell OrdinaryDiffEq that u has been modified + u_modified!(integrator, true) + + return integrator +end + +Base.resize!(a::RecursiveArrayTools.ArrayPartition, sizes::Tuple) = resize!.(a.x, sizes) + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(interval=", (cb.affect!.interval), ")") +end + +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(dt=", cb.affect!.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + refinement_cb = cb.affect! + setup = [ + "interval" => refinement_cb.interval + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + refinement_cb = cb.affect!.affect! + setup = [ + "dt" => refinement_cb.interval + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end diff --git a/src/general/corrections.jl b/src/general/corrections.jl index b1c66a06a..c63c25a27 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -220,7 +220,7 @@ function compute_correction_values!(system, kernel_correction_coefficient[particle] += volume * smoothing_kernel(system, distance) if distance > sqrt(eps()) - tmp = volume * smoothing_kernel_grad(system, pos_diff, distance) + tmp = volume * smoothing_kernel_grad(system, pos_diff, distance, particle) for i in axes(dw_gamma, 1) dw_gamma[i, particle] += tmp[i] end @@ -311,7 +311,7 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, pos_diff, distance volume = mass[neighbor] / density_fun(neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) iszero(grad_kernel) && return @@ -349,7 +349,7 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance, smoothing_length, system, particle) - return smoothing_kernel_grad(system, pos_diff, distance) + return smoothing_kernel_grad(system, pos_diff, distance, particle) end # Compute gradient of corrected kernel diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 2a3d45ea0..e15eb33e0 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -66,7 +66,7 @@ function summation_density!(system, semi, u, u_ode, density; points=particles) do particle, neighbor, pos_diff, distance mass = hydrodynamic_mass(neighbor_system, neighbor) - density[particle] += mass * smoothing_kernel(system, distance) + density[particle] += mass * smoothing_kernel(system, distance, particle) end end end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..6c0c7e37b 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -134,8 +134,8 @@ function create_neighborhood_search(neighborhood_search, system, neighbor) end @inline function compact_support(system, neighbor) - (; smoothing_kernel, smoothing_length) = system - return compact_support(smoothing_kernel, smoothing_length) + (; smoothing_kernel) = system + return compact_support(smoothing_kernel, smoothing_length(system, 1)) # TODO end @inline function compact_support(system::OpenBoundarySPHSystem, neighbor) @@ -407,7 +407,9 @@ end function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) (; systems) = semi - return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) + return 1.0 + # TODO + # return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) end function drift!(du_ode, v_ode, u_ode, semi, t) diff --git a/src/general/system.jl b/src/general/system.jl index 48d556eba..88b9e7bf2 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -106,30 +106,23 @@ end @inline set_particle_pressure!(v, system, particle, pressure) = v -@inline function smoothing_kernel(system, distance) - (; smoothing_kernel, smoothing_length) = system - return kernel(smoothing_kernel, distance, smoothing_length) +@inline function smoothing_kernel(system, distance, particle) + (; smoothing_kernel) = system + return kernel(smoothing_kernel, distance, smoothing_length(system, particle)) end -@inline function smoothing_kernel_deriv(system, distance) - (; smoothing_kernel, smoothing_length) = system - return kernel_deriv(smoothing_kernel, distance, smoothing_length) -end - -@inline function smoothing_kernel_grad(system, pos_diff, distance) - return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) -end - -@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance) +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) (; smoothing_kernel, smoothing_length) = system.boundary_model return kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) end @inline function smoothing_kernel_grad(system, pos_diff, distance, particle) - return corrected_kernel_grad(system.smoothing_kernel, pos_diff, distance, - system.smoothing_length, system.correction, system, - particle) + (; smoothing_kernel) = system + + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length(system, particle), + system.correction, system, particle) end # System update orders. This can be dispatched if needed. diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl new file mode 100644 index 000000000..fe1f90b10 --- /dev/null +++ b/src/refinement/refinement.jl @@ -0,0 +1,138 @@ +include("refinement_criteria.jl") + +struct ParticleRefinement{RC, ELTYPE} + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles +end + +function ParticleRefinement(; max_spacing_ratio, + refinement_criteria=SpatialRefinementCriterion()) + mass_ref = Vector{eltype(max_spacing_ratio)}() + + if !(refinement_criteria isa Tuple) + refinement_criteria = (refinement_criteria,) + end + return ParticleRefinement(refinement_criteria, max_spacing_ratio, mass_ref) +end + +resize_refinement!(system) = system + +function resize_refinement!(system::FluidSystem) + resize_refinement!(system, system.particle_refinement) +end + +resize_refinement!(system, ::Nothing) = system + +function resize_refinement!(system, particle_refinement) + resize!(particle_refinement.mass_ref, nparticles(system)) + + return system +end + +function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + check_refinement_criteria!(semi, v_ode, u_ode) + + # Update the spacing of particles (Algorthm 1) + update_particle_spacing(semi, v_ode, u_ode) + + # Split the particles (Algorithm 2) + + # Merge the particles (Algorithm 3) + + # Shift the particles + + # Correct the particles + + # Update smoothing lengths + + # Resize neighborhood search + + return semi +end + +function update_particle_spacing(semi, v_ode, u_ode) + foreach_system(semi) do system + update_particle_spacing(system, v_ode, u_ode, semi) + end +end + +# The methods for the `FluidSystem` are defined in `src/schemes/fluid/fluid.jl` +@inline update_particle_spacing(system, v_ode, u_ode, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, v_ode, u_ode, semi) + update_particle_spacing(system, system.particle_refinement, v_ode, u_ode, semi) +end + +@inline update_particle_spacing(system::FluidSystem, ::Nothing, v_ode, u_ode, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, particle_refinement, + v_ode, u_ode, semi) + (; smoothing_length, smoothing_length_factor) = system.cache + (; mass_ref, max_spacing_ratio) = particle_refinement + + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, semi) + + system_coords = current_coordinates(u, system) + + for particle in eachparticle(system) + dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, u_ode, system_coords, + particle) + + if dp_max / dp_min < max_spacing_ratio^3 + new_spacing = min(dp_max, max_spacing_ratio * dp_min) + else + new_spacing = dp_avg + end + + smoothing_length[particle] = smoothing_length_factor * new_spacing + mass_ref[particle] = particle_density(v, system, particle) * + new_spacing^(ndims(system)) + end + + return system +end + + +@inline function min_max_avg_spacing(system, semi, u_ode, system_coords, particle) + dp_min = Inf + dp_max = zero(eltype(system)) + dp_avg = zero(eltype(system)) + counter_neighbors = 0 + + foreach_system(semi) do neighbor_system + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) + + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + PointNeighbors.foreach_neighbor(system_coords, neighbor_coords, neighborhood_search, + particle) do particle, neighbor, pos_diff, distance + dp_neighbor = particle_spacing(neighbor_system, neighbor) + + dp_min = min(dp_min, dp_neighbor) + dp_max = max(dp_max, dp_neighbor) + dp_avg += dp_neighbor + + counter_neighbors += 1 + end + end + + dp_avg / counter_neighbors + + return dp_min, dp_max, dp_avg +end + +@inline particle_spacing(system, particle) = system.initial_condition.particle_spacing + +@inline function particle_spacing(system::FluidSystem, particle) + return particle_spacing(system, system.particle_refinement, particle) +end + +@inline particle_spacing(system, ::Nothing, _) = system.initial_condition.particle_spacing + +@inline function particle_spacing(system, refinement, particle) + (; smoothing_length_factor) = system.cache + return smoothing_length(system, particle) / smoothing_length_factor +end diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl new file mode 100644 index 000000000..627ffb619 --- /dev/null +++ b/src/refinement/refinement_criteria.jl @@ -0,0 +1,69 @@ +abstract type RefinementCriteria end + +struct SpatialRefinementCriterion <: RefinementCriteria end + +struct SolutionRefinementCriterion <: RefinementCriteria end + +function check_refinement_criteria!(semi, v_ode, u_ode) + foreach_system(semi) do system + check_refinement_criteria!(system, v_ode, u_ode, semi) + end +end + +@inline check_refinement_criteria!(system, v_ode, u_ode, semi) = system + +@inline function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + check_refinement_criteria!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) +end + +@inline function check_refinement_criteria!(system::FluidSystem, ::Nothing, + v, u, v_ode, u_ode, semi) + system +end + +@inline function check_refinement_criteria!(system::FluidSystem, refinement, + v, u, v_ode, u_ode, semi) + (; refinement_criteria) = system.particle_refinement + for criterion in refinement_criteria + criterion(system, v, u, v_ode, u_ode, semi) + end +end + +@inline function (criterion::SpatialRefinementCriterion)(system, v, u, v_ode, u_ode, semi) + system_coords = current_coordinates(u, system) + + foreach_system(semi) do neighbor_system + set_particle_spacing!(system, neighbor_system, system_coords, v_ode, u_ode, semi) + end + return system +end + +@inline set_particle_spacing!(system, _, _, _, _, _) = system + +@inline function set_particle_spacing!(particle_system, + neighbor_system::Union{BoundarySystem, SolidSystem}, + system_coords, v_ode, u_ode, semi) + (; smoothing_length, smoothing_length_factor) = particle_system.cache + + neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, semi) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + foreach_point_neighbor(particle_system, neighbor_system, + system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + + dp_neighbor = particle_spacing(neighbor_system, neighbor) + dp_particle = smoothing_length[particle] / smoothing_length_factor + + smoothing_length[particle] = smoothing_length_factor * min(dp_neighbor, dp_particle) + end + + return particle_system +end diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 8b7e507ed..ece912eeb 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -69,6 +69,10 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, smoothing_length, viscosity, correction, cache) end +function smoothing_length(boundary_model::BoundaryModelDummyParticles, particle) + return boundary_model.smoothing_length +end + @doc raw""" AdamiPressureExtrapolation(; pressure_offset=0.0) @@ -480,7 +484,7 @@ end resulting_acceleration = neighbor_system.acceleration - current_acceleration(system, particle) - kernel_weight = smoothing_kernel(boundary_model, distance) + kernel_weight = smoothing_kernel(boundary_model, distance, particle) pressure[particle] += (pressure_offset + @@ -581,12 +585,13 @@ end return density end -@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) - (; smoothing_kernel, smoothing_length, correction) = system.boundary_model +# TODO +# @inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) +# (; smoothing_kernel, smoothing_length, correction) = system.boundary_model - return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, - smoothing_length, correction, system, particle) -end +# return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, +# smoothing_length, correction, system, particle) +# end @inline function correction_matrix(system::BoundarySystem, particle) extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle) diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 59f634092..d90d2a291 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -29,7 +29,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_diff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance, particle) continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff, grad_kernel, particle) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 3644c3630..b075813f2 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -30,9 +30,10 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - dv_pressure = pressure_acceleration(particle_system, neighbor_system, neighbor, + dv_pressure = pressure_acceleration(particle_system, neighbor_system, + particle, neighbor, m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) @@ -44,20 +45,29 @@ function interact!(dv, v_particle_system, u_particle_system, # Add convection term when using `TransportVelocityAdami` dv_convection = momentum_convection(particle_system, neighbor_system, + pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) + dv_velocity_correction = velocity_correction(particle_system, neighbor_system, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + + dv_velocity_correction[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) - pressure_evolution!(dv, particle_system, v_diff, grad_kernel, - particle, pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) + pressure_evolution!(dv, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, v_diff, grad_kernel, + particle, neighbor, pos_diff, distance, sound_speed, + m_a, m_b, p_a, p_b, rho_a, rho_b) transport_velocity!(dv, particle_system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) @@ -69,11 +79,47 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end -@inline function pressure_evolution!(dv, particle_system, v_diff, grad_kernel, particle, - pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) - (; smoothing_length) = particle_system +@inline function pressure_evolution!(dv, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + v_diff, grad_kernel, particle, neighbor, + pos_diff, distance, sound_speed, + m_a, m_b, p_a, p_b, rho_a, rho_b) + (; particle_refinement) = particle_system + + # This is basically the continuity equation times `sound_speed^2` + artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) + + beta_inv_a = beta_correction(particle_system, particle) + beta_inv_b = beta_correction(neighbor_system, neighbor) + + dv[end, particle] += beta_inv_a * artificial_eos + + pressure_damping_term(particle_system, neighbor_system, + particle_refinement, particle, neighbor, + pos_diff, distance, beta_inv_a, m_a, m_b, + p_a, p_b, rho_b, rho_b, grad_kernel) + + pressure_reduction(particle_system, neighbor_system, + particle_refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_inv_a, beta_inv_b) + + return dv +end +@inline beta_correction(system, particle) = one(eltype(system)) + +@inline function beta_correction(system::FluidSystem, particle) + beta_correction(system, system.particle_refinement, particle) +end + +@inline beta_correction(particle_system, ::Nothing, particle) = one(eltype(particle_system)) + +@inline function beta_correction(particle_system, refinement, particle) + return inv(particle_system.cache.beta[particle]) +end +function pressure_damping_term(particle_system, neighbor_system, ::Nothing, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel) volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -81,15 +127,13 @@ end # EDAC pressure evolution pressure_diff = p_a - p_b - # This is basically the continuity equation times `sound_speed^2` - artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) - eta_a = rho_a * particle_system.nu_edac eta_b = rho_b * particle_system.nu_edac eta_tilde = 2 * eta_a * eta_b / (eta_a + eta_b) - # TODO For variable smoothing length use average smoothing length - tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length^2) + smoothing_length_average = 0.5 * (smoothing_length(particle_system, particle) + + smoothing_length(particle_system, particle)) + tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length_average^2) # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 # They argued that the formulation is more flexible because of the possibility to formulate @@ -102,11 +146,99 @@ end # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 # # This is similar to density diffusion in WCSPH - damping_term = volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) + return volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) +end - dv[end, particle] += artificial_eos + damping_term +function pressure_damping_term(particle_system, neighbor_system, refinement, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel_a) + # EDAC pressure evolution + pressure_diff = p_a - p_b - return dv + # Haftu et al. (2022) uses `alpha_edac = 1.5` in all their simulations + alpha_edac = 1.5 + + # TODO: Haftu et al. (2022) use `8` but I think it depeneds on the dimension (see Monaghan, 2005) + tmp = 2 * ndims(particle_system) + 4 + + nu_edac_a = alpha_edac * sound_speed * smoothing_length(particle_system, particle) / tmp + nu_edac_a = alpha_edac * sound_speed * smoothing_length(neighbor_system, particle) / tmp + + nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) + + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + grad_W_avg = 0.5 * (grad_kernel_a + grad_kernel_b) + + return beta_inv_a * nu_edac_ab * pressure_diff * dot(pos_diff, grad_W_avg) * m_b / rho_b +end + +function pressure_reduction(particle_system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_a, beta_b) + return zero(eltype(particle_system)) +end + +function pressure_reduction(particle_system, neighbor_system, refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_a, beta_b) + beta_inv_a = beta_correction(particle_system, particle) + beta_inv_b = beta_correction(neighbor_system, neighbor) + + p_a_avg = average_pressure(particle_system, particle) + p_b_avg = average_pressure(neighbor_system, neighbor) + + P_a = (p_a - p_a_avg) / (rho_a^2 * beta_inv_a) + P_b = (p_b - p_b_avg) / (rho_b^2 * beta_inv_b) + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + v_diff = advection_velocity(v_particle_system, particle_system, particle) - + current_velocity(v_particle_system, particle_system, particle) + + return m_b * (dot(v_diff, P_a * grad_kernel_a + P_b * grad_kernel_b)) +end + +@inline function velocity_correction(particle_system, neighbor_system, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + velocity_correction(particle_system, neighbor_system, + particle_system.particle_refinement, + pos_diff, distance, v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) +end + +@inline function velocity_correction(particle_system, neighbor_system, ::Nothing, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + return zero(pos_diff) +end + +@inline function velocity_correction(particle_system, neighbor_system, + particle_refinement, pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, + grad_kernel) + momentum_velocity_a = current_velocity(v_particle_system, system, particle) + advection_velocity_a = advection_velocity(v_particle_system, system, particle) + + momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + + v_diff = momentum_velocity_a - momentum_velocity_b + v_diff_tilde = advection_velocity_a - advection_velocity_b + + beta_inv_a = beta_correction(particle_system, particle) + + return -m_b * beta_inv_a * + (dot(v_diff_tilde - v_diff, grad_kernel) * momentum_velocity_a) / rho_b end # We need a separate method for EDAC since the density is stored in `v[end-1,:]`. @@ -125,3 +257,39 @@ end grad_kernel) return dv end + +# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + neighbor_system, particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, correction) + return pressure_acceleration(particle_system, particle_system.particle_refinement, + neighbor_system, particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, + correction) +end + +# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + ::Nothing, particle, neighbor_system, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, correction) + (; pressure_acceleration_formulation) = particle_system + return pressure_acceleration_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) +end + +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + particle_refinement, neighbor_system, particle, + neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, + pos_diff, distance, W_a, correction) + p_a_avg = average_pressure(particle_system, particle) + p_b_avg = average_pressure(neighbor_system, neighbor) + + P_a = beta_correction(particle_system, particle) * (p_a - p_a_avg) / rho_a^2 + P_b = beta_correction(neighbor_system, neighbor) * (p_b - p_b_avg) / rho_b^2 + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + return -m_b * (P_a * grad_kernel_a + P_b * grad_kernel_b) +end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -44,12 +44,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K - smoothing_length :: ELTYPE sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE @@ -59,55 +58,57 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, transport_velocity :: TV source_terms :: ST buffer :: B + particle_refinement :: PR cache :: C +end - function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, - smoothing_length, sound_speed; - pressure_acceleration=inter_particle_averaged_pressure, - density_calculator=SummationDensity(), - transport_velocity=nothing, - alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, - ndims(smoothing_kernel)), - source_terms=nothing, buffer_size=nothing) - buffer = isnothing(buffer_size) ? nothing : - SystemBuffer(nparticles(initial_condition), buffer_size) +function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, sound_speed; + pressure_acceleration=inter_particle_averaged_pressure, + density_calculator=SummationDensity(), + transport_velocity=nothing, + alpha=0.5, viscosity=nothing, + acceleration=ntuple(_ -> 0.0, + ndims(smoothing_kernel)), + particle_refinement=nothing, + source_terms=nothing, buffer_size=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) - initial_condition = allocate_buffer(initial_condition, buffer) + initial_condition = allocate_buffer(initial_condition, buffer) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) - mass = copy(initial_condition.mass) + mass = copy(initial_condition.mass) - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end + if ndims(smoothing_kernel) != NDIMS + throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) + end - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) + end - pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, - density_calculator, - NDIMS, ELTYPE, - nothing) + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, + density_calculator, + NDIMS, ELTYPE, + nothing) - nu_edac = (alpha * smoothing_length * sound_speed) / 8 + nu_edac = (alpha * smoothing_length * sound_speed) / 8 - cache = create_cache_density(initial_condition, density_calculator) - cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = create_cache_density(initial_condition, density_calculator) + cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) - new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), - typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), - typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms), - typeof(buffer), - typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, transport_velocity, source_terms, - buffer, cache) - end + return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, + smoothing_kernel, sound_speed, viscosity, nu_edac, + acceleration_, nothing, pressure_acceleration, + transport_velocity, source_terms, buffer, + particle_refinement, cache) end function Base.show(io::IO, system::EntropicallyDampedSPHSystem) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fd0a9c3a5..f287cabba 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -13,8 +13,30 @@ function create_cache_density(ic, ::ContinuityDensity) return (;) end +function create_cache_refinement(initial_condition, ::Nothing, smoothing_length) + return (; smoothing_length) +end + +function create_cache_refinement(initial_condition, refinement, smoothing_length) + smoothng_length_factor = smoothing_length / initial_condition.particle_spacing + return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)), + smoothing_length_factor=smoothng_length_factor) +end + @propagate_inbounds hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] +function smoothing_length(system::FluidSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::FluidSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::FluidSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + function write_u0!(u0, system::FluidSystem) (; initial_condition) = system @@ -89,7 +111,7 @@ include("entropically_damped_sph/entropically_damped_sph.jl") add_velocity!(du, v, particle, system, system.transport_velocity) end -@inline function momentum_convection(system, neighbor_system, +@inline function momentum_convection(system, neighbor_system, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) return zero(grad_kernel) @@ -98,9 +120,11 @@ end @inline function momentum_convection(system, neighbor_system::Union{EntropicallyDampedSPHSystem, WeaklyCompressibleSPHSystem}, + pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) momentum_convection(system, neighbor_system, system.transport_velocity, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) + system.particle_refinement, pos_diff, distance, v_particle_system, + v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, + grad_kernel) end diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl index 5b0b82e11..62105cf6c 100644 --- a/src/schemes/fluid/pressure_acceleration.jl +++ b/src/schemes/fluid/pressure_acceleration.jl @@ -107,7 +107,7 @@ function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing end # Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. -@inline function pressure_acceleration(particle_system, neighbor_system, neighbor, +@inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, correction) (; pressure_acceleration_formulation) = particle_system @@ -118,7 +118,7 @@ end end # Formulation using asymmetric gradient formulation for corrections depending on local neighborhood. -@inline function pressure_acceleration(particle_system, neighbor_system, neighbor, +@inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, correction::Union{KernelCorrection, diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 66140775e..9b736795e 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -105,10 +105,8 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance m_b = hydrodynamic_mass(neighbor_system, neighbor) - density_neighbor = particle_density(v_neighbor_system, - neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, - particle) + density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] * smoothing_length diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bc8fdf005..5270d1210 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -76,13 +76,14 @@ end return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) end -@inline function momentum_convection(system, neighbor_system, ::Nothing, +@inline function momentum_convection(system, neighbor_system, ::Nothing, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) end @inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + ::Nothing, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) volume_a = m_a / rho_a @@ -101,6 +102,28 @@ end return volume_term * (0.5 * (A_a + A_b)) * grad_kernel end +@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + refinement, pos_diff, distance, + v_particle_system, v_neighbor_system, rho_a, rho_b, + m_a, m_b, particle, neighbor, grad_kernel) + momentum_velocity_a = current_velocity(v_particle_system, system, particle) + advection_velocity_a = advection_velocity(v_particle_system, system, particle) + + momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + + factor_a = beta_correction(particle_system, particle) / rho_a + factor_b = beta_correction(neighbor_system, neighbor) / rho_b + + A_a = factor_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' + A_b = factor_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' + + grad_kernel_a = smoothing_kernel_grad(system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + return -m_b * (A_a * grad_kernel_a + A_b * grad_kernel_b) +end + @inline transport_velocity!(dv, system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) = dv @inline function transport_velocity!(dv, system::FluidSystem, diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1a0665a71..222fa6fd8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -112,24 +112,28 @@ end v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, - m_a, m_b, rho_a, rho_b, - grad_kernel) - (; smoothing_length) = particle_system - - rho_mean = (rho_a + rho_b) / 2 + sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + rho_mean = 0.5 * (rho_a + rho_b) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + + # TODO is that correct with the smoothing lengths? nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system)) + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system)) + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor) + # TODO do we need the average smoothing length here? pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, - smoothing_length, grad_kernel, nu_a, nu_b) + smoothing_length_particle, grad_kernel, nu_a, nu_b) return m_b * pi_ab end @@ -170,12 +174,12 @@ end # Joseph J. Monaghan. "Smoothed Particle Hydrodynamics". # In: Reports on Progress in Physics (2005), pages 1703-1759. # [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01) -function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan) - (; smoothing_length) = system +function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan, + smoothing_length_particle) (; alpha) = viscosity sound_speed = system_sound_speed(system) - return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4) + return alpha * smoothing_length_particle * sound_speed / (2 * ndims(system) + 4) end @doc raw""" @@ -213,13 +217,50 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - (; smoothing_length) = particle_system + viscosity(particle_system, neighbor_system, particle_system.particle_refinement, + v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +end + +@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, + particle_refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + beta_inv_a = beta_correction(particle_system, particle_refinement, particle) - epsilon = viscosity.epsilon nu_a = kinematic_viscosity(particle_system, viscosity_model(neighbor_system, particle_system)) + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + grad_W_avg = 0.5 * (grad_kernel_a + grad_kernel_b) + + tmp = (distance^2 + 0.001 * smoothing_length(particle_system, particle)^2) + + return m_b * beta_inv_a * 4 * nu_a * dot(pos_diff, grad_W_avg) * v_diff / + (tmp * (rho_a + rho_b)) +end + +@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + epsilon = viscosity.epsilon + + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + + # TODO is that correct with the smoothing lengths? + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system)) + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) @@ -230,8 +271,8 @@ end eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) - # TODO For variable smoothing_length use average smoothing length - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length^2) + smoothing_length_average = 0.5 * (smoothing_length_particle + smoothing_length_neighbor) + tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) volume_a = m_a / rho_a volume_b = m_b / rho_b @@ -250,7 +291,7 @@ end return visc .* v_diff end -function kinematic_viscosity(system, viscosity::ViscosityAdami) +function kinematic_viscosity(system, viscosity::ViscosityAdami, particle) return viscosity.nu end diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c10740ab8..33a1b7787 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -75,9 +75,9 @@ end @inline function density_diffusion_psi(::DensityDiffusionFerrari, rho_a, rho_b, pos_diff, distance, system, particle, neighbor) - (; smoothing_length) = system - - return ((rho_a - rho_b) / 2smoothing_length) * pos_diff / distance + # TODO should we use the average smoothing length here? + smoothing_length_particle = smoothing_length(system, particle) + return ((rho_a - rho_b) / (2 * smoothing_length_particle)) * pos_diff / distance end @doc raw""" @@ -207,7 +207,7 @@ end distance < sqrt(eps(typeof(distance))) && return (; delta) = density_diffusion - (; smoothing_length, state_equation) = particle_system + (; state_equation) = particle_system (; sound_speed) = state_equation volume_b = m_b / rho_b @@ -216,7 +216,10 @@ end particle_system, particle, neighbor) density_diffusion_term = dot(psi, grad_kernel) * volume_b - dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term + # TODO should we use the average smoothing length here? + smoothing_length_particle = smoothing_length(particle_system, particle) + dv[end, particle] += delta * smoothing_length_particle * sound_speed * + density_diffusion_term end # Density diffusion `nothing` or interaction other than fluid-fluid diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0fe58ab21..b874737bf 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -53,7 +53,8 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor) dv_pressure = pressure_correction * - pressure_acceleration(particle_system, neighbor_system, neighbor, + pressure_acceleration(particle_system, neighbor_system, + particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 6d66ba06a..3ea8352f7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -44,15 +44,14 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ -struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} +struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, + PF, ST, B, SRFT, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} density_calculator :: DC state_equation :: SE smoothing_kernel :: K - smoothing_length :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V density_diffusion :: DD @@ -62,6 +61,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, source_terms :: ST surface_tension :: SRFT buffer :: B + particle_refinement :: PR # TODO cache :: C end @@ -76,6 +76,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, + particle_refinement=nothing, surface_tension=nothing) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -116,14 +117,18 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = (; create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, + smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, nothing, - source_terms, surface_tension, buffer, cache) + source_terms, surface_tension, buffer, + particle_refinement, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4da0f7783..114f72aed 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -32,7 +32,7 @@ end rho_b = neighbor_system.material_density[neighbor] grad_kernel = smoothing_kernel_grad(particle_system, initial_pos_diff, - initial_distance) + initial_distance, particle) m_a = particle_system.mass[particle] m_b = neighbor_system.mass[neighbor] @@ -88,7 +88,7 @@ function interact!(dv, v_particle_system, u_particle_system, # Use kernel from the fluid system in order to get the same force here in # solid-fluid interaction as for fluid-solid interaction. - grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance, particle) # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles # corresponding to the chosen boundary model. @@ -99,7 +99,8 @@ function interact!(dv, v_particle_system, u_particle_system, # are switched in the following two calls. # This way, we obtain the exact same force as for the fluid-solid interaction, # but with a flipped sign (because `pos_diff` is flipped compared to fluid-solid). - dv_boundary = pressure_acceleration(neighbor_system, particle_system, particle, + dv_boundary = pressure_acceleration(neighbor_system, particle_system, + neighbor, particle, m_b, m_a, p_b, p_a, rho_b, rho_a, pos_diff, distance, grad_kernel, neighbor_system.correction) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index dfd71a2d6..8275da22b 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -286,7 +286,7 @@ end pos_diff = current_coords(system, particle) - current_coords(system, neighbor) grad_kernel = smoothing_kernel_grad(system, initial_pos_diff, - initial_distance) + initial_distance, particle) result = volume * pos_diff * grad_kernel' diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index ef5f68cc6..90c5fa435 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -252,7 +252,7 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["viscosity"] = type2string(system.viscosity) write2vtk!(vtk, system.viscosity) vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) - vtk["smoothing_length"] = system.smoothing_length + # vtk["smoothing_length"] = system.smoothing_length TODO vtk["density_calculator"] = type2string(system.density_calculator) if system isa WeaklyCompressibleSPHSystem