From a78183d871a5396fae540babba719225211d19a0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 17:34:41 +0100 Subject: [PATCH 01/64] add callback --- src/callbacks/callbacks.jl | 1 + src/callbacks/refinement.jl | 123 ++++++++++++++++++++++++++++++++++++ 2 files changed, 124 insertions(+) create mode 100644 src/callbacks/refinement.jl 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..03f7f6fe7 --- /dev/null +++ b/src/callbacks/refinement.jl @@ -0,0 +1,123 @@ +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) + 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 From 152e0a5a09faaa118fb1afd58e06b7252d05e71f Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 12:38:05 +0100 Subject: [PATCH 02/64] First prototype of variable smoothing length --- src/general/corrections.jl | 6 +-- src/general/density_calculators.jl | 2 +- src/general/system.jl | 25 +++++------- src/schemes/boundary/rhs.jl | 2 +- .../fluid/entropically_damped_sph/rhs.jl | 9 ++--- src/schemes/fluid/surface_tension.jl | 6 +-- src/schemes/fluid/viscosity.jl | 38 ++++++++++++------- .../density_diffusion.jl | 13 ++++--- .../fluid/weakly_compressible_sph/system.jl | 31 +++++++++++++-- src/schemes/solid/total_lagrangian_sph/rhs.jl | 4 +- .../solid/total_lagrangian_sph/system.jl | 2 +- 11 files changed, 84 insertions(+), 54 deletions(-) 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/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/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..862471ed1 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -30,7 +30,7 @@ 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, m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, @@ -72,8 +72,6 @@ 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 - volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -88,8 +86,9 @@ end 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 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/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1a0665a71..aa8bc2c0f 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -123,13 +123,20 @@ end 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 +177,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 +220,18 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - (; smoothing_length) = particle_system - 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)) + 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 +242,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 +262,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/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 6d66ba06a..df1d65d5d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -45,14 +45,13 @@ 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} + 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,6 +117,9 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = (; create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) + cache = (; + create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, n_particles)..., + cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, @@ -123,7 +127,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, 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) = (;) @@ -155,6 +160,26 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) return (; surface_normal) end +function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles) + return (; smoothing_length) +end + +function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles) + return (; smoothing_length=smoothing_length * ones(n_particles)) +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4da0f7783..ed1f4f745 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. 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' From a7fc80b5105356b2168a78eb7dbcd97a54d7d8e5 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 03/64] Make it work with EDAC --- src/general/semidiscretization.jl | 8 ++-- .../dummy_particles/dummy_particles.jl | 6 ++- .../fluid/entropically_damped_sph/system.jl | 39 ++++++++++++------- .../fluid/weakly_compressible_sph/system.jl | 2 +- src/visualization/write2vtk.jl | 2 +- 5 files changed, 38 insertions(+), 19 deletions(-) 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/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 8b7e507ed..345e3c807 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 + diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..7b230ddc3 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,9 +58,11 @@ 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, +function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), @@ -69,6 +70,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, 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) @@ -98,17 +100,16 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, cache = create_cache_density(initial_condition, density_calculator) cache = (; create_cache_edac(initial_condition, transport_velocity)..., 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) + cache = (; + create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, + nparticles(initial_condition))..., + cache...) + + 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 -end function Base.show(io::IO, system::EntropicallyDampedSPHSystem) @nospecialize system # reduce precompilation time @@ -146,6 +147,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index df1d65d5d..5a4a0d103 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -123,7 +123,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, + smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, nothing, 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 From 5afab8cedc2813a8121e12ac3fdf517206f0bd69 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 04/64] add h-factor --- .../fluid/entropically_damped_sph/rhs.jl | 2 +- .../fluid/entropically_damped_sph/system.jl | 96 ++++++++----------- src/schemes/fluid/fluid.jl | 22 +++++ .../fluid/weakly_compressible_sph/system.jl | 28 +----- 4 files changed, 69 insertions(+), 79 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 862471ed1..71a639974 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -87,7 +87,7 @@ end eta_tilde = 2 * eta_a * eta_b / (eta_a + eta_b) smoothing_length_average = 0.5 * (smoothing_length(particle_system, particle) + - 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 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 7b230ddc3..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -63,53 +63,53 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, 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)), - 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) - - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - - mass = copy(initial_condition.mass) - - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end + 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) + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + mass = copy(initial_condition.mass) + + 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_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, - nparticles(initial_condition))..., - 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...) - 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 + 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) @nospecialize system # reduce precompilation time @@ -147,18 +147,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fd0a9c3a5..1b0413010 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 diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 5a4a0d103..3ea8352f7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -44,8 +44,8 @@ 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, PR, 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} @@ -118,8 +118,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) cache = (; - create_cache_refinement(particle_refinement, smoothing_length, ELTYPE, NDIMS, n_particles)..., - cache...) + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, @@ -160,26 +160,6 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) return (; surface_normal) end -function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles) - return (; smoothing_length) -end - -function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles) - return (; smoothing_length=smoothing_length * ones(n_particles)) -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time From a219d3c08290e86fa87640deea00b0b37b315ac4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 13:04:47 +0100 Subject: [PATCH 05/64] adapt rhs --- .../fluid/entropically_damped_sph/rhs.jl | 210 ++++++++++++++++-- src/schemes/fluid/fluid.jl | 7 +- src/schemes/fluid/pressure_acceleration.jl | 4 +- src/schemes/fluid/transport_velocity.jl | 25 ++- .../fluid/weakly_compressible_sph/rhs.jl | 3 +- src/schemes/solid/total_lagrangian_sph/rhs.jl | 3 +- 6 files changed, 230 insertions(+), 22 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 71a639974..481030ec9 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -32,7 +32,8 @@ function interact!(dv, v_particle_system, u_particle_system, 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,28 @@ 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_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,9 +78,45 @@ 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) +@inline function pressure_evolution!(dv, particle_system, 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_system.particle_refinement, + particle) + beta_inv_b = beta_correction(neighbor_system, neighbor_system.particle_refinement, + 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) + + 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 + +function beta_correction(particle_system, ::Nothing, particle) + return one(eltype(particle_system)) +end + +function beta_correction(particle_system, ::ParticleRefinement, 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_b, rho_b) volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -79,9 +124,6 @@ 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) @@ -101,11 +143,109 @@ 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, ::ParticleRefinement, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_b, rho_b) + # 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 * particle_system.smoothing_length[particle] / tmp + nu_edac_a = alpha_edac * sound_speed * particle_system.smoothing_length[neighbor] / tmp + + nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) + + # TODO: Use wrapped version + grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, + particle_system.smoothing_length[particle]) + grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, + neighbor_system.smoothing_length[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_damping_term(particle_system, neighbor_system, ::ParticleRefinement, + particle, neighbor, pos_diff, distance, 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, ::ParticleRefinement, + 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_refinement, particle) + beta_inv_b = beta_correction(particle_system, particle_refinement, 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 = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, + particle_system.smoothing_length[particle]) + grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, + neighbor_system.smoothing_length[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(eltype(particle_system)) +end + +@inline function velocity_correction(particle_system, neighbor_system, + particle_refinement::ParticleRefinement, + 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_refinement, 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,:]`. @@ -124,3 +264,45 @@ 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::ParticleRefinement, + 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_refinement, particle) * + (p_a - p_a_avg) / rho_a^2 + P_b = beta_correction(particle_system, particle_refinement, neighbor) * + (p_b - p_b_avg) / rho_b^2 + + # TODO: Use wrapped version + grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, + particle_system.smoothing_length[particle]) + grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, + neighbor_system.smoothing_length[neighbor]) + + return -m_b * (P_a * grad_kernel_a + P_b * grad_kernel_b) +end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 1b0413010..9b8896c3e 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -111,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) @@ -120,9 +120,10 @@ 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) + 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/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bc8fdf005..90b744f1d 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, + ::ParticleRefinement, 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_refinement, particle) / rho_a + factor_b = beta_correction(neighbor_system, particle_refinement, 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/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/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index ed1f4f745..114f72aed 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -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) From 847566ad92370ec0cab0b0e25c82ef2b89eb803c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 13:07:20 +0100 Subject: [PATCH 06/64] =?UTF-8?q?ada=C3=BCt=20viscosity?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../dummy_particles/dummy_particles.jl | 11 +-- .../fluid/entropically_damped_sph/rhs.jl | 81 ++++++++----------- src/schemes/fluid/fluid.jl | 5 +- src/schemes/fluid/transport_velocity.jl | 6 +- src/schemes/fluid/viscosity.jl | 41 ++++++++-- 5 files changed, 81 insertions(+), 63 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 345e3c807..ece912eeb 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -585,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/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 481030ec9..b075813f2 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -64,7 +64,8 @@ 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) - pressure_evolution!(dv, particle_system, neighbor_system, v_diff, grad_kernel, + 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) @@ -79,6 +80,7 @@ function interact!(dv, v_particle_system, u_particle_system, end @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) @@ -87,16 +89,14 @@ end # 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_system.particle_refinement, - particle) - beta_inv_b = beta_correction(neighbor_system, neighbor_system.particle_refinement, - neighbor) + 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) + + 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, @@ -105,18 +105,21 @@ end return dv end +@inline beta_correction(system, particle) = one(eltype(system)) -function beta_correction(particle_system, ::Nothing, particle) - return one(eltype(particle_system)) +@inline function beta_correction(system::FluidSystem, particle) + beta_correction(system, system.particle_refinement, particle) end -function beta_correction(particle_system, ::ParticleRefinement, particle) +@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_b, rho_b) + 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 @@ -146,9 +149,9 @@ function pressure_damping_term(particle_system, neighbor_system, ::Nothing, return volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) end -function pressure_damping_term(particle_system, neighbor_system, ::ParticleRefinement, +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_b, rho_b) + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel_a) # EDAC pressure evolution pressure_diff = p_a - p_b @@ -158,24 +161,17 @@ function pressure_damping_term(particle_system, neighbor_system, ::ParticleRefin # 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 * particle_system.smoothing_length[particle] / tmp - nu_edac_a = alpha_edac * sound_speed * particle_system.smoothing_length[neighbor] / tmp + 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) - # TODO: Use wrapped version - grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, - particle_system.smoothing_length[particle]) - grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, - neighbor_system.smoothing_length[neighbor]) + 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_damping_term(particle_system, neighbor_system, ::ParticleRefinement, - particle, neighbor, pos_diff, distance, m_b, rho_b) -end function pressure_reduction(particle_system, neighbor_system, ::Nothing, v_particle_system, v_neighbor_system, @@ -184,12 +180,12 @@ function pressure_reduction(particle_system, neighbor_system, ::Nothing, return zero(eltype(particle_system)) end -function pressure_reduction(particle_system, neighbor_system, ::ParticleRefinement, +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_refinement, particle) - beta_inv_b = beta_correction(particle_system, particle_refinement, neighbor) + 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) @@ -197,10 +193,8 @@ function pressure_reduction(particle_system, neighbor_system, ::ParticleRefineme 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 = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, - particle_system.smoothing_length[particle]) - grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, - neighbor_system.smoothing_length[neighbor]) + 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) @@ -224,12 +218,11 @@ end v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) - return zero(eltype(particle_system)) + return zero(pos_diff) end @inline function velocity_correction(particle_system, neighbor_system, - particle_refinement::ParticleRefinement, - 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) @@ -242,7 +235,7 @@ end 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_refinement, particle) + 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 @@ -286,23 +279,17 @@ end end @inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, - particle_refinement::ParticleRefinement, - neighbor_system, particle, neighbor, m_a, m_b, p_a, - p_b, - rho_a, rho_b, pos_diff, distance, W_a, correction) + 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_refinement, particle) * - (p_a - p_a_avg) / rho_a^2 - P_b = beta_correction(particle_system, particle_refinement, neighbor) * - (p_b - p_b_avg) / rho_b^2 + 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 - # TODO: Use wrapped version - grad_kernel_a = kernel_grad(particle_system.smoothing_kernel, pos_diff, distance, - particle_system.smoothing_length[particle]) - grad_kernel_b = kernel_grad(neighbor_system.smoothing_kernel, pos_diff, distance, - neighbor_system.smoothing_length[neighbor]) + 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/fluid.jl b/src/schemes/fluid/fluid.jl index 9b8896c3e..f287cabba 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -124,6 +124,7 @@ end 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, - pos_diff, distance, 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/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 90b744f1d..5270d1210 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -103,7 +103,7 @@ end end @inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, - ::ParticleRefinement, pos_diff, distance, + 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) @@ -112,8 +112,8 @@ 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_refinement, particle) / rho_a - factor_b = beta_correction(neighbor_system, particle_refinement, neighbor) / rho_b + 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)' diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index aa8bc2c0f..222fa6fd8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -112,12 +112,9 @@ 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) @@ -220,6 +217,38 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + 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) + + 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) From 9a668d765c2fee9da451cee04e6c01426c887444 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:06:13 +0100 Subject: [PATCH 07/64] add `ParticleRefinement` --- src/TrixiParticles.jl | 3 +- src/refinement/refinement.jl | 133 ++++++++++++++++++++++++++ src/refinement/refinement_criteria.jl | 64 +++++++++++++ 3 files changed, 199 insertions(+), 1 deletion(-) create mode 100644 src/refinement/refinement.jl create mode 100644 src/refinement/refinement_criteria.jl 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/refinement/refinement.jl b/src/refinement/refinement.jl new file mode 100644 index 000000000..ae435e481 --- /dev/null +++ b/src/refinement/refinement.jl @@ -0,0 +1,133 @@ +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, 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, u_ode) + foreach_system(semi) do system + u = wrap_u(u_ode, system, semi) + update_particle_spacing(system, u, semi) + end +end + +# The methods for the `FluidSystem` are defined in `src/schemes/fluid/fluid.jl` +@inline update_particle_spacing(system, u, semi) = systemerror + +@inline function update_particle_spacing(system::FluidSystem, u, semi) + update_particle_spacing(system, system.particle_refinement, u, semi) +end + +@inline update_particle_spacing(system, ::Nothing, u, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, particle_refinement, u, semi) + (; smoothing_length, smoothing_length_factor) = system.cache + (; mass_ref) = particle_refinement + + system_coords = current_coordinates(u, system) + + for particle in eachparticle(system) + dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, 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] = system.density[particle] * new_spacing^(ndims(system)) + end + + return system +end + +@inline function min_max_avg_spacing(system, semi, 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(particle_system, neighbor_system, + semi) + + u_neighbor_system = wrap_u(u_ode, neighbor, 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..de4bab43a --- /dev/null +++ b/src/refinement/refinement_criteria.jl @@ -0,0 +1,64 @@ +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, 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 From c89285ce12c8b1d5d41d695994f31fe59e8fd1c5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:08:11 +0100 Subject: [PATCH 08/64] initial refinement: resize --- src/callbacks/refinement.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/callbacks/refinement.jl b/src/callbacks/refinement.jl index 03f7f6fe7..130de8756 100644 --- a/src/callbacks/refinement.jl +++ b/src/callbacks/refinement.jl @@ -41,6 +41,12 @@ function initial_refinement!(cb, 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 From 6ff4fa68130904854894fac036eb62dab9473649 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:04:53 +0100 Subject: [PATCH 09/64] fix tuple check --- src/refinement/refinement.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index ae435e481..3cbcfa1d5 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -10,7 +10,7 @@ function ParticleRefinement(; max_spacing_ratio, refinement_criteria=SpatialRefinementCriterion()) mass_ref = Vector{eltype(max_spacing_ratio)}() - if refinement_criteria isa Tuple + if !(refinement_criteria isa Tuple) refinement_criteria = (refinement_criteria,) end return ParticleRefinement(refinement_criteria, max_spacing_ratio, mass_ref) From 6ff9f037e11d98a12c5ab85e8055e68ad8c0ee81 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:18:31 +0100 Subject: [PATCH 10/64] fix update --- src/refinement/refinement.jl | 37 +++++++++++++++------------ src/refinement/refinement_criteria.jl | 2 +- 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 3cbcfa1d5..19a58229d 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -34,7 +34,7 @@ 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, u_ode) + update_particle_spacing(semi, v_ode, u_ode) # Split the particles (Algorithm 2) @@ -51,30 +51,34 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) return semi end -function update_particle_spacing(semi, u_ode) +function update_particle_spacing(semi, v_ode, u_ode) foreach_system(semi) do system - u = wrap_u(u_ode, system, semi) - update_particle_spacing(system, u, semi) + 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, u, semi) = systemerror +@inline update_particle_spacing(system, v_ode, u_ode, semi) = system -@inline function update_particle_spacing(system::FluidSystem, u, semi) - update_particle_spacing(system, system.particle_refinement, u, semi) +@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, ::Nothing, u, semi) = system +@inline update_particle_spacing(system, ::Nothing, v_ode, u_ode, semi) = system -@inline function update_particle_spacing(system::FluidSystem, particle_refinement, u, semi) +@inline function update_particle_spacing(system::FluidSystem, particle_refinement, + v_ode, u_ode, semi) (; smoothing_length, smoothing_length_factor) = system.cache - (; mass_ref) = particle_refinement + (; 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, system_coords, particle) + 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) @@ -83,23 +87,24 @@ end end smoothing_length[particle] = smoothing_length_factor * new_spacing - mass_ref[particle] = system.density[particle] * new_spacing^(ndims(system)) + mass_ref[particle] = particle_density(v, system, particle) * + new_spacing^(ndims(system)) end return system end -@inline function min_max_avg_spacing(system, semi, system_coords, particle) + +@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(particle_system, neighbor_system, - semi) + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) - u_neighbor_system = wrap_u(u_ode, neighbor, 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, diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl index de4bab43a..bb8a95017 100644 --- a/src/refinement/refinement_criteria.jl +++ b/src/refinement/refinement_criteria.jl @@ -36,7 +36,7 @@ end return system end -@inline set_particle_spacing!(system, _, _, _, _) = system +@inline set_particle_spacing!(system, _, _, _, _, _) = system @inline function set_particle_spacing!(particle_system, neighbor_system::Union{BoundarySystem, SolidSystem}, From 39d6a9b721ab4f3f2ada70700fdee1a9849049cc Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 16:02:14 +0100 Subject: [PATCH 11/64] fix dispatch --- src/refinement/refinement_criteria.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl index bb8a95017..48c3c8b5c 100644 --- a/src/refinement/refinement_criteria.jl +++ b/src/refinement/refinement_criteria.jl @@ -19,6 +19,11 @@ end check_refinement_criteria!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) end +@inlin 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 From 9a4daa4cc8e77b1af905f045f62c23a59bbc3c04 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 16:06:51 +0100 Subject: [PATCH 12/64] fix typo --- src/refinement/refinement_criteria.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl index 48c3c8b5c..627ffb619 100644 --- a/src/refinement/refinement_criteria.jl +++ b/src/refinement/refinement_criteria.jl @@ -19,7 +19,7 @@ end check_refinement_criteria!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) end -@inlin function check_refinement_criteria!(system::FluidSystem, ::Nothing, +@inline function check_refinement_criteria!(system::FluidSystem, ::Nothing, v, u, v_ode, u_ode, semi) system end From 486ba70ba45e5ae0c2c62c79e711141f56f45b1b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 16:10:46 +0100 Subject: [PATCH 13/64] fix ambiguous method --- src/refinement/refinement.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 19a58229d..fe1f90b10 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -64,7 +64,7 @@ end update_particle_spacing(system, system.particle_refinement, v_ode, u_ode, semi) end -@inline update_particle_spacing(system, ::Nothing, v_ode, u_ode, semi) = system +@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) From 5a209e77e6836e3f6b92c1df679b5cdd5a52a5c7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Nov 2024 16:38:58 +0100 Subject: [PATCH 14/64] add resize --- src/TrixiParticles.jl | 1 + src/multi_resolution/multi_resolution.jl | 1 + src/multi_resolution/resize.jl | 170 +++++++++++++++++++++++ 3 files changed, 172 insertions(+) create mode 100644 src/multi_resolution/multi_resolution.jl create mode 100644 src/multi_resolution/resize.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 8193ec9cf..9f77ffa82 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 new file mode 100644 index 000000000..648a9fb1f --- /dev/null +++ b/src/multi_resolution/multi_resolution.jl @@ -0,0 +1 @@ +include("resize.jl") diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl new file mode 100644 index 000000000..23203e5c4 --- /dev/null +++ b/src/multi_resolution/resize.jl @@ -0,0 +1,170 @@ +function resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Resize all systems + foreach_system(semi) do system + resize!(system, capacity(system)) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Delete at specific indices + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + deleteat!(system, v, u) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) + copyto!(_v_ode, v_ode) + copyto!(_u_ode, u_ode) + + # Get ranges after resizing the systems + ranges_v_new, ranges_u_new = ranges_vu(semi.systems) + + ranges_v_old = semi.ranges_v + ranges_u_old = semi.ranges_u + + # Set ranges after resizing the systems + for i in 1:length(semi.systems) + semi.ranges_v[i] = ranges_v_new[i] + semi.ranges_u[i] = ranges_u_new[i] + end + + for i in eachindex(ranges_u_old) + length_u = min(length(ranges_u_old[i]), length(ranges_u_new[i])) + for j in 1:length_u + u_ode[ranges_u_new[i][1] + j] = _u_ode[ranges_u_old[i][1] + j] + end + + length_v = min(length(ranges_v_old[i]), length(ranges_v_new[i])) + for j in 1:length_v + v_ode[ranges_v_new[i][1] + j] = _v_ode[ranges_v_old[i][1] + j] + end + end + + capacity_global = sum(system -> nparticles(system), semi.systems) + + resize!(v_ode, capacity_global) + resize!(u_ode, capacity_global) + + resize!(_v_ode, capacity_global) + resize!(_u_ode, capacity_global) + + # TODO: Do the following in the callback + # resize!(integrator, (length(v_ode), length(u_ode))) + + # # Tell OrdinaryDiffEq that u has been modified + # u_modified!(integrator, true) + return v_ode +end + +resize!(system, capacity_system) = system + +function resize!(system::FluidSystem, capacity_system) + return resize!(system, system.particle_refinement, capacity_system) +end + +resize!(system, ::Nothing, capacity_system) = system + +function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) + (; mass, pressure, cache, density_calculator) = system + + refinement.n_particles_before_resize = nparticles(system) + + resize!(mass, capacity_system) + resize!(pressure, capacity_system) + resize_density!(system, capacity_system, density_calculator) + resize_cache!(system, cache, n) +end + +function resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) + (; mass, cache, density_calculator) = system + + refinement.n_particles_before_resize = nparticles(system) + + resize!(mass, capacity_system) + resize_density!(system, capacity_system, density_calculator) + resize_cache!(system, capacity_system) + + return system +end + +resize_density!(system, n::Int, ::SummationDensity) = resize!(system.cache.density, n) +resize_density!(system, n::Int, ::ContinuityDensity) = system + +function resize_cache!(system, n::Int) + resize!(system.cache.smoothing_length, n) + + return system +end + +function resize_cache!(system::EntropicallyDampedSPHSystem, n) + resize!(system.cache.smoothing_length, n) + resize!(system.cache.pressure_average, n) + resize!(system.cache.neighbor_counter, n) + + return system +end + +deleteat!(system, v, u) = system + +function deleteat!(system::FluidSystem, v, u) + return deleteat!(system, system.particle_refinement, v, u) +end + +deleteat!(system, ::Nothing, v, u) = system + +function deleteat!(system::FluidSystem, refinement, v, u) + (; delete_candidates) = refinement + + delete_counter = 0 + + for particle in eachparticle(system) + if !iszero(delete_candidates[particle]) + # swap particles (keep -> delete) + dump_id = nparticles(system) - delete_counter + + vel_keep = current_velocity(v, system, dump_id) + pos_keep = current_coords(u, system, dump_id) + + mass_keep = hydrodynamic_mass(system, dump_id) + density_keep = particle_density(system, v, dump_id) + pressure_keep = particle_pressure(system, v, dump_id) + smoothing_length_keep = smoothing_length(system, dump_id) + + system.mass[particle] = mass_keep + system.cache.smoothing_length[particle] = smoothing_length_keep + + set_particle_pressure!(v, system, particle, pressure_keep) + + set_particle_density!(v, system, particle, density_keep) + + for dim in 1:ndims(system) + v[dim, particle] = vel_keep[dim] + u[dim, particle] = pos_keep[dim] + end + + delete_counter += 1 + end + end + + resize!(system, nparticles(system) - delete_counter) + + return system +end + +@inline capacity(system) = capacity(system, system.particle_refinement) + +@inline capacity(system, ::Nothing) = nparticles(system) + +@inline function capacity(system, particle_refinement) + return particle_refinement.n_new_particles + nparticles(system) +end From cba5e8fa7102f87e56a4c93307afb5b57d18c3ea Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Nov 2024 17:15:23 +0100 Subject: [PATCH 15/64] adapt for tests --- src/TrixiParticles.jl | 2 +- src/multi_resolution/multi_resolution.jl | 1 + src/multi_resolution/particle_refinement.jl | 8 ++++++++ src/multi_resolution/resize.jl | 12 +++++++----- 4 files changed, 17 insertions(+), 6 deletions(-) create mode 100644 src/multi_resolution/particle_refinement.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 9f77ffa82..a23708c7f 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -41,7 +41,6 @@ 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") @@ -54,6 +53,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") +include("multi_resolution/multi_resolution.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/multi_resolution/multi_resolution.jl b/src/multi_resolution/multi_resolution.jl index 648a9fb1f..547ee3c35 100644 --- a/src/multi_resolution/multi_resolution.jl +++ b/src/multi_resolution/multi_resolution.jl @@ -1 +1,2 @@ include("resize.jl") +include("particle_refinement.jl") diff --git a/src/multi_resolution/particle_refinement.jl b/src/multi_resolution/particle_refinement.jl new file mode 100644 index 000000000..b5d6f8340 --- /dev/null +++ b/src/multi_resolution/particle_refinement.jl @@ -0,0 +1,8 @@ +struct ParticleRefinement + n_particles_before_resize :: Ref{Int} + n_new_particles :: Ref{Int} +end + +function ParticleRefinement() + return ParticleRefinement(Ref(0), Ref(0)) +end diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl index 23203e5c4..a32556c30 100644 --- a/src/multi_resolution/resize.jl +++ b/src/multi_resolution/resize.jl @@ -77,22 +77,24 @@ resize!(system, ::Nothing, capacity_system) = system function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) (; mass, pressure, cache, density_calculator) = system - refinement.n_particles_before_resize = nparticles(system) + refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) resize!(pressure, capacity_system) resize_density!(system, capacity_system, density_calculator) - resize_cache!(system, cache, n) + # TODO + # resize_cache!(system, cache, n) end function resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) (; mass, cache, density_calculator) = system - refinement.n_particles_before_resize = nparticles(system) + refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) resize_density!(system, capacity_system, density_calculator) - resize_cache!(system, capacity_system) + # TODO + # resize_cache!(system, capacity_system) return system end @@ -166,5 +168,5 @@ end @inline capacity(system, ::Nothing) = nparticles(system) @inline function capacity(system, particle_refinement) - return particle_refinement.n_new_particles + nparticles(system) + return particle_refinement.n_new_particles[] + nparticles(system) end From 9d7f7286b5d760d6e18da35c8f8df5493e066f04 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 15 Nov 2024 14:44:25 +0100 Subject: [PATCH 16/64] fix `resize!(system, n)` --- src/general/semidiscretization.jl | 19 +++++++++++-------- src/multi_resolution/resize.jl | 26 ++++++++++++++------------ 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 6c0c7e37b..f734a6191 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -73,14 +73,7 @@ function Semidiscretization(systems...; # Other checks might be added here later. check_configuration(systems) - sizes_u = [u_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) - for i in eachindex(sizes_u)) - sizes_v = [v_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) - for i in eachindex(sizes_v)) + ranges_v, ranges_u = ranges_vu(systems) # Create a tuple of n neighborhood searches for each of the n systems. # We will need one neighborhood search for each pair of systems. @@ -92,6 +85,16 @@ function Semidiscretization(systems...; return Semidiscretization(systems, ranges_u, ranges_v, searches) end +function ranges_vu(systems) + sizes_u = [u_nvariables(system) * n_moving_particles(system) for system in systems] + ranges_u = [(sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) for i in eachindex(sizes_u)] + + sizes_v = [v_nvariables(system) * n_moving_particles(system) for system in systems] + ranges_v = [(sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)] + + return ranges_v, ranges_u +end + # Inline show function e.g. Semidiscretization(neighborhood_search=...) function Base.show(io::IO, semi::Semidiscretization) @nospecialize semi # reduce precompilation time diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl index a32556c30..2c7fa9405 100644 --- a/src/multi_resolution/resize.jl +++ b/src/multi_resolution/resize.jl @@ -1,4 +1,4 @@ -function resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) +function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) # Resize all systems foreach_system(semi) do system resize!(system, capacity(system)) @@ -22,15 +22,15 @@ function deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) return semi end -function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) +function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) copyto!(_v_ode, v_ode) copyto!(_u_ode, u_ode) # Get ranges after resizing the systems ranges_v_new, ranges_u_new = ranges_vu(semi.systems) - ranges_v_old = semi.ranges_v - ranges_u_old = semi.ranges_u + ranges_v_old = copy(semi.ranges_v) + ranges_u_old = copy(semi.ranges_u) # Set ranges after resizing the systems for i in 1:length(semi.systems) @@ -40,12 +40,12 @@ function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) for i in eachindex(ranges_u_old) length_u = min(length(ranges_u_old[i]), length(ranges_u_new[i])) - for j in 1:length_u + for j in 0:(length_u - 1) u_ode[ranges_u_new[i][1] + j] = _u_ode[ranges_u_old[i][1] + j] end length_v = min(length(ranges_v_old[i]), length(ranges_v_new[i])) - for j in 1:length_v + for j in 0:(length_v - 1) v_ode[ranges_v_new[i][1] + j] = _v_ode[ranges_v_old[i][1] + j] end end @@ -66,15 +66,15 @@ function resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) return v_ode end -resize!(system, capacity_system) = system +Base.resize!(system::System, capacity_system) = system -function resize!(system::FluidSystem, capacity_system) +function Base.resize!(system::FluidSystem, capacity_system) return resize!(system, system.particle_refinement, capacity_system) end -resize!(system, ::Nothing, capacity_system) = system +Base.resize!(system, ::Nothing, capacity_system) = system -function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) +function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) (; mass, pressure, cache, density_calculator) = system refinement.n_particles_before_resize[] = nparticles(system) @@ -86,7 +86,7 @@ function resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_syste # resize_cache!(system, cache, n) end -function resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) +function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) (; mass, cache, density_calculator) = system refinement.n_particles_before_resize[] = nparticles(system) @@ -163,7 +163,9 @@ function deleteat!(system::FluidSystem, refinement, v, u) return system end -@inline capacity(system) = capacity(system, system.particle_refinement) +@inline capacity(system) = nparticles(system) + +@inline capacity(system::FluidSystem) = capacity(system, system.particle_refinement) @inline capacity(system, ::Nothing) = nparticles(system) From 0acad857a81c2f0b4a5ba6273c7c9fe56affdb29 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 15 Nov 2024 14:56:15 +0100 Subject: [PATCH 17/64] fix `deleteat!` --- src/multi_resolution/particle_refinement.jl | 3 ++- src/multi_resolution/resize.jl | 19 ++++++++++--------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/multi_resolution/particle_refinement.jl b/src/multi_resolution/particle_refinement.jl index b5d6f8340..bf886ac12 100644 --- a/src/multi_resolution/particle_refinement.jl +++ b/src/multi_resolution/particle_refinement.jl @@ -1,8 +1,9 @@ struct ParticleRefinement n_particles_before_resize :: Ref{Int} n_new_particles :: Ref{Int} + delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles end function ParticleRefinement() - return ParticleRefinement(Ref(0), Ref(0)) + return ParticleRefinement(Ref(0), Ref(0), Bool[]) end diff --git a/src/multi_resolution/resize.jl b/src/multi_resolution/resize.jl index 2c7fa9405..87dc97fa5 100644 --- a/src/multi_resolution/resize.jl +++ b/src/multi_resolution/resize.jl @@ -9,7 +9,7 @@ function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) return semi end -function deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) +function Base.deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) # Delete at specific indices foreach_system(semi) do system v = wrap_v(v_ode, system, semi) @@ -116,15 +116,15 @@ function resize_cache!(system::EntropicallyDampedSPHSystem, n) return system end -deleteat!(system, v, u) = system +Base.deleteat!(system::System, v, u) = system -function deleteat!(system::FluidSystem, v, u) +function Base.deleteat!(system::FluidSystem, v, u) return deleteat!(system, system.particle_refinement, v, u) end -deleteat!(system, ::Nothing, v, u) = system +Base.deleteat!(system::FluidSystem, ::Nothing, v, u) = system -function deleteat!(system::FluidSystem, refinement, v, u) +function Base.deleteat!(system::FluidSystem, refinement, v, u) (; delete_candidates) = refinement delete_counter = 0 @@ -138,12 +138,13 @@ function deleteat!(system::FluidSystem, refinement, v, u) pos_keep = current_coords(u, system, dump_id) mass_keep = hydrodynamic_mass(system, dump_id) - density_keep = particle_density(system, v, dump_id) - pressure_keep = particle_pressure(system, v, dump_id) - smoothing_length_keep = smoothing_length(system, dump_id) + density_keep = particle_density(v, system, dump_id) + pressure_keep = particle_pressure(v, system,dump_id) + #TODO + # smoothing_length_keep = smoothing_length(system, dump_id) + # system.cache.smoothing_length[particle] = smoothing_length_keep system.mass[particle] = mass_keep - system.cache.smoothing_length[particle] = smoothing_length_keep set_particle_pressure!(v, system, particle, pressure_keep) From bd9b14b2ca4357e6c7f12a06d1161b1a10a495b1 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:24:47 +0100 Subject: [PATCH 18/64] restructure --- src/TrixiParticles.jl | 2 +- src/multi_resolution/multi_resolution.jl | 2 -- src/multi_resolution/particle_refinement.jl | 9 --------- src/{multi_resolution => refinement}/resize.jl | 0 4 files changed, 1 insertion(+), 12 deletions(-) delete mode 100644 src/multi_resolution/multi_resolution.jl delete mode 100644 src/multi_resolution/particle_refinement.jl rename src/{multi_resolution => refinement}/resize.jl (100%) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index a23708c7f..e445e43c3 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -53,7 +53,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") -include("multi_resolution/multi_resolution.jl") +include("refinement/refinement.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/multi_resolution/multi_resolution.jl b/src/multi_resolution/multi_resolution.jl deleted file mode 100644 index 547ee3c35..000000000 --- a/src/multi_resolution/multi_resolution.jl +++ /dev/null @@ -1,2 +0,0 @@ -include("resize.jl") -include("particle_refinement.jl") diff --git a/src/multi_resolution/particle_refinement.jl b/src/multi_resolution/particle_refinement.jl deleted file mode 100644 index bf886ac12..000000000 --- a/src/multi_resolution/particle_refinement.jl +++ /dev/null @@ -1,9 +0,0 @@ -struct ParticleRefinement - n_particles_before_resize :: Ref{Int} - n_new_particles :: Ref{Int} - delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles -end - -function ParticleRefinement() - return ParticleRefinement(Ref(0), Ref(0), Bool[]) -end diff --git a/src/multi_resolution/resize.jl b/src/refinement/resize.jl similarity index 100% rename from src/multi_resolution/resize.jl rename to src/refinement/resize.jl From d443541042be37f2308f1437443b2318d670d28d Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:33:43 +0100 Subject: [PATCH 19/64] include order --- src/TrixiParticles.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index e445e43c3..d87a03796 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -53,7 +53,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") -include("refinement/refinement.jl") +include("refinement/resize.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition From 131ab25b2825fe29967a2ccfad6d1e12cca9ea11 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:38:40 +0100 Subject: [PATCH 20/64] apply formatter --- src/refinement/resize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index 87dc97fa5..a7e70fa47 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -139,7 +139,7 @@ function Base.deleteat!(system::FluidSystem, refinement, v, u) mass_keep = hydrodynamic_mass(system, dump_id) density_keep = particle_density(v, system, dump_id) - pressure_keep = particle_pressure(v, system,dump_id) + pressure_keep = particle_pressure(v, system, dump_id) #TODO # smoothing_length_keep = smoothing_length(system, dump_id) # system.cache.smoothing_length[particle] = smoothing_length_keep From bf1be31f52c77fd26357fa523d5e6c813cd48b6b Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 21/64] Make it work with EDAC --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 75c0066c3..f4dc04806 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,6 +147,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 46906d8375603a90af6fa266eb2ae495b4b5b3c4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 22/64] add h-factor --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f4dc04806..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,18 +147,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 9140847c4a678a83b3da632cef0d8b97ad5a2dfa Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:12:06 +0100 Subject: [PATCH 23/64] add split pattern --- src/refinement/refinement_pattern.jl | 93 ++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 src/refinement/refinement_pattern.jl diff --git a/src/refinement/refinement_pattern.jl b/src/refinement/refinement_pattern.jl new file mode 100644 index 000000000..759fe8262 --- /dev/null +++ b/src/refinement/refinement_pattern.jl @@ -0,0 +1,93 @@ +struct CubicSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + function CubicSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) + ELTYPE = typeof(epsilon) + + direction_1 = [1 / sqrt(2), 1 / sqrt(2)] + direction_2 = [1 / sqrt(2), -1 / sqrt(2)] + direction_3 = -direction_1 + direction_4 = -direction_2 + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3), + SVector{2, ELTYPE}(epsilon * direction_4)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + end +end + +@inline function nchilds(system, refinement_pattern::CubicSplitting) + return 4 + refinement_pattern.center_particle +end + +struct TriangularSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + function TriangularSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) + ELTYPE = typeof(epsilon) + direction_1 = [0.0, 1.0] + direction_2 = [-sqrt(3) / 2, -0.5] + direction_3 = [sqrt(3) / 2, -0.5] + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + end +end + +@inline function nchilds(system::System{2}, refinement_pattern::TriangularSplitting) + return 3 + refinement_pattern.center_particle +end + +struct HexagonalSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + function HexagonalSplitting(; epsilon=0.4, alpha=0.9, center_particle=true) + ELTYPE = typeof(epsilon) + + direction_1 = [1.0, 0.0] + direction_2 = [-1.0, 0.0] + direction_3 = [0.5, sqrt(3) / 2] + direction_4 = [0.5, -sqrt(3) / 2] + direction_5 = [-0.5, sqrt(3) / 2] + direction_6 = [-0.5, -sqrt(3) / 2] + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3), + SVector{2, ELTYPE}(epsilon * direction_4), + SVector{2, ELTYPE}(epsilon * direction_5), + SVector{2, ELTYPE}(epsilon * direction_6)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + end +end + +@inline function nchilds(system::System{2}, refinement_pattern::HexagonalSplitting) + return 6 + refinement_pattern.center_particle +end From 7af99f4444758c8729962af9d92b30759bdd6bb7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:19:32 +0100 Subject: [PATCH 24/64] add split method --- src/refinement/refinement.jl | 18 ++++-- src/refinement/split.jl | 103 +++++++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 6 deletions(-) create mode 100644 src/refinement/split.jl diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index fe1f90b10..44b2482c3 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,19 +1,24 @@ include("refinement_criteria.jl") -struct ParticleRefinement{RC, ELTYPE} - refinement_criteria :: RC - max_spacing_ratio :: ELTYPE - mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles +struct ParticleRefinement{SP, RC, ELTYPE} + splitting_pattern :: SP + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + n_particles_before_resize :: Int + n_new_particles :: Int end -function ParticleRefinement(; max_spacing_ratio, +function ParticleRefinement(; splitting_pattern, 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) + + return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, + mass_ref, 0, 0) end resize_refinement!(system) = system @@ -37,6 +42,7 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) update_particle_spacing(semi, v_ode, u_ode) # Split the particles (Algorithm 2) + split_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) # Merge the particles (Algorithm 3) diff --git a/src/refinement/split.jl b/src/refinement/split.jl new file mode 100644 index 000000000..a04000b4e --- /dev/null +++ b/src/refinement/split.jl @@ -0,0 +1,103 @@ +function split_particles!(semi, v_ode, u_ode, v_ode_tmp, u_ode_tmp) + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + collect_split_candidates!(system, v, u) + end + + resize!(semi, v_ode, u_ode, v_ode_tmp, u_ode_tmp) + + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + split_particles!(system, v, u) + end + + return semi +end + +@inline collect_split_candidates!(system, v, u) = System + +@inline function collect_split_candidates!(system::FluidSystem, v, u) + return collect_split_candidates!(system, system.particle_refinement, v, u) +end + +@inline collect_split_candidates!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function collect_split_candidates!(system::FluidSystem, particle_refinement, v, u) + (; mass_ref, max_spacing_ratio, refinement_pattern) = particle_refinement + + n_split_candidates = 0 + + for particle in eachparticle(system) + m_a = hydrodynamic_mass(system, particle) + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a > m_max + n_split_candidates += 1 + end + end + + n_childs_exclude_center = nchilds(system, refinement_pattern) - 1 + + particle_refinement.n_new_particles[] = n_split_candidates * n_childs_exclude_center + + return system +end + +@inline split_particles!(system, v, u) = System + +@inline function split_particles!(system::FluidSystem, v, u) + return split_particles!(system, system.particle_refinement, v, u) +end + +@inline split_particles!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function split_particles!(system::FluidSystem, particle_refinement, v, u) + (; smoothing_length) = system.cache + (; mass_ref, max_spacing_ratio, refinement_pattern, n_particles_before_resize) = particle_refinement + (; alpha, relative_position) = refinement_pattern + + for particle in eachparticle(system) + m_a = hydrodynamic_mass(system, particle) + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a > m_max + smoothing_length_old = smoothing_length[particle] + mass_old = system.mass[particle] + + system.mass[particle] = mass_old / nchilds(system, refinement_pattern) + + set_particle_pressure!(v, system, partice, pressure) + + set_particle_density!(v, system, partice, density) + + smoothing_length[particle] = alpha * smoothing_length_old + + pos_center = current_coords(system, u, particle) + vel_center = current_velocity(system, v, particle) + + for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) + child = n_particles_before_resize + child_id_local + + system.mass[child] = mass_old / nchilds(system, refinement_pattern) + + set_particle_pressure!(v, system, child, pressure) + + set_particle_density!(v, system, child, density) + + smoothing_length[child] = alpha * smoothing_length_old + + rel_pos = smoothing_length_old * relative_position[child] + new_pos = pos_center + rel_pos + + for dim in 1:ndims(system) + u[dim, child] = new_pos[dim] + v[dim, child] = vel_center[dim] + end + end + end + end + + return system +end From df0554bbdc061f92a69cd9f7ff95fc4e6aca605e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:27:56 +0100 Subject: [PATCH 25/64] include --- src/TrixiParticles.jl | 3 ++- src/refinement/refinement.jl | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d87a03796..54bfb7957 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -22,7 +22,8 @@ using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, - get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem + get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, + RecursiveArrayTools @reexport using StaticArrays: SVector using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 44b2482c3..6577fd7a6 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,4 +1,6 @@ include("refinement_criteria.jl") +include("refinement_pattern.jl") +include("split.jl") struct ParticleRefinement{SP, RC, ELTYPE} splitting_pattern :: SP From 1a5dbe73e768a7265f493323546af043cf15a467 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:44:03 +0100 Subject: [PATCH 26/64] add merge candidates --- src/refinement/merge.jl | 113 +++++++++++++++++++++++++++++++++++ src/refinement/refinement.jl | 7 ++- 2 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 src/refinement/merge.jl diff --git a/src/refinement/merge.jl b/src/refinement/merge.jl new file mode 100644 index 000000000..0a25d824f --- /dev/null +++ b/src/refinement/merge.jl @@ -0,0 +1,113 @@ +function merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) + foreach_system(semi) do system + (; delete_candidates) = system.particle_refinement + delete_candidates .= false + + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + for _ in 1:3 + merge_particles!(system, v, u) + end + end + + deleteat!(semi, v_ode, u_ode, v_tmp, u_tmp) + + return semi +end + +@inline merge_particles!(system, v, u) = System + +@inline function merge_particles!(system::FluidSystem, v, u) + return merge_particles!(system, system.particle_refinement, v, u) +end + +@inline merge_particles!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function merge_particles!(system::FluidSystem, particle_refinement, v, u) + (; smoothing_kernel, cache) = system + (; mass_ref, max_spacing_ratio, merge_candidates, delete_candidates) = particle_refinement + + set_zero!(merge_candidates) + + system_coords = current_coordinates(u, particle_system) + neighborhood_search = get_neighborhood_search(system, semi) + + # Collect merge candidates + foreach_point_neighbor(system, system, system_coords, system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + particle == neighbor && return + delete_candidates[particle] && return + + m_a = hydrodynamic_mass(system, particle) + m_b = hydrodynamic_mass(system, neighbor) + + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a <= m_max + m_merge = m_a + m_b + m_max_min = min(m_max, max_spacing_ratio * mass_ref[neighbor]) + if m_merge < m_max_min + if merge_candidates[particle] == 0 + merge_candidates[particle] = neighbor + else + stored_neighbor = current_coords(u, system, merge_candidates[particle]) + pos_diff_stored = stored_neighbor - current_coords(u, system, particle) + + if distance < norm(pos_diff_stored) + merge_candidates[particle] = neighbor + end + end + end + end + end + + # Merge and delete particles + for particle in merge_candidates + iszero(merge_candidates[particle]) && continue + + candidate = merge_candidates[particle] + + if particle == merge_candidates[candidate] + if particle < candidate + m_a = hydrodynamic_mass(system, particle) + m_b = hydrodynamic_mass(system, candidate) + + m_merge = m_a + m_b + + pos_a = current_coords(u, system, particle) + pos_b = current_coords(u, system, candidate) + + vel_a = current_velocity(v, system, particle) + vel_b = current_velocity(v, system, candidate) + + pos_merge = (m_a * pos_a + m_b * pos_b) / m_merge + vel_merge = (m_a * vel_a + m_b * vel_b) / m_merge + + for dim in 1:ndims(system) + u[dim, particle] = pos_merge[dim] + v[dim, particle] = vel_merge[dim] + end + + # update smoothing length + h_a = smoothing_length(system, particle) + h_b = smoothing_length(system, candidate) + + # TODO: + # Check normalization_factor = smoothing_kernel(zero(pos_diff), one(h_a)) + tmp_m = m_merge * normalization_factor(smoothing_kernel, h_a) + tmp_a = m_a * smoothing_kernel(pos_merge - pos_a, h_a) + tmp_b = m_b * smoothing_kernel(pos_merge - pos_b, h_b) + + cache.smoothing_length[particle] = (tmp_m / (tmp_a + tmp_b))^(1 / + ndims(system)) + + system.mass[particl] = m_merge + else + delete_candidates[candidate] = true + end + end + end + + return system +end diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 6577fd7a6..b29807125 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -1,12 +1,14 @@ include("refinement_criteria.jl") include("refinement_pattern.jl") include("split.jl") +include("merge.jl") struct ParticleRefinement{SP, RC, ELTYPE} splitting_pattern :: SP refinement_criteria :: RC max_spacing_ratio :: ELTYPE mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles n_particles_before_resize :: Int n_new_particles :: Int end @@ -14,13 +16,14 @@ end function ParticleRefinement(; splitting_pattern, max_spacing_ratio, refinement_criteria=SpatialRefinementCriterion()) mass_ref = Vector{eltype(max_spacing_ratio)}() + delete_candidates = Vector{Bool}() if !(refinement_criteria isa Tuple) refinement_criteria = (refinement_criteria,) end return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, 0, 0) + mass_ref, delete_candidates, 0, 0) end resize_refinement!(system) = system @@ -33,6 +36,7 @@ resize_refinement!(system, ::Nothing) = system function resize_refinement!(system, particle_refinement) resize!(particle_refinement.mass_ref, nparticles(system)) + resize!(particle_refinement.delete_candidates, nparticles(system)) return system end @@ -47,6 +51,7 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) split_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) # Merge the particles (Algorithm 3) + merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) # Shift the particles From cc1921e85ba3b7c6bc28ec784c3f7937f0563439 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:02:48 +0100 Subject: [PATCH 27/64] add vector for merge candidates --- src/refinement/refinement.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index b29807125..dc01d85e2 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -8,6 +8,7 @@ struct ParticleRefinement{SP, RC, ELTYPE} refinement_criteria :: RC max_spacing_ratio :: ELTYPE mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles n_particles_before_resize :: Int n_new_particles :: Int @@ -23,7 +24,7 @@ function ParticleRefinement(; splitting_pattern, max_spacing_ratio, end return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, delete_candidates, 0, 0) + mass_ref, Int[], delete_candidates, 0, 0) end resize_refinement!(system) = system @@ -37,6 +38,7 @@ resize_refinement!(system, ::Nothing) = system function resize_refinement!(system, particle_refinement) resize!(particle_refinement.mass_ref, nparticles(system)) resize!(particle_refinement.delete_candidates, nparticles(system)) + resize!(particle_refinement.merge_candidates, nparticles(system)) return system end From 7b50acdce5cd3d6600c38258a8430f166a57ef16 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:22:44 +0100 Subject: [PATCH 28/64] fix --- src/refinement/merge.jl | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/refinement/merge.jl b/src/refinement/merge.jl index 0a25d824f..99b44c562 100644 --- a/src/refinement/merge.jl +++ b/src/refinement/merge.jl @@ -1,14 +1,9 @@ function merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) foreach_system(semi) do system - (; delete_candidates) = system.particle_refinement - delete_candidates .= false - v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - for _ in 1:3 - merge_particles!(system, v, u) - end + merge_particles!(system, semi, v, u) end deleteat!(semi, v_ode, u_ode, v_tmp, u_tmp) @@ -16,21 +11,34 @@ function merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) return semi end -@inline merge_particles!(system, v, u) = System +@inline merge_particles!(system, semi, v, u) = system -@inline function merge_particles!(system::FluidSystem, v, u) - return merge_particles!(system, system.particle_refinement, v, u) +@inline function merge_particles!(system::FluidSystem, semi, v, u) + return merge_particles!(system, system.particle_refinement, semi, v, u) end -@inline merge_particles!(system::FluidSystem, ::Nothing, v, u) = system +@inline merge_particles!(system::FluidSystem, ::Nothing, semi, v, u) = system + +@inline function merge_particles!(system::FluidSystem, particle_refinement, semi, v, u) + (; delete_candidates) = system.particle_refinement + + delete_candidates .= false + + # Merge particles iteratively + for _ in 1:3 + merge_particles_inner!(system, particle_refinement, semi, v, u) + end + + return system +end -@inline function merge_particles!(system::FluidSystem, particle_refinement, v, u) +function merge_particles_inner!(system, particle_refinement, semi, v, u) (; smoothing_kernel, cache) = system (; mass_ref, max_spacing_ratio, merge_candidates, delete_candidates) = particle_refinement set_zero!(merge_candidates) - system_coords = current_coordinates(u, particle_system) + system_coords = current_coordinates(u, system) neighborhood_search = get_neighborhood_search(system, semi) # Collect merge candidates @@ -64,7 +72,7 @@ end # Merge and delete particles for particle in merge_candidates - iszero(merge_candidates[particle]) && continue + iszero(particle) && continue candidate = merge_candidates[particle] From b5abc331b2599deda4902a51fb4a8009167fce17 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:52:02 +0100 Subject: [PATCH 29/64] fix iterator of merge candidates --- src/refinement/merge.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/refinement/merge.jl b/src/refinement/merge.jl index 99b44c562..ea55bb0bc 100644 --- a/src/refinement/merge.jl +++ b/src/refinement/merge.jl @@ -71,8 +71,7 @@ function merge_particles_inner!(system, particle_refinement, semi, v, u) end # Merge and delete particles - for particle in merge_candidates - iszero(particle) && continue + for particle in findall(!iszero, merge_candidates) candidate = merge_candidates[particle] @@ -104,13 +103,14 @@ function merge_particles_inner!(system, particle_refinement, semi, v, u) # TODO: # Check normalization_factor = smoothing_kernel(zero(pos_diff), one(h_a)) tmp_m = m_merge * normalization_factor(smoothing_kernel, h_a) - tmp_a = m_a * smoothing_kernel(pos_merge - pos_a, h_a) - tmp_b = m_b * smoothing_kernel(pos_merge - pos_b, h_b) + + tmp_a = m_a * kernel(smoothing_kernel, norm(pos_merge - pos_a), h_a) + tmp_b = m_b * kernel(smoothing_kernel, norm(pos_merge - pos_b), h_b) cache.smoothing_length[particle] = (tmp_m / (tmp_a + tmp_b))^(1 / ndims(system)) - system.mass[particl] = m_merge + system.mass[particle] = m_merge else delete_candidates[candidate] = true end From b9d763b3f806426f3160f29086d8bd1e4ce6af4c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:23:30 +0100 Subject: [PATCH 30/64] rename to refinement_pattern --- src/refinement/refinement.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index dc01d85e2..0197e12d5 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -4,7 +4,7 @@ include("split.jl") include("merge.jl") struct ParticleRefinement{SP, RC, ELTYPE} - splitting_pattern :: SP + refinement_pattern :: SP refinement_criteria :: RC max_spacing_ratio :: ELTYPE mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles @@ -14,7 +14,7 @@ struct ParticleRefinement{SP, RC, ELTYPE} n_new_particles :: Int end -function ParticleRefinement(; splitting_pattern, max_spacing_ratio, +function ParticleRefinement(; refinement_pattern, max_spacing_ratio, refinement_criteria=SpatialRefinementCriterion()) mass_ref = Vector{eltype(max_spacing_ratio)}() delete_candidates = Vector{Bool}() @@ -23,7 +23,7 @@ function ParticleRefinement(; splitting_pattern, max_spacing_ratio, refinement_criteria = (refinement_criteria,) end - return ParticleRefinement(splitting_pattern, refinement_criteria, max_spacing_ratio, + return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, mass_ref, Int[], delete_candidates, 0, 0) end From 812ecbc0a4782123228016f78f3616273ff10eea Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:33:38 +0100 Subject: [PATCH 31/64] fix `Ref` --- src/refinement/refinement.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 0197e12d5..b7d086d97 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -4,14 +4,14 @@ include("split.jl") include("merge.jl") struct ParticleRefinement{SP, RC, ELTYPE} - refinement_pattern :: SP + refinement_pattern :: SP refinement_criteria :: RC max_spacing_ratio :: ELTYPE mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles - n_particles_before_resize :: Int - n_new_particles :: Int + n_particles_before_resize :: Ref{Int} + n_new_particles :: Ref{Int} end function ParticleRefinement(; refinement_pattern, max_spacing_ratio, @@ -24,7 +24,7 @@ function ParticleRefinement(; refinement_pattern, max_spacing_ratio, end return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, Int[], delete_candidates, 0, 0) + mass_ref, Int[], delete_candidates, Ref(0), Ref(0)) end resize_refinement!(system) = system From c64e870bb5e6a83317c2fd0667c4b89a4b5aca8b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 14:47:45 +0100 Subject: [PATCH 32/64] add vector for split candidates --- src/refinement/refinement.jl | 3 +- src/refinement/split.jl | 56 +++++++++++++++++------------------- 2 files changed, 28 insertions(+), 31 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index b7d086d97..bf4d375ee 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -10,6 +10,7 @@ struct ParticleRefinement{SP, RC, ELTYPE} mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles + split_candidates :: Vector{Int} n_particles_before_resize :: Ref{Int} n_new_particles :: Ref{Int} end @@ -24,7 +25,7 @@ function ParticleRefinement(; refinement_pattern, max_spacing_ratio, end return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, Int[], delete_candidates, Ref(0), Ref(0)) + mass_ref, Int[], delete_candidates, Int[], Ref(0), Ref(0)) end resize_refinement!(system) = system diff --git a/src/refinement/split.jl b/src/refinement/split.jl index a04000b4e..39bef2a98 100644 --- a/src/refinement/split.jl +++ b/src/refinement/split.jl @@ -25,22 +25,23 @@ end @inline collect_split_candidates!(system::FluidSystem, ::Nothing, v, u) = system @inline function collect_split_candidates!(system::FluidSystem, particle_refinement, v, u) - (; mass_ref, max_spacing_ratio, refinement_pattern) = particle_refinement + (; mass_ref, max_spacing_ratio, split_candidates, refinement_pattern) = particle_refinement - n_split_candidates = 0 + resize!(split_candidates, 0) for particle in eachparticle(system) m_a = hydrodynamic_mass(system, particle) m_max = max_spacing_ratio * mass_ref[particle] if m_a > m_max - n_split_candidates += 1 + push!(split_candidates, particle) end end n_childs_exclude_center = nchilds(system, refinement_pattern) - 1 - particle_refinement.n_new_particles[] = n_split_candidates * n_childs_exclude_center + particle_refinement.n_new_particles[] = length(split_candidates) * + n_childs_exclude_center return system end @@ -55,46 +56,41 @@ end @inline function split_particles!(system::FluidSystem, particle_refinement, v, u) (; smoothing_length) = system.cache - (; mass_ref, max_spacing_ratio, refinement_pattern, n_particles_before_resize) = particle_refinement + (; split_candidates, refinement_pattern, n_particles_before_resize) = particle_refinement (; alpha, relative_position) = refinement_pattern - for particle in eachparticle(system) - m_a = hydrodynamic_mass(system, particle) - m_max = max_spacing_ratio * mass_ref[particle] - - if m_a > m_max - smoothing_length_old = smoothing_length[particle] - mass_old = system.mass[particle] + for particle in split_candidates + smoothing_length_old = smoothing_length[particle] + mass_old = system.mass[particle] - system.mass[particle] = mass_old / nchilds(system, refinement_pattern) + system.mass[particle] = mass_old / nchilds(system, refinement_pattern) - set_particle_pressure!(v, system, partice, pressure) + set_particle_pressure!(v, system, partice, pressure) - set_particle_density!(v, system, partice, density) + set_particle_density!(v, system, partice, density) - smoothing_length[particle] = alpha * smoothing_length_old + smoothing_length[particle] = alpha * smoothing_length_old - pos_center = current_coords(system, u, particle) - vel_center = current_velocity(system, v, particle) + pos_center = current_coords(system, u, particle) + vel_center = current_velocity(system, v, particle) - for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) - child = n_particles_before_resize + child_id_local + for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) + child = n_particles_before_resize + child_id_local - system.mass[child] = mass_old / nchilds(system, refinement_pattern) + system.mass[child] = mass_old / nchilds(system, refinement_pattern) - set_particle_pressure!(v, system, child, pressure) + set_particle_pressure!(v, system, child, pressure) - set_particle_density!(v, system, child, density) + set_particle_density!(v, system, child, density) - smoothing_length[child] = alpha * smoothing_length_old + smoothing_length[child] = alpha * smoothing_length_old - rel_pos = smoothing_length_old * relative_position[child] - new_pos = pos_center + rel_pos + rel_pos = smoothing_length_old * relative_position[child] + new_pos = pos_center + rel_pos - for dim in 1:ndims(system) - u[dim, child] = new_pos[dim] - v[dim, child] = vel_center[dim] - end + for dim in 1:ndims(system) + u[dim, child] = new_pos[dim] + v[dim, child] = vel_center[dim] end end end From 3c53c8b4abd4df2c4b75cb426e4af9ab8cafe4c5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:23:36 +0100 Subject: [PATCH 33/64] fix typo --- src/refinement/split.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/refinement/split.jl b/src/refinement/split.jl index 39bef2a98..88851ae90 100644 --- a/src/refinement/split.jl +++ b/src/refinement/split.jl @@ -65,9 +65,9 @@ end system.mass[particle] = mass_old / nchilds(system, refinement_pattern) - set_particle_pressure!(v, system, partice, pressure) + set_particle_pressure!(v, system, particle, pressure) - set_particle_density!(v, system, partice, density) + set_particle_density!(v, system, particle, density) smoothing_length[particle] = alpha * smoothing_length_old From cfd06eb2429a0b8dca7c24a14934f1fc8fbeadc7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:37:10 +0100 Subject: [PATCH 34/64] fix ids --- src/refinement/split.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/refinement/split.jl b/src/refinement/split.jl index 88851ae90..47c83fef1 100644 --- a/src/refinement/split.jl +++ b/src/refinement/split.jl @@ -65,27 +65,26 @@ end system.mass[particle] = mass_old / nchilds(system, refinement_pattern) - set_particle_pressure!(v, system, particle, pressure) - - set_particle_density!(v, system, particle, density) + p_a = particle_pressure(v, system, particle) + rho_a = particle_density(v, system, particle) smoothing_length[particle] = alpha * smoothing_length_old - pos_center = current_coords(system, u, particle) - vel_center = current_velocity(system, v, particle) + pos_center = current_coords(u, system, particle) + vel_center = current_velocity(v, system, particle) for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) - child = n_particles_before_resize + child_id_local + child = n_particles_before_resize[] + child_id_local system.mass[child] = mass_old / nchilds(system, refinement_pattern) - set_particle_pressure!(v, system, child, pressure) + set_particle_pressure!(v, system, child, p_a) - set_particle_density!(v, system, child, density) + set_particle_density!(v, system, child, rho_a) smoothing_length[child] = alpha * smoothing_length_old - rel_pos = smoothing_length_old * relative_position[child] + rel_pos = smoothing_length_old * relative_position[child_id_local] new_pos = pos_center + rel_pos for dim in 1:ndims(system) From dd2dfdce7f5ec0a2e4c45326296a486001d4cf32 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 14:47:45 +0100 Subject: [PATCH 35/64] add vector for split candidates --- src/refinement/refinement.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index bf4d375ee..c611c68ed 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -10,7 +10,7 @@ struct ParticleRefinement{SP, RC, ELTYPE} mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles - split_candidates :: Vector{Int} + split_candidates :: Vector{Int} # variable length n_particles_before_resize :: Ref{Int} n_new_particles :: Ref{Int} end From c39505e4452a95dcd78a839b8ac6f612690a4697 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 14:39:22 +0100 Subject: [PATCH 36/64] fix global capacity --- src/refinement/resize.jl | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index a7e70fa47..eede0fc00 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -1,7 +1,7 @@ function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) # Resize all systems foreach_system(semi) do system - resize!(system, capacity(system)) + capacity(system) > nparticles(system) && resize!(system, capacity(system)) end resize!(v_ode, u_ode, _v_ode, _u_ode, semi) @@ -23,17 +23,19 @@ function Base.deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) end function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) + (; systems) = semi + copyto!(_v_ode, v_ode) copyto!(_u_ode, u_ode) # Get ranges after resizing the systems - ranges_v_new, ranges_u_new = ranges_vu(semi.systems) + ranges_v_new, ranges_u_new = ranges_vu(systems) ranges_v_old = copy(semi.ranges_v) ranges_u_old = copy(semi.ranges_u) # Set ranges after resizing the systems - for i in 1:length(semi.systems) + for i in 1:length(systems) semi.ranges_v[i] = ranges_v_new[i] semi.ranges_u[i] = ranges_u_new[i] end @@ -50,13 +52,14 @@ function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) end end - capacity_global = sum(system -> nparticles(system), semi.systems) + sizes_u = sum(u_nvariables(system) * n_moving_particles(system) for system in systems) + sizes_v = sum(v_nvariables(system) * n_moving_particles(system) for system in systems) - resize!(v_ode, capacity_global) - resize!(u_ode, capacity_global) + resize!(v_ode, sizes_v) + resize!(u_ode, sizes_u) - resize!(_v_ode, capacity_global) - resize!(_u_ode, capacity_global) + resize!(_v_ode, sizes_v) + resize!(_u_ode, sizes_u) # TODO: Do the following in the callback # resize!(integrator, (length(v_ode), length(u_ode))) From 3967e9257b4374ef2ba1f3c54fd7eaf01c5a82b1 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 37/64] Make it work with EDAC --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 75c0066c3..f4dc04806 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,6 +147,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 52bc334e8ef3a9f55ef24e53f5099bbea15a3191 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 38/64] add h-factor --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f4dc04806..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,18 +147,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 375989c0af26e13f5a823c671d4f5fe54449565e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:06:13 +0100 Subject: [PATCH 39/64] add `ParticleRefinement` --- examples/fluid/particle_refinement_2d.jl | 76 ++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 examples/fluid/particle_refinement_2d.jl diff --git a/examples/fluid/particle_refinement_2d.jl b/examples/fluid/particle_refinement_2d.jl new file mode 100644 index 000000000..0675e63a8 --- /dev/null +++ b/examples/fluid/particle_refinement_2d.jl @@ -0,0 +1,76 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (1.0, 1.0) +tank_size = (1.0, 1.0) + +fluid_density = 1000.0 +sound_speed = 10.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, acceleration=(0.0, -gravity), + state_equation=state_equation) + +# ========================================================================================== +# ==== Fluid + +smoothing_length = 1.2 * fluid_particle_spacing +smoothing_kernel = SchoenbergQuarticSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() +pressure = 1000.0 +particle_refinement = ParticleRefinement(; + refinement_pattern=TrixiParticles.HexagonalSplitting(), + max_spacing_ratio=1.2) +fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, smoothing_length, + sound_speed; viscosity=ViscosityAdami(; nu=1e-4), + particle_refinement=particle_refinement, + transport_velocity=TransportVelocityAdami(pressure), + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary + +# This is to set another boundary density calculation with `trixi_include` +boundary_density_calculator = AdamiPressureExtrapolation() + +# This is to set wall viscosity with `trixi_include` +viscosity_wall = nothing +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=nothing) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +# This is to easily add a new callback with `trixi_include` +extra_callback = nothing + +refinement_callback = ParticleRefinementCallback(interval=1) + +callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback(), refinement_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); From eb838bc7179cf376fc0b7ff5376a920c59403d7b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 16:02:05 +0100 Subject: [PATCH 40/64] export `ParticleRefinementCallback` --- src/TrixiParticles.jl | 2 +- src/callbacks/refinement.jl | 7 ++----- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 54bfb7957..b84d9fdfb 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -62,7 +62,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow, OutFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback, UpdateCallback + PostprocessCallback, StepsizeCallback, UpdateCallback, ParticleRefinementCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/refinement.jl b/src/callbacks/refinement.jl index 130de8756..8f04bed40 100644 --- a/src/callbacks/refinement.jl +++ b/src/callbacks/refinement.jl @@ -21,12 +21,12 @@ function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) if dt > 0 # Add a `tstop` every `dt`, and save the final solution. return PeriodicCallback(refinement_callback, dt, - initialize=initial_update!, + initialize=initial_refinement!, save_positions=(false, false)) else # The first one is the `condition`, the second the `affect!` return DiscreteCallback(refinement_callback, refinement_callback, - initialize=initial_update!, + initialize=initial_refinement!, save_positions=(false, false)) end end @@ -63,9 +63,6 @@ function (refinement_callback::ParticleRefinementCallback)(integrator) 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 From 04a5dc7b7135455edc747fa01e4f7053813aaff0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 15:18:47 +0100 Subject: [PATCH 41/64] add beta correction --- src/schemes/fluid/fluid.jl | 53 ++++++++++++++++++++++++- src/schemes/fluid/transport_velocity.jl | 6 --- 2 files changed, 52 insertions(+), 7 deletions(-) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index f287cabba..ab880a2c4 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -19,8 +19,11 @@ end function create_cache_refinement(initial_condition, refinement, smoothing_length) smoothng_length_factor = smoothing_length / initial_condition.particle_spacing + + beta = Vector{eltype(initial_condition)}(undef, nparticles(initial_condition)) + return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)), - smoothing_length_factor=smoothng_length_factor) + smoothing_length_factor=smoothng_length_factor, beta=beta) end @propagate_inbounds hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] @@ -128,3 +131,51 @@ end v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) end + +function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + # Check if `UpdateCallback` is used when simulating with TVF + update_final!(system, system.transport_velocity, + v, u, v_ode, u_ode, semi, t; update_from_callback) + + # Compute correction factor when using particle refinement + compute_beta_correction!(system, system.particle_refinement, v_ode, u_ode, semi) + + return system +end + +compute_beta_correction!(system, ::Nothing, v_ode, u_ode, semi) = system + +function compute_beta_correction!(system, refinement, v_ode, u_ode, semi) + (; beta) = system.cache + + set_zero!(beta) + + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, semi) + + system_coords = current_coordinates(u, system) + + foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, + pos_diff, distance + rho_a = particle_density(v, system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + W_deriv = kernel_deriv(smoothing_kernel, distance, + smoothing_length(system, particle)) + + beta[particle] -= m_b * distance * W_deriv * (rho_a * ndims(system)) + end + end + + return system +end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 5270d1210..acd6050dd 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -177,12 +177,6 @@ function update_callback_used!(system, transport_velocity) system.cache.update_callback_used[] = true end -function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) - update_final!(system, system.transport_velocity, - v, u, v_ode, u_ode, semi, t; update_from_callback) -end - function update_final!(system::FluidSystem, ::Nothing, v, u, v_ode, u_ode, semi, t; update_from_callback=false) return system From 9b218aac81d989a267d8b208e5769ecf35e27f73 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 15:35:43 +0100 Subject: [PATCH 42/64] fix kinematic viscosity --- src/schemes/fluid/fluid.jl | 2 +- src/schemes/fluid/viscosity.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index ab880a2c4..4a87c8646 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -170,7 +170,7 @@ function compute_beta_correction!(system, refinement, v_ode, u_ode, semi) rho_a = particle_density(v, system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - W_deriv = kernel_deriv(smoothing_kernel, distance, + W_deriv = kernel_deriv(system.smoothing_kernel, distance, smoothing_length(system, particle)) beta[particle] -= m_b * distance * W_deriv * (rho_a * ndims(system)) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 222fa6fd8..c43cd32d8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -231,7 +231,8 @@ end beta_inv_a = beta_correction(particle_system, particle_refinement, particle) nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system)) + viscosity_model(neighbor_system, particle_system), + smoothing_length(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) @@ -291,7 +292,7 @@ end return visc .* v_diff end -function kinematic_viscosity(system, viscosity::ViscosityAdami, particle) +function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length) return viscosity.nu end From 46ad110e9d752d952a19346a2f1688950db608f8 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 15:40:12 +0100 Subject: [PATCH 43/64] viscosty fix --- src/schemes/fluid/viscosity.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index c43cd32d8..390dffb28 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -234,6 +234,10 @@ end viscosity_model(neighbor_system, particle_system), smoothing_length(particle_system, particle)) + 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 + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) From 0fd34f1d9bf42a58de8553da73867aab4d5ddb98 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 15:53:33 +0100 Subject: [PATCH 44/64] final fix --- .../boundary/dummy_particles/dummy_particles.jl | 4 ++++ src/schemes/boundary/system.jl | 4 ++++ .../fluid/entropically_damped_sph/rhs.jl | 17 ++++++++++++++--- src/schemes/fluid/transport_velocity.jl | 2 +- 4 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index ece912eeb..dcf3d291a 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(system, boundary_model::BoundaryModelDummyParticles, particle) + return boundary_model.smoothing_length +end + function smoothing_length(boundary_model::BoundaryModelDummyParticles, particle) return boundary_model.smoothing_length end diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index a3f47f611..d37fed82e 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -328,6 +328,10 @@ end return kernel(smoothing_kernel, distance, smoothing_length) end +@inline function smoothing_length(system::BoundarySPHSystem, particle) + return smoothing_length(system, system.boundary_model, particle) +end + function update_positions!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t) (; movement) = system diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index b075813f2..7e81ded7a 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -152,6 +152,8 @@ end 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) + (; sound_speed) = particle_system + # EDAC pressure evolution pressure_diff = p_a - p_b @@ -162,7 +164,7 @@ function pressure_damping_term(particle_system, neighbor_system, refinement, 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_b = alpha_edac * sound_speed * smoothing_length(neighbor_system, neighbor) / tmp nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) @@ -207,6 +209,15 @@ end 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::EntropicallyDampedSPHSystem, + 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, @@ -226,8 +237,8 @@ end 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/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index acd6050dd..374c81146 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -112,7 +112,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)' From 65c3cf77768e8e1c4c0067b82edf1619bd61cca9 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Nov 2024 17:15:23 +0100 Subject: [PATCH 45/64] adapt for tests --- src/refinement/particle_refinement.jl | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 src/refinement/particle_refinement.jl diff --git a/src/refinement/particle_refinement.jl b/src/refinement/particle_refinement.jl new file mode 100644 index 000000000..b5d6f8340 --- /dev/null +++ b/src/refinement/particle_refinement.jl @@ -0,0 +1,8 @@ +struct ParticleRefinement + n_particles_before_resize :: Ref{Int} + n_new_particles :: Ref{Int} +end + +function ParticleRefinement() + return ParticleRefinement(Ref(0), Ref(0)) +end From 6363c7574e1192c2c18096a12dcb30c20fa8ba43 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:24:47 +0100 Subject: [PATCH 46/64] restructure --- src/refinement/particle_refinement.jl | 8 -------- 1 file changed, 8 deletions(-) delete mode 100644 src/refinement/particle_refinement.jl diff --git a/src/refinement/particle_refinement.jl b/src/refinement/particle_refinement.jl deleted file mode 100644 index b5d6f8340..000000000 --- a/src/refinement/particle_refinement.jl +++ /dev/null @@ -1,8 +0,0 @@ -struct ParticleRefinement - n_particles_before_resize :: Ref{Int} - n_new_particles :: Ref{Int} -end - -function ParticleRefinement() - return ParticleRefinement(Ref(0), Ref(0)) -end From 4e8691302bbc58a68fa317ee65d9609ca8541929 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 12:38:05 +0100 Subject: [PATCH 47/64] First prototype of variable smoothing length --- .../fluid/weakly_compressible_sph/system.jl | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 3ea8352f7..24676c8b5 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -160,6 +160,26 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) return (; surface_normal) end +function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles) + return (; smoothing_length) +end + +function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles) + return (; smoothing_length=smoothing_length * ones(n_particles)) +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time From 3a5f60573243c5e80e8fde9eec1e7cfe4cfc1fb3 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 48/64] Make it work with EDAC --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 75c0066c3..f4dc04806 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,6 +147,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 18f45e6bccc0e547beaebf836501aebfe0de36e8 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 49/64] add h-factor --- .../fluid/entropically_damped_sph/system.jl | 12 ----------- .../fluid/weakly_compressible_sph/system.jl | 20 ------------------- 2 files changed, 32 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f4dc04806..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,18 +147,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 24676c8b5..3ea8352f7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -160,26 +160,6 @@ function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) return (; surface_normal) end -function create_cache_refinement(::Nothing, smoothing_length, ELTYPE, NDIMS, n_particles) - return (; smoothing_length) -end - -function create_cache_refinement(refinement, smoothing_length, ELTYPE, NDIMS, n_particles) - return (; smoothing_length=smoothing_length * ones(n_particles)) -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::WeaklyCompressibleSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time From 1ebcef5919abc94b450270c77fba67c1371c0221 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 13:07:20 +0100 Subject: [PATCH 50/64] =?UTF-8?q?ada=C3=BCt=20viscosity?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/fluid/refinement_2d.jl | 72 +++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 examples/fluid/refinement_2d.jl diff --git a/examples/fluid/refinement_2d.jl b/examples/fluid/refinement_2d.jl new file mode 100644 index 000000000..9dc77f906 --- /dev/null +++ b/examples/fluid/refinement_2d.jl @@ -0,0 +1,72 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (1.0, 1.0) +tank_size = (1.0, 1.0) + +fluid_density = 1000.0 +sound_speed = 10.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, acceleration=(0.0, -gravity), + state_equation=state_equation) + +# ========================================================================================== +# ==== Fluid + +smoothing_length = 1.2 * fluid_particle_spacing +smoothing_kernel = SchoenbergQuarticSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() +pressure = 1000.0 +particle_refinement = ParticleRefinement(; splitting_pattern=TrixiParticles.HexagonalSplitting(), max_spacing_ratio=1.2) +fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, smoothing_length, + sound_speed; viscosity=ViscosityAdami(; nu=1e-4), + particle_refinement=particle_refinement, + transport_velocity=TransportVelocityAdami(pressure), + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary + +# This is to set another boundary density calculation with `trixi_include` +boundary_density_calculator = AdamiPressureExtrapolation() + +# This is to set wall viscosity with `trixi_include` +viscosity_wall = nothing +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=nothing) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +# This is to easily add a new callback with `trixi_include` +extra_callback = nothing + +callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback()) + +# Use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); From 93984919eca8a6ea7b4b3eae020b8f6a836d0a88 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Nov 2024 18:48:07 +0100 Subject: [PATCH 51/64] update h --- src/refinement/refinement.jl | 48 ++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index c611c68ed..1ce7f96c1 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -61,6 +61,7 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) # Correct the particles # Update smoothing lengths + upate_smoothing_lengths!(semi, u_ode) # Resize neighborhood search @@ -110,6 +111,53 @@ end return system end +function upate_smoothing_lengths!(semi, u_ode) + foreach_system(semi) do system + u = wrap_u(u_ode, system, semi) + upate_smoothing_lengths!(system, semi, u) + end + + return semi +end + +upate_smoothing_lengths!(system, semi, u) = system + +function upate_smoothing_lengths!(system::FluidSystem, semi, u) + upate_smoothing_lengths!(system, system.particle_refinement, semi, u) +end + +upate_smoothing_lengths!(system, ::Nothing, semi, u) = system + +function upate_smoothing_lengths!(system::FluidSystem, refinement, semi, u) + (; smoothing_length, smoothing_length_factor) = system.cache + + neighborhood_search = get_neighborhood_search(system, semi) + + system_coords = current_coordinates(u, system) + + # TODO: Think about a solution when density is not constant in `initial_condition` + # Problem: Resized systems + rho_0 = first(system.initial_condition.density) + + for particle in each_moving_particle(system) + mass_avg = zero(eltype(system)) + counter_neighbors = 0 + + PointNeighbors.foreach_neighbor(system_coords, system_coords, neighborhood_search, + particle) do particle, neighbor, pos_diff, distance + mass_avg += system.mass[neighbor] + + counter_neighbors += 1 + end + + mass_avg /= counter_neighbors + + smoothing_length[particle] = smoothing_length_factor * + (mass_avg / rho_0)^(1 / ndims(system)) + end + + return system +end @inline function min_max_avg_spacing(system, semi, u_ode, system_coords, particle) dp_min = Inf From b2c4d1268f5274736ee18aaa0851f7b951742f2c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 52/64] Make it work with EDAC --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 75c0066c3..f4dc04806 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,6 +147,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 0efc4e1f4cf34cb9fa91cd3671fb66fbe155a243 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 53/64] add h-factor --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f4dc04806..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -147,18 +147,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 57a9f130a95f1aa76a84372d83dc85dda348a31b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Nov 2024 17:15:23 +0100 Subject: [PATCH 54/64] adapt for tests --- src/schemes/fluid/entropically_damped_sph/system.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 75c0066c3..b769014be 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -100,15 +100,12 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, 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...) 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) + smoothing_kernel, smoothing_length, 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) From 87f6eabc9f9dea6d44a31ffe2b4d8ded0da1d6c7 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 6 Nov 2024 15:57:35 +0100 Subject: [PATCH 55/64] Make it work with EDAC --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index b769014be..399be8cb2 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -144,6 +144,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +function smoothing_length(system::EntropicallyDampedSPHSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end + create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From 0b05aae7a9ef5977c1e9ed8a4e187c60c7deecc7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Nov 2024 10:53:40 +0100 Subject: [PATCH 56/64] add h-factor --- src/schemes/fluid/entropically_damped_sph/system.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 399be8cb2..b769014be 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -144,18 +144,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -function smoothing_length(system::EntropicallyDampedSPHSystem, particle) - return smoothing_length(system, system.particle_refinement, particle) -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, ::Nothing, particle) - return system.cache.smoothing_length -end - -function smoothing_length(system::EntropicallyDampedSPHSystem, refinement, particle) - return system.cache.smoothing_length[particle] -end - create_cache_edac(initial_condition, ::Nothing) = (;) function create_cache_edac(initial_condition, ::TransportVelocityAdami) From d5abc14b9d844c676c35e92d04de2acf71e0b490 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 18 Nov 2024 14:41:30 +0100 Subject: [PATCH 57/64] Make simulations run with `Float32` (#662) * Replace hardcoded Float64 values * Avoid Float64 `eps()` * Remove the last Float64 occurrence in fluid-fluid kernel * Fix tests * Make all smoothing kernels preserve types and add test for this * Make eps(typeof(h)) consistent * Reformat code * Fix tests --- src/schemes/fluid/viscosity.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 390dffb28..7cef990a5 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -112,9 +112,12 @@ end v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, - rho_a, rho_b, grad_kernel) - rho_mean = 0.5 * (rho_a + rho_b) + sound_speed, + m_a, m_b, rho_a, rho_b, + grad_kernel) + (; smoothing_length) = particle_system + + rho_mean = (rho_a + rho_b) / 2 v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) From 6235fc2453cf7bddd39b96dfd2373700d8e373a6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Nov 2024 11:11:04 +0100 Subject: [PATCH 58/64] add cache --- examples/fluid/refinement_2d.jl | 72 ------------------- .../fluid/entropically_damped_sph/system.jl | 11 +-- 2 files changed, 7 insertions(+), 76 deletions(-) delete mode 100644 examples/fluid/refinement_2d.jl diff --git a/examples/fluid/refinement_2d.jl b/examples/fluid/refinement_2d.jl deleted file mode 100644 index 9dc77f906..000000000 --- a/examples/fluid/refinement_2d.jl +++ /dev/null @@ -1,72 +0,0 @@ -using TrixiParticles -using OrdinaryDiffEq - -# ========================================================================================== -# ==== Resolution -fluid_particle_spacing = 0.05 - -# Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 3 - -# ========================================================================================== -# ==== Experiment Setup -gravity = 9.81 -tspan = (0.0, 1.0) - -# Boundary geometry and initial fluid particle positions -initial_fluid_size = (1.0, 1.0) -tank_size = (1.0, 1.0) - -fluid_density = 1000.0 -sound_speed = 10.0 -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7, clip_negative_pressure=false) - -tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, - n_layers=boundary_layers, acceleration=(0.0, -gravity), - state_equation=state_equation) - -# ========================================================================================== -# ==== Fluid - -smoothing_length = 1.2 * fluid_particle_spacing -smoothing_kernel = SchoenbergQuarticSplineKernel{2}() - -fluid_density_calculator = ContinuityDensity() -pressure = 1000.0 -particle_refinement = ParticleRefinement(; splitting_pattern=TrixiParticles.HexagonalSplitting(), max_spacing_ratio=1.2) -fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, smoothing_length, - sound_speed; viscosity=ViscosityAdami(; nu=1e-4), - particle_refinement=particle_refinement, - transport_velocity=TransportVelocityAdami(pressure), - acceleration=(0.0, -gravity)) - -# ========================================================================================== -# ==== Boundary - -# This is to set another boundary density calculation with `trixi_include` -boundary_density_calculator = AdamiPressureExtrapolation() - -# This is to set wall viscosity with `trixi_include` -viscosity_wall = nothing -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - boundary_density_calculator, - smoothing_kernel, smoothing_length, - viscosity=viscosity_wall) -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=nothing) - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system) -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=50) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") - -# This is to easily add a new callback with `trixi_include` -extra_callback = nothing - -callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback()) - -# Use a Runge-Kutta method with automatic (error based) time step size control -sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index b769014be..75c0066c3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -100,12 +100,15 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, 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...) return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, - smoothing_kernel, smoothing_length, sound_speed, - viscosity, nu_edac, acceleration_, nothing, - pressure_acceleration, transport_velocity, - source_terms, buffer, particle_refinement, cache) + 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) From 0b5202dd6c0962315aac4188312e8a57961a3e0d Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:38:47 +0100 Subject: [PATCH 59/64] resize cache --- src/refinement/resize.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index eede0fc00..6e5f48971 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -78,26 +78,26 @@ end Base.resize!(system, ::Nothing, capacity_system) = system function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) - (; mass, pressure, cache, density_calculator) = system + (; mass, pressure, density_calculator) = system refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) resize!(pressure, capacity_system) resize_density!(system, capacity_system, density_calculator) - # TODO - # resize_cache!(system, cache, n) + + resize_cache!(system, capacity_system) end function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) - (; mass, cache, density_calculator) = system + (; mass, density_calculator) = system refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) resize_density!(system, capacity_system, density_calculator) - # TODO - # resize_cache!(system, capacity_system) + + resize_cache!(system, capacity_system) return system end @@ -105,14 +105,15 @@ end resize_density!(system, n::Int, ::SummationDensity) = resize!(system.cache.density, n) resize_density!(system, n::Int, ::ContinuityDensity) = system -function resize_cache!(system, n::Int) +function resize_cache!(system::WeaklyCompressibleSPHSystem, n::Int) resize!(system.cache.smoothing_length, n) return system end -function resize_cache!(system::EntropicallyDampedSPHSystem, n) +function resize_cache!(system::EntropicallyDampedSPHSystem, n::Int) resize!(system.cache.smoothing_length, n) + resize!(system.cache.beta, n) resize!(system.cache.pressure_average, n) resize!(system.cache.neighbor_counter, n) From f4f5d389caefc28973c5d0622229189b424a9b19 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 15:43:41 +0100 Subject: [PATCH 60/64] resize refinement --- src/refinement/resize.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index 6e5f48971..4bfc60e28 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -84,9 +84,14 @@ function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_ resize!(mass, capacity_system) resize!(pressure, capacity_system) + resize_density!(system, capacity_system, density_calculator) resize_cache!(system, capacity_system) + + resize_refinement!(system) + + return system end function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) @@ -95,10 +100,13 @@ function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_ refinement.n_particles_before_resize[] = nparticles(system) resize!(mass, capacity_system) + resize_density!(system, capacity_system, density_calculator) resize_cache!(system, capacity_system) + resize_refinement!(system) + return system end From 9446c4129ff23b3550104e35e622455d05bcf09f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Nov 2024 16:15:40 +0100 Subject: [PATCH 61/64] fix ambiguous method --- src/refinement/refinement.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 1ce7f96c1..616d3e98d 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -126,7 +126,7 @@ function upate_smoothing_lengths!(system::FluidSystem, semi, u) upate_smoothing_lengths!(system, system.particle_refinement, semi, u) end -upate_smoothing_lengths!(system, ::Nothing, semi, u) = system +upate_smoothing_lengths!(system::FluidSystem, ::Nothing, semi, u) = system function upate_smoothing_lengths!(system::FluidSystem, refinement, semi, u) (; smoothing_length, smoothing_length_factor) = system.cache From e9da8f5d20a2210dbda84cf0f6bc1627e925dae2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 21 Nov 2024 15:32:37 +0100 Subject: [PATCH 62/64] fix resize --- src/refinement/resize.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index 4bfc60e28..a63dbabd8 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -80,8 +80,6 @@ Base.resize!(system, ::Nothing, capacity_system) = system function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) (; mass, pressure, density_calculator) = system - refinement.n_particles_before_resize[] = nparticles(system) - resize!(mass, capacity_system) resize!(pressure, capacity_system) @@ -97,8 +95,6 @@ end function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) (; mass, density_calculator) = system - refinement.n_particles_before_resize[] = nparticles(system) - resize!(mass, capacity_system) resize_density!(system, capacity_system, density_calculator) From 36468f650e332a1b4f3ece2e8f3f06f661afb822 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 21 Nov 2024 15:35:38 +0100 Subject: [PATCH 63/64] fix splitting --- src/refinement/refinement.jl | 43 ++++++++++++++++------------ src/refinement/refinement_pattern.jl | 21 ++++---------- src/refinement/split.jl | 28 +++++++++--------- 3 files changed, 44 insertions(+), 48 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 616d3e98d..dd6dd9a64 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -4,15 +4,14 @@ include("split.jl") include("merge.jl") struct ParticleRefinement{SP, RC, ELTYPE} - refinement_pattern :: SP - refinement_criteria :: RC - max_spacing_ratio :: ELTYPE - mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles - merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles - delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles - split_candidates :: Vector{Int} # variable length - n_particles_before_resize :: Ref{Int} - n_new_particles :: Ref{Int} + refinement_pattern :: SP + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles + delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles + split_candidates :: Vector{Int} # variable length + n_new_particles :: Ref{Int} end function ParticleRefinement(; refinement_pattern, max_spacing_ratio, @@ -25,7 +24,7 @@ function ParticleRefinement(; refinement_pattern, max_spacing_ratio, end return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, Int[], delete_candidates, Int[], Ref(0), Ref(0)) + mass_ref, Int[], delete_candidates, Int[], Ref(0)) end resize_refinement!(system) = system @@ -41,29 +40,35 @@ function resize_refinement!(system, particle_refinement) resize!(particle_refinement.delete_candidates, nparticles(system)) resize!(particle_refinement.merge_candidates, nparticles(system)) + particle_refinement.mass_ref .= system.mass + return system end function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + # Check refinement criteria check_refinement_criteria!(semi, v_ode, u_ode) - # Update the spacing of particles (Algorthm 1) + # Update spacing of particles (Algorthm 1) update_particle_spacing(semi, v_ode, u_ode) - # Split the particles (Algorithm 2) + # Split particles (Algorithm 2) split_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) - # Merge the particles (Algorithm 3) + # Merge particles (Algorithm 3) merge_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) - # Shift the particles + # Shift particles + # TODO - # Correct the particles + # Correct particles + # TODO # Update smoothing lengths upate_smoothing_lengths!(semi, u_ode) # Resize neighborhood search + # TODO return semi end @@ -160,9 +165,9 @@ function upate_smoothing_lengths!(system::FluidSystem, refinement, semi, u) 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)) + dp_min = particle_spacing(system, particle) + dp_max = particle_spacing(system, particle) + dp_avg = particle_spacing(system, particle) counter_neighbors = 0 foreach_system(semi) do neighbor_system @@ -183,7 +188,7 @@ end end end - dp_avg / counter_neighbors + dp_avg = counter_neighbors == 0 ? dp_avg : dp_avg / counter_neighbors return dp_min, dp_max, dp_avg end diff --git a/src/refinement/refinement_pattern.jl b/src/refinement/refinement_pattern.jl index 759fe8262..c01022314 100644 --- a/src/refinement/refinement_pattern.jl +++ b/src/refinement/refinement_pattern.jl @@ -3,6 +3,7 @@ struct CubicSplitting{ELTYPE} alpha :: ELTYPE center_particle :: Bool relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int function CubicSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) ELTYPE = typeof(epsilon) @@ -21,19 +22,16 @@ struct CubicSplitting{ELTYPE} push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end - new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 4 + center_particle) end end -@inline function nchilds(system, refinement_pattern::CubicSplitting) - return 4 + refinement_pattern.center_particle -end - struct TriangularSplitting{ELTYPE} epsilon :: ELTYPE alpha :: ELTYPE center_particle :: Bool relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int function TriangularSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) ELTYPE = typeof(epsilon) @@ -49,19 +47,16 @@ struct TriangularSplitting{ELTYPE} push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end - new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 3 + center_particle) end end -@inline function nchilds(system::System{2}, refinement_pattern::TriangularSplitting) - return 3 + refinement_pattern.center_particle -end - struct HexagonalSplitting{ELTYPE} epsilon :: ELTYPE alpha :: ELTYPE center_particle :: Bool relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int function HexagonalSplitting(; epsilon=0.4, alpha=0.9, center_particle=true) ELTYPE = typeof(epsilon) @@ -84,10 +79,6 @@ struct HexagonalSplitting{ELTYPE} push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end - new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 6 + center_particle) end end - -@inline function nchilds(system::System{2}, refinement_pattern::HexagonalSplitting) - return 6 + refinement_pattern.center_particle -end diff --git a/src/refinement/split.jl b/src/refinement/split.jl index 47c83fef1..b7743d7d9 100644 --- a/src/refinement/split.jl +++ b/src/refinement/split.jl @@ -26,6 +26,7 @@ end @inline function collect_split_candidates!(system::FluidSystem, particle_refinement, v, u) (; mass_ref, max_spacing_ratio, split_candidates, refinement_pattern) = particle_refinement + (; n_children, center_particle) = refinement_pattern resize!(split_candidates, 0) @@ -38,10 +39,8 @@ end end end - n_childs_exclude_center = nchilds(system, refinement_pattern) - 1 - particle_refinement.n_new_particles[] = length(split_candidates) * - n_childs_exclude_center + (n_children - center_particle) return system end @@ -56,14 +55,15 @@ end @inline function split_particles!(system::FluidSystem, particle_refinement, v, u) (; smoothing_length) = system.cache - (; split_candidates, refinement_pattern, n_particles_before_resize) = particle_refinement - (; alpha, relative_position) = refinement_pattern + (; split_candidates, refinement_pattern, n_new_particles) = particle_refinement + (; alpha, relative_position, n_children, center_particle) = refinement_pattern + child_id_global = nparticles(system) - n_new_particles[] for particle in split_candidates smoothing_length_old = smoothing_length[particle] mass_old = system.mass[particle] - system.mass[particle] = mass_old / nchilds(system, refinement_pattern) + system.mass[particle] = mass_old / n_children p_a = particle_pressure(v, system, particle) rho_a = particle_density(v, system, particle) @@ -73,23 +73,23 @@ end pos_center = current_coords(u, system, particle) vel_center = current_velocity(v, system, particle) - for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) - child = n_particles_before_resize[] + child_id_local + for child_id_local in 1:(n_children - center_particle) + child_id_global += 1 - system.mass[child] = mass_old / nchilds(system, refinement_pattern) + system.mass[child_id_global] = mass_old / n_children - set_particle_pressure!(v, system, child, p_a) + set_particle_pressure!(v, system, child_id_global, p_a) - set_particle_density!(v, system, child, rho_a) + set_particle_density!(v, system, child_id_global, rho_a) - smoothing_length[child] = alpha * smoothing_length_old + smoothing_length[child_id_global] = alpha * smoothing_length_old rel_pos = smoothing_length_old * relative_position[child_id_local] new_pos = pos_center + rel_pos for dim in 1:ndims(system) - u[dim, child] = new_pos[dim] - v[dim, child] = vel_center[dim] + u[dim, child_id_global] = new_pos[dim] + v[dim, child_id_global] = vel_center[dim] end end end From 7e7477932c5ddb4f0d41a6f8d3ce256667956b0a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 22 Nov 2024 12:04:05 +0100 Subject: [PATCH 64/64] fix minor --- examples/fluid/particle_refinement_2d.jl | 26 ++++++++++++++++-------- src/refinement/refinement.jl | 2 +- src/refinement/refinement_criteria.jl | 2 +- src/refinement/resize.jl | 7 +++---- 4 files changed, 23 insertions(+), 14 deletions(-) diff --git a/examples/fluid/particle_refinement_2d.jl b/examples/fluid/particle_refinement_2d.jl index 0675e63a8..d873d43f0 100644 --- a/examples/fluid/particle_refinement_2d.jl +++ b/examples/fluid/particle_refinement_2d.jl @@ -14,8 +14,8 @@ gravity = 9.81 tspan = (0.0, 1.0) # Boundary geometry and initial fluid particle positions -initial_fluid_size = (1.0, 1.0) -tank_size = (1.0, 1.0) +initial_fluid_size = (2.0, 2.0) +tank_size = (2.0, 2.0) fluid_density = 1000.0 sound_speed = 10.0 @@ -23,21 +23,24 @@ state_equation = StateEquationCole(; sound_speed, reference_density=fluid_densit exponent=7, clip_negative_pressure=false) tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, - n_layers=boundary_layers, acceleration=(0.0, -gravity), - state_equation=state_equation) + n_layers=boundary_layers, spacing_ratio=2.0) +sphere = SphereShape(0.5 * fluid_particle_spacing, 0.3, (1.0, 1.0), + fluid_density, sphere_type=RoundSphere()) # ========================================================================================== # ==== Fluid +fluid = setdiff(tank.fluid, sphere) + smoothing_length = 1.2 * fluid_particle_spacing -smoothing_kernel = SchoenbergQuarticSplineKernel{2}() +smoothing_kernel = SchoenbergQuinticSplineKernel{2}() fluid_density_calculator = ContinuityDensity() pressure = 1000.0 particle_refinement = ParticleRefinement(; refinement_pattern=TrixiParticles.HexagonalSplitting(), max_spacing_ratio=1.2) -fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, smoothing_length, +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed; viscosity=ViscosityAdami(; nu=1e-4), particle_refinement=particle_refinement, transport_velocity=TransportVelocityAdami(pressure), @@ -57,9 +60,15 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar viscosity=viscosity_wall) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=nothing) +boundary_model_sphere = BoundaryModelDummyParticles(sphere.density, sphere.mass, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) +boundary_system_sphere = BoundarySPHSystem(sphere, boundary_model_sphere, movement=nothing) + # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system) +semi = Semidiscretization(fluid_system, boundary_system, boundary_system_sphere) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) @@ -70,7 +79,8 @@ extra_callback = nothing refinement_callback = ParticleRefinementCallback(interval=1) -callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback(), refinement_callback) +callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback(), + refinement_callback) # Use a Runge-Kutta method with automatic (error based) time step size control sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index dd6dd9a64..399a42bab 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -150,7 +150,7 @@ function upate_smoothing_lengths!(system::FluidSystem, refinement, semi, u) PointNeighbors.foreach_neighbor(system_coords, system_coords, neighborhood_search, particle) do particle, neighbor, pos_diff, distance - mass_avg += system.mass[neighbor] + mass_avg += hydrodynamic_mass(system, neighbor) counter_neighbors += 1 end diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl index 627ffb619..b8106cbb8 100644 --- a/src/refinement/refinement_criteria.jl +++ b/src/refinement/refinement_criteria.jl @@ -59,8 +59,8 @@ end # Only consider particles with a distance > 0. distance < sqrt(eps()) && return + dp_particle = particle_spacing(particle_system, particle) 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 diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl index a63dbabd8..224517282 100644 --- a/src/refinement/resize.jl +++ b/src/refinement/resize.jl @@ -138,7 +138,7 @@ function Base.deleteat!(system::FluidSystem, refinement, v, u) delete_counter = 0 for particle in eachparticle(system) - if !iszero(delete_candidates[particle]) + if delete_candidates[particle] # swap particles (keep -> delete) dump_id = nparticles(system) - delete_counter @@ -148,9 +148,8 @@ function Base.deleteat!(system::FluidSystem, refinement, v, u) mass_keep = hydrodynamic_mass(system, dump_id) density_keep = particle_density(v, system, dump_id) pressure_keep = particle_pressure(v, system, dump_id) - #TODO - # smoothing_length_keep = smoothing_length(system, dump_id) - # system.cache.smoothing_length[particle] = smoothing_length_keep + + system.cache.smoothing_length[particle] = smoothing_length(system, dump_id) system.mass[particle] = mass_keep