From 0974472035d5ac451ad57614c17b3fb3e917bf18 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 8 Nov 2024 15:20:16 +0100 Subject: [PATCH] fix --- src/TrixiParticles.jl | 1 + src/multi_resolution/multi_resolution.jl | 1 + src/multi_resolution/particle_refinement.jl | 19 ++++++++++++++----- src/multi_resolution/refinement_pattern.jl | 6 +++--- .../fluid/entropically_damped_sph/rhs.jl | 19 +++++++++++++++---- .../fluid/entropically_damped_sph/system.jl | 16 +++++++++++----- src/schemes/fluid/fluid.jl | 4 +++- src/schemes/fluid/transport_velocity.jl | 2 +- src/schemes/fluid/viscosity.jl | 12 +++++++----- 9 files changed, 56 insertions(+), 24 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 43be6e010..8eff154db 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -41,6 +41,7 @@ include("general/system.jl") # `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded` include("util.jl") include("preprocessing/preprocessing.jl") +include("multi_resolution/multi_resolution.jl") include("callbacks/callbacks.jl") include("general/general.jl") include("setups/setups.jl") diff --git a/src/multi_resolution/multi_resolution.jl b/src/multi_resolution/multi_resolution.jl index c7c0641ca..7a261a34c 100644 --- a/src/multi_resolution/multi_resolution.jl +++ b/src/multi_resolution/multi_resolution.jl @@ -1,3 +1,4 @@ +include("refinement_criteria.jl") include("particle_refinement.jl") include("refinement_pattern.jl") diff --git a/src/multi_resolution/particle_refinement.jl b/src/multi_resolution/particle_refinement.jl index 8fb10d32e..85471f0cd 100644 --- a/src/multi_resolution/particle_refinement.jl +++ b/src/multi_resolution/particle_refinement.jl @@ -19,7 +19,10 @@ function ParticleRefinement(; splitting_pattern, max_spacing_ratio, split_candidates = Vector{ELTYPE}() delete_candidates = Vector{Bool}() - return ParticleRefinement(splitting_pattern, Tuple(refinement_criteria), + if refinement_criteria isa Tuple + refinement_criteria = (refinement_criteria,) + end + return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, mass_ref, merge_candidates, split_candidates, delete_candidates, 0, 0) end @@ -45,12 +48,18 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) # Resize neighborhood search foreach_system(semi) do system - foreach_system(semi) do neighbor_system - search = get_neighborhood_search(system, neighbor_system, semi) - u_neighbor = wrap_u(u_ode, neighbor_system, semi) + foreach_system(semi) do neighbor + nhs_old = get_neighborhood_search(system, neighbor, semi) + u_neighbor = wrap_u(u_ode, neighbor, semi) # TODO - resize_nhs!(search, system, neighbor_system, u_neighbor) + resize_nhs!(nhs_old, system, neighbor, u_neighbor) + nhs = copy_neighborhood_search(nhs_old, compact_support(system, neighbor), + nparticles(neighbor)) + + system_coords = current_coordinates(u_ode, system) + neighbor_coords = current_coordinates(u_ode, neighbor) + PointNeighbors.initialize!(nhs, system_coords, neighbor_coords) end end diff --git a/src/multi_resolution/refinement_pattern.jl b/src/multi_resolution/refinement_pattern.jl index c3230572b..759fe8262 100644 --- a/src/multi_resolution/refinement_pattern.jl +++ b/src/multi_resolution/refinement_pattern.jl @@ -17,7 +17,7 @@ struct CubicSplitting{ELTYPE} SVector{2, ELTYPE}(epsilon * direction_3), SVector{2, ELTYPE}(epsilon * direction_4)] - if refinement_pattern.center_particle + if center_particle push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end @@ -45,7 +45,7 @@ struct TriangularSplitting{ELTYPE} SVector{2, ELTYPE}(epsilon * direction_2), SVector{2, ELTYPE}(epsilon * direction_3)] - if refinement_pattern.center_particle + if center_particle push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end @@ -80,7 +80,7 @@ struct HexagonalSplitting{ELTYPE} SVector{2, ELTYPE}(epsilon * direction_5), SVector{2, ELTYPE}(epsilon * direction_6)] - if refinement_pattern.center_particle + if center_particle push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index aadd8d28e..e8d6679b1 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -157,6 +157,8 @@ end function 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_a, rho_b) + (; sound_speed) = particle_system + # EDAC pressure evolution pressure_diff = p_a - p_b @@ -227,14 +229,23 @@ end return zero(pos_diff) end +#TODO @inline function velocity_correction(particle_system, neighbor_system, - particle_refinement, - pos_diff, distance, + particle_refinement, 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) +end + +@inline function velocity_correction(particle_system, neighbor_system::FluidSystem, + 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_a = current_velocity(v_particle_system, particle_system, particle) + advection_velocity_a = advection_velocity(v_particle_system, particle_system, particle) momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 15fa2941a..9b99b71f0 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -209,16 +209,22 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) compute_density!(system, u, u_ode, semi, system.density_calculator) update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi) - factor_for_variable_smoothing_length!(system, system.particle_refinement, u_ode, semi) + factor_for_variable_smoothing_length!(system, system.particle_refinement, + v_ode, u_ode, semi) end -factor_for_variable_smoothing_length!(system, ::Nothing, u_ode, semi) = system +factor_for_variable_smoothing_length!(system, ::Nothing, v_ode, u_ode, semi) = system -function factor_for_variable_smoothing_length!(system, particle_refinement, u_ode, semi) +function factor_for_variable_smoothing_length!(system, particle_refinement, + v_ode, u_ode, semi) (; smoothing_kernel) = system - u = wrap_u(u_ode, system, semi) + (; beta) = system.cache + NDIMS = ndims(system) + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, semi) + # Use all other systems for the average pressure foreach_system(semi) do neighbor_system u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) @@ -238,7 +244,7 @@ function factor_for_variable_smoothing_length!(system, particle_refinement, u_od h_a = smoothing_length(system, particle) w_deriv = kernel_deriv(smoothing_kernel, distance, h_a) - beta[partcle] += -m_a * distance * w_deriv / (NDIMS * rho_a) + beta[particle] += -m_a * distance * w_deriv / (NDIMS * rho_a) end end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index e2e659470..3a99ef10f 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -19,8 +19,10 @@ 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_lengt_factor=smoothng_length_factor) + smoothing_lengt_factor=smoothng_length_factor, + beta=zeros(nparticles(initial_condition))) end @inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index ef5a540be..8e9923a41 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -113,7 +113,7 @@ end 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_a = beta_correction(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)' diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 015033156..35dd0f75b 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -228,7 +228,7 @@ end beta_inv_a = beta_correction(particle_system, particle) nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system)) + viscosity_model(neighbor_system, particle_system), particle) grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) @@ -237,6 +237,10 @@ end tmp = (distance^2 + 0.001 * smoothing_length(particle_system, particle)^2) + 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 + return m_b * beta_inv_a * 4 * nu_a * dot(pos_diff, grad_W_avg) * v_diff / (tmp * (rho_a + rho_b)) end @@ -253,11 +257,9 @@ end # TODO is that correct with the smoothing lengths? nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system), - smoothing_length_particle) + viscosity_model(neighbor_system, particle_system), particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system), - smoothing_length_neighbor) + viscosity_model(particle_system, neighbor_system), particle) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)