From 23f5e14086c58561025e4895f28d766cb9fbe8d5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 15:36:58 +0200 Subject: [PATCH 01/62] add average pressure --- .../fluid/entropically_damped_sph/rhs.jl | 5 +- .../fluid/entropically_damped_sph/system.jl | 85 +++++++++++++++++-- 2 files changed, 84 insertions(+), 6 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index cac6f3ed6..019d0e878 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -30,7 +30,8 @@ function interact!(dv, v_particle_system, u_particle_system, volume_term = (volume_a^2 + volume_b^2) / m_a # Inter-particle averaged pressure - pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b) + p_avg = average_pressure(particle_system, particle) # Only used with TVF + pressure_tilde = (rho_b * (p_a - p_avg) + rho_a * (p_b - p_avg)) / (rho_a + rho_b) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) @@ -80,5 +81,7 @@ end return dv end +@inline average_pressure(system, particle) = 0.0 + @inline viscosity_function(system) = system.viscosity @inline viscosity_function(system::BoundarySPHSystem) = system.boundary_model.viscosity diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 8b9d076e5..240681484 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -31,7 +31,7 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). In: Computers and Fluids 179 (2019), pages 579-594. [doi: 10.1016/j.compfluid.2018.11.023](https://doi.org/10.1016/j.compfluid.2018.11.023) """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF} <: +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, TV, C} <: FluidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] @@ -44,11 +44,14 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF} <: nu_edac :: ELTYPE initial_pressure_function :: PF acceleration :: SVector{NDIMS, ELTYPE} + transport_velocity :: TV + cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; alpha=0.5, viscosity=NoViscosity(), initial_pressure_function=nothing, + transport_velocity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel))) NDIMS = ndims(initial_condition) @@ -71,11 +74,15 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF} <: density_calculator = SummationDensity() + cache = create_cache(initial_condition, transport_velocity) + new{NDIMS, ELTYPE, typeof(density_calculator), typeof(smoothing_kernel), - typeof(viscosity), typeof(initial_pressure_function) - }(initial_condition, mass, density, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, initial_pressure_function, - acceleration_) + typeof(viscosity), typeof(initial_pressure_function), + typeof(initial_velocity_function), typeof(transport_velocity), + typeof(cache)}(initial_condition, mass, density, density_calculator, + smoothing_kernel, smoothing_length, sound_speed, viscosity, + nu_edac, initial_pressure_function, initial_velocity_function, + acceleration_, transport_velocity, cache) end end @@ -105,6 +112,16 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end +create_cache(initial_condition, ::Nothing) = (;) + +function create_cache(initial_condition, ::TransportVelocityAdami) + pressure_average = copy(initial_condition.pressure) + neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) + advection_velocity = copy(initial_condition.velocity) + + return (; pressure_average, neighbor_counter, advection_velocity) +end + @inline function particle_density(v, system::EntropicallyDampedSPHSystem, particle) return system.density[particle] end @@ -113,11 +130,69 @@ end return v[end, particle] end +@inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) + average_pressure(system, system.transport_velocity, particle) +end + +@inline function average_pressure(system, ::TransportVelocityAdami, particle) + return system.cache.pressure_average[particle] +end + +@inline average_pressure(system, ::Nothing, particle) = 0.0 + @inline v_nvariables(system::EntropicallyDampedSPHSystem) = ndims(system) + 1 function update_quantities!(system::EntropicallyDampedSPHSystem, system_index, v, u, v_ode, u_ode, semi, t) summation_density!(system, system_index, semi, u, u_ode, system.density) + update_average_pressure!(system, system.transport_velocity, system_index, + v_ode, u_ode, semi) +end + +function update_average_pressure!(system, ::Nothing, system_index, v_ode, u_ode, semi) + return system +end + +function update_average_pressure!(system, ::TransportVelocityAdami, system_index, + v_ode, u_ode, semi) + (; systems, neighborhood_searches) = semi + (; cache) = system + (; pressure_average, neighbor_counter) = cache + + set_zero!(pressure_average) + set_zero!(neighbor_counter) + + u = wrap_u(u_ode, system_index, system, semi) + + # Use all other systems for the average pressure + @trixi_timeit timer() "compute average pressure" foreach_enumerate(systems) do (neighbor_system_index, + neighbor_system) + u_neighbor_system = wrap_u(u_ode, neighbor_system_index, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system_index, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + neighborhood_search = neighborhood_searches[system_index][neighbor_system_index] + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + for_particle_neighbor(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 + + pressure_average[particle] += particle_pressure(v_neighbor_system, + neighbor_system, + neighbor) + neighbor_counter[particle] += 1 + end + end + + for particle in eachparticle(system) + if neighbor_counter[particle] > 0 + pressure_average[particle] /= neighbor_counter[particle] + end + end end function write_v0!(v0, system::EntropicallyDampedSPHSystem) From 5906f80e582aa580c476ed353234c435688b808c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 15:45:15 +0200 Subject: [PATCH 02/62] add transport velocity --- src/TrixiParticles.jl | 2 +- .../fluid/entropically_damped_sph/system.jl | 74 ++++++++++++++++++- src/schemes/fluid/fluid.jl | 1 + src/schemes/fluid/transport_velocity.jl | 30 ++++++++ 4 files changed, 102 insertions(+), 5 deletions(-) create mode 100644 src/schemes/fluid/transport_velocity.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 698835972..c01f9f36d 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -36,7 +36,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem export InfoCallback, SolutionSavingCallback export ContinuityDensity, SummationDensity -export PenaltyForceGanzenmueller +export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel export StateEquationIdealGas, StateEquationCole diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 240681484..bc5c12001 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -31,7 +31,7 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). In: Computers and Fluids 179 (2019), pages 579-594. [doi: 10.1016/j.compfluid.2018.11.023](https://doi.org/10.1016/j.compfluid.2018.11.023) """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, TV, C} <: +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, VF, TV, C} <: FluidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] @@ -43,6 +43,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, TV, C} < viscosity :: V nu_edac :: ELTYPE initial_pressure_function :: PF + initial_velocity_function :: VF acceleration :: SVector{NDIMS, ELTYPE} transport_velocity :: TV cache :: C @@ -51,6 +52,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, TV, C} < smoothing_length, sound_speed; alpha=0.5, viscosity=NoViscosity(), initial_pressure_function=nothing, + initial_velocity_function=nothing, transport_velocity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel))) @@ -140,7 +142,39 @@ end @inline average_pressure(system, ::Nothing, particle) = 0.0 -@inline v_nvariables(system::EntropicallyDampedSPHSystem) = ndims(system) + 1 +@inline function v_nvariables(system::EntropicallyDampedSPHSystem) + v_nvariables(system, system.transport_velocity) +end + +@inline v_nvariables(system, ::Nothing) = ndims(system) + 1 +@inline v_nvariables(system, ::TransportVelocityAdami) = ndims(system) * 2 + 1 + +@inline function add_velocity!(du, v, particle, system::EntropicallyDampedSPHSystem) + add_velocity!(du, v, particle, system, system.transport_velocity) +end + +# Add momentum velocity. +@inline function add_velocity!(du, v, particle, system, ::Nothing) + for i in 1:ndims(system) + du[i, particle] = v[i, particle] + end + + return du +end + +# Add advection velocity. +@inline function add_velocity!(du, v, particle, system, ::TransportVelocityAdami) + for i in 1:ndims(system) + du[i, particle] = v[ndims(system) + i, particle] + end + + return du +end + +@inline function advection_velocity(v, system::EntropicallyDampedSPHSystem, particle) + (; cache) = system + extract_svector(cache.advection_velocity, system, particle) +end function update_quantities!(system::EntropicallyDampedSPHSystem, system_index, v, u, v_ode, u_ode, semi, t) @@ -196,12 +230,29 @@ function update_average_pressure!(system, ::TransportVelocityAdami, system_index end function write_v0!(v0, system::EntropicallyDampedSPHSystem) - (; initial_condition) = system + write_v0!(v0, system, system.transport_velocity) +end +function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::Nothing) for particle in eachparticle(system) # Write particle velocities + v_init = initial_velocity(system, particle) for dim in 1:ndims(system) - v0[dim, particle] = initial_condition.velocity[dim, particle] + v0[dim, particle] = v_init[dim] + end + v0[end, particle] = initial_pressure(system, particle) + end + + return v0 +end + +function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::TransportVelocityAdami) + for particle in eachparticle(system) + # Write particle velocities + v_init = initial_velocity(system, particle) + for dim in 1:ndims(system) + v0[dim, particle] = v_init[dim] + v0[ndims(system) + dim, particle] = v_init[dim] end v0[end, particle] = initial_pressure(system, particle) end @@ -229,3 +280,18 @@ end particle_position = initial_coords(system, particle) return initial_pressure_function(particle_position) end + +@inline function initial_velocity(system, particle) + initial_velocity(system, particle, system.initial_velocity_function) +end + +@inline function initial_velocity(system, particle, ::Nothing) + return extract_svector(system.initial_condition.velocity, system, particle) +end + +@inline function initial_velocity(system, particle, init_velocity_function) + position = initial_coords(system, particle) + v_init = SVector(ntuple(i -> init_velocity_function[i](position), Val(ndims(system)))) + + return v_init +end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 258cb19f3..1085a0117 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -18,5 +18,6 @@ function write_u0!(u0, system::FluidSystem) end include("viscosity.jl") +include("transport_velocity.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl new file mode 100644 index 000000000..e604ad856 --- /dev/null +++ b/src/schemes/fluid/transport_velocity.jl @@ -0,0 +1,30 @@ +struct TransportVelocityAdami{ELTYPE} + background_pressure::ELTYPE + function TransportVelocityAdami(background_pressure) + new{typeof(background_pressure)}(background_pressure) + end +end + +@inline function update_transport_velocity!(system::FluidSystem, system_index, + v_ode, u_ode, semi) + update_transport_velocity!(system, system_index, v_ode, u_ode, semi, + system.transport_velocity) +end + +@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi, + ::TransportVelocityAdami) + (; cache) = system + v = wrap_v(v_ode, system_index, system, semi) + + for particle in each_moving_particle(system) + for i in 1:ndims(system) + cache.advection_velocity[i, particle] = v[ndims(system) + i, particle] + v[ndims(system) + i, particle] = v[i, particle] + end + end +end + +@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi, + ::Nothing) + return system +end From 69099c38379b39384a05914de8e0771738a5df7a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 15:54:48 +0200 Subject: [PATCH 03/62] add tvf in rhs --- .../fluid/entropically_damped_sph/rhs.jl | 68 ++++++++++++++++++- 1 file changed, 67 insertions(+), 1 deletion(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 019d0e878..392ec7cbf 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -42,8 +42,13 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b) + dv_convection = momentum_convection(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + rho_a, rho_b, particle, neighbor, grad_kernel, + volume_term) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity[i] + dv[i, particle] += dv_pressure[i] + dv_viscosity[i] + dv_convection[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - @@ -52,6 +57,8 @@ function interact!(dv, v_particle_system, u_particle_system, pressure_evolution!(dv, particle_system, v_diff, grad_kernel, particle, pos_diff, distance, sound_speed, volume_term, m_b, p_a, p_b, rho_a, rho_b) + + transport_velocity!(dv, particle_system, volume_term, grad_kernel, particle) end return dv @@ -81,6 +88,65 @@ end return dv end +@inline function momentum_convection(system, neighbor_system, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) +end + +@inline function momentum_convection(system::EntropicallyDampedSPHSystem, + neighbor_system::EntropicallyDampedSPHSystem, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + momentum_convection(system, neighbor_system, system.transport_velocity, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) +end + +@inline function momentum_convection(system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) +end + +@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + 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) + + A_a = rho_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' + A_b = rho_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' + + return volume_term * (0.5 * (A_a + A_b)) * grad_kernel +end + +@inline function transport_velocity!(dv, system::EntropicallyDampedSPHSystem, volume_term, + grad_kernel, particle) + transport_velocity!(dv, system, system.transport_velocity, volume_term, grad_kernel, + particle) +end + +@inline transport_velocity!(dv, system, volume_term, grad_kernel, particle) = dv + +@inline transport_velocity!(dv, system, ::Nothing, volume_term, grad_kernel, particle) = dv + +@inline function transport_velocity!(dv, system, ::TransportVelocityAdami, volume_term, + grad_kernel, particle) + (; transport_velocity) = system + (; background_pressure) = transport_velocity + n_dims = ndims(system) + + for dim in 1:n_dims + dv[n_dims + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] + end + + return dv +end + @inline average_pressure(system, particle) = 0.0 @inline viscosity_function(system) = system.viscosity From 189628811cfd87e633e3c327d25fbe869ef1c7d2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 16:03:38 +0200 Subject: [PATCH 04/62] add update callback --- src/TrixiParticles.jl | 2 +- src/callbacks/callbacks.jl | 1 + src/callbacks/update.jl | 25 +++++++++++++++++++++++++ 3 files changed, 27 insertions(+), 1 deletion(-) create mode 100644 src/callbacks/update.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c01f9f36d..22b381bfa 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -34,7 +34,7 @@ export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, BoundarySPHSystem -export InfoCallback, SolutionSavingCallback +export InfoCallback, SolutionSavingCallback, UpdateAfterTimeStep export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 20da95052..741f6ec70 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -1,2 +1,3 @@ include("info.jl") include("solution_saving.jl") +include("update.jl") diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl new file mode 100644 index 000000000..be6c75f6e --- /dev/null +++ b/src/callbacks/update.jl @@ -0,0 +1,25 @@ +# TODO: naming +function UpdateAfterTimeStep() + return DiscreteCallback(update_condition, update_after_dt!, initialize=initial_update!, + save_positions=(false, false)) +end + +# condition +update_condition(u, t, integrator) = true + +# affect +function update_after_dt!(integrator) + semi = integrator.p + v_ode, u_ode = integrator.u.x + + foreach_enumerate(semi.systems) do (system_index, system) + update_transport_velocity!(system, system_index, v_ode, u_ode, semi) + end + + return integrator +end + +# initialize +initial_update!(cb, u, t, integrator) = update_after_dt!(integrator) + +update_transport_velocity!(system, system_index, v_ode, u_ode, semi) = system From 5a4cf21d6b3baac5bdeb42a9f63ef5a5d97ffa31 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 16:45:49 +0200 Subject: [PATCH 05/62] use `step_limiter!` instead of using a callback --- src/TrixiParticles.jl | 2 +- src/callbacks/callbacks.jl | 1 - src/callbacks/update.jl | 25 ------------------- .../fluid/entropically_damped_sph/system.jl | 1 + src/schemes/fluid/transport_velocity.jl | 11 ++++++++ 5 files changed, 13 insertions(+), 27 deletions(-) delete mode 100644 src/callbacks/update.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 22b381bfa..c01f9f36d 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -34,7 +34,7 @@ export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, BoundarySPHSystem -export InfoCallback, SolutionSavingCallback, UpdateAfterTimeStep +export InfoCallback, SolutionSavingCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 741f6ec70..20da95052 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -1,3 +1,2 @@ include("info.jl") include("solution_saving.jl") -include("update.jl") diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl deleted file mode 100644 index be6c75f6e..000000000 --- a/src/callbacks/update.jl +++ /dev/null @@ -1,25 +0,0 @@ -# TODO: naming -function UpdateAfterTimeStep() - return DiscreteCallback(update_condition, update_after_dt!, initialize=initial_update!, - save_positions=(false, false)) -end - -# condition -update_condition(u, t, integrator) = true - -# affect -function update_after_dt!(integrator) - semi = integrator.p - v_ode, u_ode = integrator.u.x - - foreach_enumerate(semi.systems) do (system_index, system) - update_transport_velocity!(system, system_index, v_ode, u_ode, semi) - end - - return integrator -end - -# initialize -initial_update!(cb, u, t, integrator) = update_after_dt!(integrator) - -update_transport_velocity!(system, system_index, v_ode, u_ode, semi) = system diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index bc5c12001..d04e49e95 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -166,6 +166,7 @@ end @inline function add_velocity!(du, v, particle, system, ::TransportVelocityAdami) for i in 1:ndims(system) du[i, particle] = v[ndims(system) + i, particle] + system.cache.advection_velocity[i, particle] = v[ndims(system) + i, particle] end return du diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index e604ad856..fd072b391 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -5,6 +5,17 @@ struct TransportVelocityAdami{ELTYPE} end end +function update_transport_velocity!(u, integrator, p, t) + semi = integrator.p + v_ode, u_ode = integrator.u.x + + foreach_enumerate(semi.systems) do (system_index, system) + update_transport_velocity!(system, system_index, v_ode, u_ode, semi) + end + + return integrator +end + @inline function update_transport_velocity!(system::FluidSystem, system_index, v_ode, u_ode, semi) update_transport_velocity!(system, system_index, v_ode, u_ode, semi, From af1edcb4204e9841807cf2121258b48f62077e01 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 17:34:35 +0200 Subject: [PATCH 06/62] add taylor green vortex example --- examples/fluid/taylor_green_vortex_2d.jl | 138 +++++++++++++++++++++++ src/schemes/fluid/transport_velocity.jl | 7 +- 2 files changed, 141 insertions(+), 4 deletions(-) create mode 100644 examples/fluid/taylor_green_vortex_2d.jl diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl new file mode 100644 index 000000000..7246bd0ac --- /dev/null +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -0,0 +1,138 @@ +using TrixiParticles +using OrdinaryDiffEq +using Random, LinearAlgebra, StaticArrays +Random.seed!(42); + +# ========================================================================================== +# ==== Fluid + +particle_spacing = 0.02 + +reynolds_number = 100.0 + +box_length = 1.0 + +n_particles_xy = round(Int, box_length / particle_spacing) + +U = 1.0 # m/s + +b = -8pi^2 / reynolds_number + +fluid_density = 1.0 +sound_speed = 10U + +nu = U * box_length / reynolds_number + +background_pressure = sound_speed^2 * fluid_density + +p(pos) = -U^2 * exp(2 * b * 0) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 +p(pos, t) = -U^2 * exp(2 * b * t) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 + +v_x(pos) = -U * exp(b * 0) * cos(2pi * pos[1]) * sin(2pi * pos[2]) +v_y(pos) = U * exp(b * 0) * sin(2pi * pos[1]) * cos(2pi * pos[2]) +v_x(pos, t) = -U * exp(b * t) * cos(2pi * pos[1]) * sin(2pi * pos[2]) +v_y(pos, t) = U * exp(b * t) * sin(2pi * pos[1]) * cos(2pi * pos[2]) + +smoothing_length = 1.0 * particle_spacing +smoothing_kernel = SchoenbergQuinticSplineKernel{2}() + +fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0), + fluid_density, pressure=background_pressure, + init_velocity=(0.0, 0.0)) + +# Add small random displacement to the particles to avoid stagnant streamlines. +fluid.coordinates .+= rand((-particle_spacing / 5):1e-5:(particle_spacing / 5), + size(fluid.coordinates)) + +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + sound_speed, initial_pressure_function=p, + initial_velocity_function=(v_x, v_y), + transport_velocity=TransportVelocityAdami(background_pressure), + viscosity=ViscosityAdami(nu)) +semi = Semidiscretization(fluid_system, + neighborhood_search=GridNeighborhoodSearch, + periodic_box_min_corner=[0.0, 0.0], + periodic_box_max_corner=[box_length, box_length]) + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu)) + +function compute_L1v_error(v, u, t, system) + (; initial_velocity_function) = system + + v_analytical_avg = 0.0 + L1v = 0.0 + + for particle in TrixiParticles.eachparticle(system) + position = TrixiParticles.current_coords(u, system, particle) + + v_mag = norm(TrixiParticles.current_velocity(v, system, particle)) + v_analyitcal = norm(SVector(ntuple(i -> initial_velocity_function[i](position, t), + Val(ndims(system))))) + + v_analytical_avg += abs(v_analyitcal) + L1v += abs(v_mag - v_analyitcal) + end + v_analytical_avg /= TrixiParticles.nparticles(system) + + L1v /= TrixiParticles.nparticles(system) + + return L1v /= v_analytical_avg +end + +function compute_L1p_error(v, u, t, system) + p_max_exact = 0.0 + + L1p = 0.0 + + for particle in TrixiParticles.eachparticle(system) + position = TrixiParticles.current_coords(u, system, particle) + + # compute pressure error + p_analytical = system.initial_pressure_function(position, t) + p_max_exact = max(p_max_exact, abs(p_analytical)) + + # p_computed - p_average + p_computed = TrixiParticles.particle_pressure(v, system, particle) - + TrixiParticles.average_pressure(system, particle) + L1p += abs(p_computed - p_analytical) + end + + L1p /= TrixiParticles.nparticles(system) + + return L1p /= p_max_exact +end + +function average_pressure(v, u, t, system) + return [TrixiParticles.average_pressure(system, i) for i in axes(v, 2)] +end + +function time_vector(v, u, t, system) + return t +end + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, + L1v=compute_L1v_error, + L1p=compute_L1p_error, + p_avg=average_pressure, + t=time_vector) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Enable threading of the RK method for better performance on multiple threads. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces +# become extremely large when fluid particles are very close to boundary particles, +# and the time integration method interprets this as an instability. +sol = solve(ode, + RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!), + abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=dt_max,#1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index fd072b391..ec62da8a4 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -5,15 +5,14 @@ struct TransportVelocityAdami{ELTYPE} end end -function update_transport_velocity!(u, integrator, p, t) - semi = integrator.p - v_ode, u_ode = integrator.u.x +function update_transport_velocity!(u, ode, semi, t) + v_ode, u_ode = ode.u.x foreach_enumerate(semi.systems) do (system_index, system) update_transport_velocity!(system, system_index, v_ode, u_ode, semi) end - return integrator + return ode end @inline function update_transport_velocity!(system::FluidSystem, system_index, From 91eefc2dce1e249c1e38ee30e81ada8c9cc2d5a0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 7 Sep 2023 23:00:13 +0200 Subject: [PATCH 07/62] add lid driven cavity example --- examples/fluid/lid_driven_cavity_2d.jl | 108 ++++++++++++++++++++++++ src/schemes/fluid/transport_velocity.jl | 14 +-- 2 files changed, 117 insertions(+), 5 deletions(-) create mode 100644 examples/fluid/lid_driven_cavity_2d.jl diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl new file mode 100644 index 000000000..78d297bf0 --- /dev/null +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -0,0 +1,108 @@ +using TrixiParticles +using OrdinaryDiffEq + +gravity = 0.0 +reynolds_number = 100.0 +particle_spacing = 0.02 + +# ========================================================================================== +# ==== Fluid + +water_width = 1.0 +water_height = 1.0 +water_density = 1.0 + +tank_width = 1.0 +tank_height = 1.0 + +velocity_lid = 1.0 +sound_speed = 10 * velocity_lid + +pressure = sound_speed^2 * water_density + +smoothing_length = 1.0 * particle_spacing +smoothing_kernel = SchoenbergQuinticSplineKernel{2}() + +nu = velocity_lid / reynolds_number +viscosity = ViscosityAdami(nu) + +tank = RectangularTank(particle_spacing, (water_width, water_height), + (tank_width, tank_height), water_density, + n_layers=3, spacing_ratio=1, faces=(true, true, true, false), + pressure=pressure) + +# ========================================================================================== +# ==== Lid + +lid_position = 0.0 - particle_spacing * 4 +lid_length = tank.n_particles_per_dimension[1] + 8 + +lid = RectangularShape(particle_spacing, (lid_length, 3), + (lid_position, water_height), water_density) + +f_y(t) = 0.0 +f_x(t) = velocity_lid * t + +is_moving(t) = true + +movement = BoundaryMovement((f_x, f_y), is_moving) + +# ========================================================================================== +# ==== Boundary models + +boundary_model_tank = BoundaryModelDummyParticles(tank.boundary.density, + tank.boundary.mass, + AdamiPressureExtrapolation(), + viscosity=viscosity, + smoothing_kernel, smoothing_length) + +boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass, + AdamiPressureExtrapolation(), + viscosity=viscosity, + smoothing_kernel, smoothing_length) + +# ========================================================================================== +# ==== Systems + +fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, + smoothing_length, + sound_speed, + viscosity=viscosity, + transport_velocity=TransportVelocityAdami(pressure), + acceleration=(0.0, gravity)) + +boundary_system_tank = BoundarySPHSystem(tank.boundary, boundary_model_tank) + +boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=movement) + +# ========================================================================================== +# ==== Simulation + +tspan = (0.0, 5.0) + +semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_lid, + neighborhood_search=GridNeighborhoodSearch, + periodic_box_min_corner=[0.0 - 4 * particle_spacing, -0.24], + periodic_box_max_corner=[1.0 + 4 * particle_spacing, + 1.0 + 4 * particle_spacing]) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) + +saving_callback = SolutionSavingCallback(dt=0.02) +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Enable threading of the RK method for better performance on multiple threads. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces +# become extremely large when fluid particles are very close to boundary particles, +# and the time integration method interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!), + abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index ec62da8a4..2eb9ec949 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -15,12 +15,21 @@ function update_transport_velocity!(u, ode, semi, t) return ode end +@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi) + return system +end + @inline function update_transport_velocity!(system::FluidSystem, system_index, v_ode, u_ode, semi) update_transport_velocity!(system, system_index, v_ode, u_ode, semi, system.transport_velocity) end +@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi, + ::Nothing) + return system +end + @inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi, ::TransportVelocityAdami) (; cache) = system @@ -33,8 +42,3 @@ end end end end - -@inline function update_transport_velocity!(system, system_index, v_ode, u_ode, semi, - ::Nothing) - return system -end From 2a1b7938584e36e82286845d10e29df0660760ac Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 9 Sep 2023 12:44:09 +0200 Subject: [PATCH 08/62] add periodic array of cylinder example --- examples/fluid/lid_driven_cavity_2d.jl | 7 ++ .../fluid/periodic_array_of_cylinders_2d.jl | 84 +++++++++++++++++++ examples/fluid/taylor_green_vortex_2d.jl | 7 ++ 3 files changed, 98 insertions(+) create mode 100644 examples/fluid/periodic_array_of_cylinders_2d.jl diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index 78d297bf0..72cf3bc0d 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -1,3 +1,10 @@ +# Lid-driven cavity +# +# S. Adami et al +# "A transport-velocity formulation for smoothed particle hydrodynamics". +# In: Journal of Computational Physics, Volume 241 (2013), pages 292-307. +# https://doi.org/10.1016/j.jcp.2013.01.043 + using TrixiParticles using OrdinaryDiffEq diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl new file mode 100644 index 000000000..52c37866a --- /dev/null +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -0,0 +1,84 @@ +# Channel flow through periodic array of cylinders +# +# S. Adami et al +# "A transport-velocity formulation for smoothed particle hydrodynamics". +# In: Journal of Computational Physics, Volume 241 (2013), pages 292-307. +# https://doi.org/10.1016/j.jcp.2013.01.043 + +using TrixiParticles +using OrdinaryDiffEq + +acceleration_x = 2.5e-4 + +# ========================================================================================== +# ==== Fluid + +n_particles_x = 144 + +cylinder_radius = 0.02 +box_length = 6cylinder_radius +box_height = 4cylinder_radius + +fluid_density = 1000.0 +nu = 0.1 / fluid_density + +velocity_ref = 1.2e-4 + +# Adami uses c = 0.1*sqrt(acceleration_x*cylinder_radius) but from the original setup +# from M. Ellero and N. A. Adams (https://doi.org/10.1002/nme.3088) is using c = 0.02 +sound_speed = 0.02 + +pressure = sound_speed^2 * fluid_density + +particle_spacing = box_length/n_particles_x +smoothing_length = 1.2 * particle_spacing +smoothing_kernel = SchoenbergQuarticSplineKernel{2}() + +box = RectangularTank(particle_spacing, (box_length, box_height), (box_length, box_height), + fluid_density, n_layers=3, spacing_ratio=1.0, pressure=pressure, + faces=(false, false, true, true)) + +cylinder = SphereShape(particle_spacing, cylinder_radius, (box_length / 2, box_height / 2), + fluid_density, sphere_type=RoundSphere()) + +fluid = setdiff(box.fluid, cylinder) +boundary = union(cylinder, box.boundary) + +boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass, + AdamiPressureExtrapolation(), + viscosity=ViscosityAdami(nu), + smoothing_kernel, smoothing_length) + +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=ViscosityAdami(nu), + transport_velocity=TransportVelocityAdami(pressure), + acceleration=(acceleration_x, 0.0)) + +boundary_system = BoundarySPHSystem(boundary, boundary_model) + +semi = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch, + periodic_box_min_corner=[0.0, -box_length * 2], + periodic_box_max_corner=[box_length, box_length * 2]) + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=0.02) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Enable threading of the RK method for better performance on multiple threads. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces +# become extremely large when fluid particles are very close to boundary particles, +# and the time integration method interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!), + abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 7246bd0ac..198beefdf 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -1,3 +1,10 @@ +# Taylor Green vortex +# +# P. Ramachandran, K. Puri +# "Entropically damped artificial compressibility for SPH". +# In: Computers and Fluids, Volume 179 (2019), pages 579-594. +# https://doi.org/10.1016/j.compfluid.2018.11.023 + using TrixiParticles using OrdinaryDiffEq using Random, LinearAlgebra, StaticArrays From 0d9c9321e811884013b92e79e72778540c5e6fd6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 11 Sep 2023 11:00:34 +0200 Subject: [PATCH 09/62] add tests --- .../fluid/periodic_array_of_cylinders_2d.jl | 2 +- examples/fluid/taylor_green_vortex_2d.jl | 6 ++--- .../fluid/entropically_damped_sph/system.jl | 10 ++++++-- test/examples/examples.jl | 24 +++++++++++++++++++ test/systems/edac_system.jl | 8 +++++++ 5 files changed, 44 insertions(+), 6 deletions(-) diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 52c37866a..6a49fa57b 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -30,7 +30,7 @@ sound_speed = 0.02 pressure = sound_speed^2 * fluid_density -particle_spacing = box_length/n_particles_x +particle_spacing = box_length / n_particles_x smoothing_length = 1.2 * particle_spacing smoothing_kernel = SchoenbergQuarticSplineKernel{2}() diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 198beefdf..78c48ea7c 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -76,11 +76,11 @@ function compute_L1v_error(v, u, t, system) position = TrixiParticles.current_coords(u, system, particle) v_mag = norm(TrixiParticles.current_velocity(v, system, particle)) - v_analyitcal = norm(SVector(ntuple(i -> initial_velocity_function[i](position, t), + v_analytical = norm(SVector(ntuple(i -> initial_velocity_function[i](position, t), Val(ndims(system))))) - v_analytical_avg += abs(v_analyitcal) - L1v += abs(v_mag - v_analyitcal) + v_analytical_avg += abs(v_analytical) + L1v += abs(v_mag - v_analytical) end v_analytical_avg /= TrixiParticles.nparticles(system) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index d04e49e95..2f1b15485 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -1,6 +1,10 @@ @doc raw""" - EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, - sound_speed; alpha=0.5, viscosity=NoViscosity(), + EntropicallyDampedSPHSystem(initial_condition, + smoothing_kernel, smoothing_length, sound_speed; + alpha=0.5, viscosity=NoViscosity(), + initial_pressure_function=nothing, + initial_velocity_function=nothing, + transport_velocity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS)) Entropically damped artificial compressibility (EDAC) for SPH introduced by (Ramachandran 2019). @@ -109,6 +113,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) + summary_line(io, "initial pressure function", system.initial_pressure_function) + summary_line(io, "initial velocity function", system.initial_velocity_function) summary_line(io, "acceleration", system.acceleration) summary_footer(io) end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index cae568d7e..d201ce518 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -54,6 +54,30 @@ tspan=(0.0, 0.4)) @test sol.retcode == ReturnCode.Success end + + @trixi_testset "fluid/lid_driven_cavity_2d.jl" begin + @test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "lid_driven_cavity_2d.jl"), + tspan=(0.0, 0.1)) + @test sol.retcode == ReturnCode.Success + end + + @trixi_testset "fluid/taylor_green_vortex_2d.jl" begin + @test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "taylor_green_vortex_2d.jl"), + tspan=(0.0, 0.1)) + @test sol.retcode == ReturnCode.Success + end + + @trixi_testset "fluid/periodic_array_of_cylinders_2d.jl" begin + @test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "periodic_array_of_cylinders_2d.jl"), + tspan=(0.0, 0.1)) + @test sol.retcode == ReturnCode.Success + end end @testset verbose=true "Solid" begin diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 33c1f8110..7f4e93d60 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -33,6 +33,9 @@ @test system.smoothing_kernel == smoothing_kernel @test system.smoothing_length == smoothing_length @test system.viscosity isa TrixiParticles.NoViscosity + @test system.initial_pressure_function isa Nothing + @test system.initial_velocity_function isa Nothing + @test system.transport_velocity isa Nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -86,6 +89,9 @@ @test system.smoothing_kernel == smoothing_kernel @test system.smoothing_length == smoothing_length @test system.viscosity isa TrixiParticles.NoViscosity + @test system.initial_pressure_function isa Nothing + @test system.initial_velocity_function isa Nothing + @test system.transport_velocity isa Nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @test length(system.mass) == size(setup.coordinates, 2) @@ -135,6 +141,8 @@ │ viscosity: …………………………………………………… NoViscosity │ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ + │ initial pressure function: ………… nothing │ + │ initial velocity function: ………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box From 7e564c3e15320ee5db495a38914f74863d71bc2e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 11 Sep 2023 15:40:15 +0200 Subject: [PATCH 10/62] add Random package --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index d3477f221..54d9001c6 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Morton = "2a6d852e-3fac-5a38-885c-fe708af2d09e" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" From 31d59b926046b655e7c1fc82c99d4ba65875969f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 11 Sep 2023 17:58:00 +0200 Subject: [PATCH 11/62] reexport `seed!` --- examples/fluid/taylor_green_vortex_2d.jl | 3 +-- src/TrixiParticles.jl | 4 +++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 78c48ea7c..288ff26bf 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -7,8 +7,6 @@ using TrixiParticles using OrdinaryDiffEq -using Random, LinearAlgebra, StaticArrays -Random.seed!(42); # ========================================================================================== # ==== Fluid @@ -48,6 +46,7 @@ fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0. init_velocity=(0.0, 0.0)) # Add small random displacement to the particles to avoid stagnant streamlines. +seed!(42); fluid.coordinates .+= rand((-particle_spacing / 5):1e-5:(particle_spacing / 5), size(fluid.coordinates)) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c01f9f36d..095360712 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -4,13 +4,15 @@ using Reexport: @reexport using Dates using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect -using LinearAlgebra: norm, dot, I, tr +using LinearAlgebra: dot, I, tr +@reexport using LinearAlgebra: norm using Morton: cartesian2morton using Polyester: Polyester, @batch using Printf: @printf using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, get_tmp_cache @reexport using StaticArrays: SVector +@reexport using Random: seed! using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt using ThreadingUtilities From 1ce89bb249dcbe9eeb429e6da2d733eea64e88dc Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 23 Nov 2023 15:33:48 +0100 Subject: [PATCH 12/62] remove obsolet `system_index` --- src/schemes/fluid/entropically_damped_sph/system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index ff7f35586..74d595da5 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -189,7 +189,7 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi) end -function update_average_pressure!(system, ::Nothing, system_index, v_ode, u_ode, semi) +function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi) return system end From fc2ad8ebcf47d17d5ac1397e6110bedcf6449b4f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 23 Nov 2023 15:40:54 +0100 Subject: [PATCH 13/62] skip CI --- examples/fluid/taylor_green_vortex_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 4c0cf49ab..38497b7cf 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -58,7 +58,7 @@ fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_le sound_speed, initial_pressure_function=p, initial_velocity_function=(v_x, v_y), transport_velocity=TransportVelocityAdami(background_pressure), - viscosity=ViscosityAdami(;nu)) + viscosity=ViscosityAdami(; nu)) # ========================================================================================== # ==== Simulation From 1e4d10b2dd99fa906b2a62f77ec02552025d36dc Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sun, 26 Nov 2023 12:14:23 +0100 Subject: [PATCH 14/62] add callback --- examples/fluid/taylor_green_vortex_2d.jl | 14 ++++---- src/TrixiParticles.jl | 3 +- src/callbacks/callbacks.jl | 1 + src/callbacks/update.jl | 44 ++++++++++++++++++++++++ 4 files changed, 53 insertions(+), 9 deletions(-) create mode 100644 src/callbacks/update.jl diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 38497b7cf..188346a48 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -19,9 +19,9 @@ reynolds_number = 100.0 box_length = 1.0 -n_particles_xy = round(Int, box_length / particle_spacing) - U = 1.0 # m/s +fluid_density = 1.0 +sound_speed = 10U b = -8pi^2 / reynolds_number @@ -33,11 +33,10 @@ v_y(pos) = U * exp(b * 0) * sin(2pi * pos[1]) * cos(2pi * pos[2]) v_x(pos, t) = -U * exp(b * t) * cos(2pi * pos[1]) * sin(2pi * pos[2]) v_y(pos, t) = U * exp(b * t) * sin(2pi * pos[1]) * cos(2pi * pos[2]) +n_particles_xy = round(Int, box_length / particle_spacing) + # ========================================================================================== # ==== Fluid -fluid_density = 1.0 -sound_speed = 10U - nu = U * box_length / reynolds_number background_pressure = sound_speed^2 * fluid_density @@ -133,7 +132,7 @@ saving_callback = SolutionSavingCallback(dt=0.02, p_avg=average_pressure, t=time_vector) -callbacks = CallbackSet(info_callback, saving_callback) +callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep()) # Use a Runge-Kutta method with automatic (error based) time step size control. # Enable threading of the RK method for better performance on multiple threads. @@ -143,8 +142,7 @@ callbacks = CallbackSet(info_callback, saving_callback) # Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces # become extremely large when fluid particles are very close to boundary particles, # and the time integration method interprets this as an instability. -sol = solve(ode, - RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!), +sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) dtmax=dt_max,#1e-2, # Limit stepsize to prevent crashing diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index ea3583ac6..3558a52ac 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -39,7 +39,8 @@ export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, BoundarySPHSystem -export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback +export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, + UpdateEachTimeStep export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 623479592..09b7f7437 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -8,3 +8,4 @@ end include("info.jl") include("solution_saving.jl") include("density_reinit.jl") +include("update.jl") diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl new file mode 100644 index 000000000..9bfcb067c --- /dev/null +++ b/src/callbacks/update.jl @@ -0,0 +1,44 @@ +# TODO: naming +struct UpdateEachDt end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:UpdateEachDt}) + @nospecialize cb # reduce precompilation time + print(io, "UpdateEachTimeStep()") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:UpdateEachDt}) + @nospecialize cb # reduce precompilation time + if get(io, :compact, false) + show(io, cb) + else + summary_box(io, "UpdateEachTimeStep") + end +end + +function UpdateEachTimeStep() + update_each_dt = UpdateEachDt() + return DiscreteCallback(update_each_dt, update_each_dt, initialize=initial_update!, + save_positions=(false, false)) +end + +# condition +(update_each_dt::UpdateEachDt)(u, t, integrator) = true + +# affect +function (update_each_dt::UpdateEachDt)(integrator) + semi = integrator.p + v_ode, u_ode = integrator.u.x + + foreach_system(semi) do system + update_transport_velocity!(system, v_ode, semi) + end + + # Tell OrdinaryDiffEq that u has been modified + u_modified!(integrator, true) + + return integrator +end + +# initialize +initial_update!(cb, u, t, integrator) = cb.affect!(integrator) From d7a573ab959292f150e41b47558106cc28ffaa85 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sun, 26 Nov 2023 12:35:37 +0100 Subject: [PATCH 15/62] rename edac cache --- src/schemes/fluid/entropically_damped_sph/system.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 74d595da5..94a8d9d72 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -80,7 +80,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, VF, TV, density_calculator = SummationDensity() - cache = create_cache(initial_condition, transport_velocity) + cache = create_cache_edac(initial_condition, transport_velocity) new{NDIMS, ELTYPE, typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(initial_pressure_function), @@ -120,9 +120,9 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -create_cache(initial_condition, ::Nothing) = (;) +create_cache_edac(initial_condition, ::Nothing) = (;) -function create_cache(initial_condition, ::TransportVelocityAdami) +function create_cache_edac(initial_condition, ::TransportVelocityAdami) pressure_average = copy(initial_condition.pressure) neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) advection_velocity = copy(initial_condition.velocity) From 0d956e6cc6a0573f3d662b33d2ae025e967b067e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 2 Jan 2024 16:41:31 +0100 Subject: [PATCH 16/62] move functions to `transport_velocity.jl` --- .../fluid/entropically_damped_sph/rhs.jl | 61 ------------------- .../fluid/entropically_damped_sph/system.jl | 2 + src/schemes/fluid/transport_velocity.jl | 59 ++++++++++++++++++ 3 files changed, 61 insertions(+), 61 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index ca3de307a..93190c7d5 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -88,64 +88,3 @@ end return dv end - -@inline function momentum_convection(system, neighbor_system, - v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) - return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) -end - -@inline function momentum_convection(system::EntropicallyDampedSPHSystem, - neighbor_system::EntropicallyDampedSPHSystem, - v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) - momentum_convection(system, neighbor_system, system.transport_velocity, - v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) -end - -@inline function momentum_convection(system, neighbor_system, ::Nothing, - v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) - return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) -end - -@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, - v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) - 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) - - A_a = rho_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' - A_b = rho_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' - - return volume_term * (0.5 * (A_a + A_b)) * grad_kernel -end - -@inline function transport_velocity!(dv, system::EntropicallyDampedSPHSystem, volume_term, - grad_kernel, particle) - transport_velocity!(dv, system, system.transport_velocity, volume_term, grad_kernel, - particle) -end - -@inline transport_velocity!(dv, system, volume_term, grad_kernel, particle) = dv - -@inline transport_velocity!(dv, system, ::Nothing, volume_term, grad_kernel, particle) = dv - -@inline function transport_velocity!(dv, system, ::TransportVelocityAdami, volume_term, - grad_kernel, particle) - (; transport_velocity) = system - (; background_pressure) = transport_velocity - n_dims = ndims(system) - - for dim in 1:n_dims - dv[n_dims + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] - end - - return dv -end - -@inline average_pressure(system, particle) = 0.0 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index a51f0cc10..af33b33b6 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -137,6 +137,8 @@ end return v[end, particle] end +@inline average_pressure(system, particle) = zero(eltype(system)) + @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) average_pressure(system, system.transport_velocity, particle) end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 8d2a16a4d..c744761a9 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -46,3 +46,62 @@ end @inline function advection_velocity(v, system, particle) return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) end + +@inline function momentum_convection(system, neighbor_system, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) +end + +@inline function momentum_convection(system::FluidSystem, + neighbor_system::FluidSystem, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + momentum_convection(system, neighbor_system, system.transport_velocity, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) +end + +@inline function momentum_convection(system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) +end + +@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + v_particle_system, v_neighbor_system, rho_a, rho_b, + particle, neighbor, grad_kernel, volume_term) + 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) + + A_a = rho_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' + A_b = rho_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' + + return volume_term * (0.5 * (A_a + A_b)) * grad_kernel +end + +@inline function transport_velocity!(dv, system::FluidSystem, volume_term, + grad_kernel, particle) + transport_velocity!(dv, system, system.transport_velocity, volume_term, grad_kernel, + particle) +end + +@inline transport_velocity!(dv, system, volume_term, grad_kernel, particle) = dv + +@inline transport_velocity!(dv, system, ::Nothing, volume_term, grad_kernel, particle) = dv + +@inline function transport_velocity!(dv, system, ::TransportVelocityAdami, volume_term, + grad_kernel, particle) + (; transport_velocity) = system + (; background_pressure) = transport_velocity + n_dims = ndims(system) + + for dim in 1:n_dims + dv[n_dims + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] + end + + return dv +end From edeec8df0b1c1ea8b487f1e8d4b2ade26906b8ad Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 2 Jan 2024 17:17:04 +0100 Subject: [PATCH 17/62] prepare for wcsph --- .../fluid/entropically_damped_sph/rhs.jl | 1 + .../fluid/entropically_damped_sph/system.jl | 50 +------------------ src/schemes/fluid/fluid.jl | 33 ++++++++++++ src/schemes/fluid/transport_velocity.jl | 35 ++++++++++--- 4 files changed, 65 insertions(+), 54 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 93190c7d5..00f8e1200 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -43,6 +43,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_mean) + # Add convection term when using `TransportVelocityAdami` dv_convection = momentum_convection(particle_system, neighbor_system, v_particle_system, v_neighbor_system, rho_a, rho_b, particle, neighbor, grad_kernel, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index af33b33b6..40ba090dd 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -150,14 +150,7 @@ end @inline average_pressure(system, ::Nothing, particle) = zero(eltype(system)) @inline function v_nvariables(system::EntropicallyDampedSPHSystem) - v_nvariables(system, system.transport_velocity) -end - -@inline v_nvariables(system, ::Nothing) = ndims(system) + 1 -@inline v_nvariables(system, ::TransportVelocityAdami) = ndims(system) * 2 + 1 - -@inline function add_velocity!(du, v, particle, system::EntropicallyDampedSPHSystem) - add_velocity!(du, v, particle, system, system.transport_velocity) + return ndims(system) * factor_tvf(system) + 1 end function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, @@ -192,7 +185,6 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, neighborhood_search) do particle, neighbor, pos_diff, distance - pressure_average[particle] += particle_pressure(v_neighbor_system, neighbor_system, neighbor) @@ -207,31 +199,8 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode end end -function write_v0!(v0, system::EntropicallyDampedSPHSystem) - write_v0!(v0, system, system.transport_velocity) -end - -function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::Nothing) - for particle in eachparticle(system) - # Write particle velocities - v_init = initial_velocity(system, particle) - for dim in 1:ndims(system) - v0[dim, particle] = v_init[dim] - end - v0[end, particle] = initial_pressure(system, particle) - end - - return v0 -end - -function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::TransportVelocityAdami) +function write_v0!(v0, density_calculator, system::EntropicallyDampedSPHSystem) for particle in eachparticle(system) - # Write particle velocities - v_init = initial_velocity(system, particle) - for dim in 1:ndims(system) - v0[dim, particle] = v_init[dim] - v0[ndims(system) + dim, particle] = v_init[dim] - end v0[end, particle] = initial_pressure(system, particle) end @@ -258,18 +227,3 @@ end particle_position = initial_coords(system, particle) return initial_pressure_function(particle_position) end - -@inline function initial_velocity(system, particle) - initial_velocity(system, particle, system.initial_velocity_function) -end - -@inline function initial_velocity(system, particle, ::Nothing) - return extract_svector(system.initial_condition.velocity, system, particle) -end - -@inline function initial_velocity(system, particle, init_velocity_function) - position = initial_coords(system, particle) - v_init = SVector(ntuple(i -> init_velocity_function[i](position), Val(ndims(system)))) - - return v_init -end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 1cf87cbdb..93762dc32 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -17,6 +17,39 @@ function write_u0!(u0, system::FluidSystem) return u0 end +function write_v0!(v0, system::FluidSystem) + write_v0!(v0, system, system.transport_velocity) +end + +function write_v0!(v0, system::FluidSystem, ::Nothing) + for particle in eachparticle(system) + # Write particle velocities + v_init = initial_velocity(system, particle) + for dim in 1:ndims(system) + v0[dim, particle] = v_init[dim] + end + end + + write_v0!(v0, system.density_calculator, system) + + return v0 +end + +@inline function initial_velocity(system, particle) + initial_velocity(system, particle, system.initial_velocity_function) +end + +@inline function initial_velocity(system, particle, ::Nothing) + return extract_svector(system.initial_condition.velocity, system, particle) +end + +@inline function initial_velocity(system, particle, init_velocity_function) + position = initial_coords(system, particle) + v_init = SVector(ntuple(i -> init_velocity_function[i](position), Val(ndims(system)))) + + return v_init +end + @inline viscosity_model(system::FluidSystem) = system.viscosity include("viscosity.jl") diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index c744761a9..cb14083ab 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -5,6 +5,11 @@ struct TransportVelocityAdami{ELTYPE} end end +# Calculate `v_nvariables` appropriately +@inline factor_tvf(system::FluidSystem) = factor_tvf(system, system.transport_velocity) +@inline factor_tvf(system, ::Nothing) = 1 +@inline factor_tvf(system, ::TransportVelocityAdami) = 2 + @inline update_transport_velocity!(system, v_ode, semi) = system @inline function update_transport_velocity!(system::FluidSystem, v_ode, semi) @@ -15,7 +20,6 @@ end @inline function update_transport_velocity!(system, v_ode, semi, ::TransportVelocityAdami) v = wrap_v(v_ode, system, semi) - for particle in each_moving_particle(system) for i in 1:ndims(system) v[ndims(system) + i, particle] = v[i, particle] @@ -25,6 +29,25 @@ end return system end +function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami) + for particle in eachparticle(system) + # Write particle velocities + v_init = initial_velocity(system, particle) + for dim in 1:ndims(system) + v0[dim, particle] = v_init[dim] + v0[ndims(system) + dim, particle] = v_init[dim] + end + end + + write_v0!(v0, system.density_calculator, system) + + return v0 +end + +@inline function add_velocity!(du, v, particle, system::FluidSystem) + add_velocity!(du, v, particle, system, system.transport_velocity) +end + # Add momentum velocity. @inline function add_velocity!(du, v, particle, system, ::Nothing) for i in 1:ndims(system) @@ -83,24 +106,24 @@ end return volume_term * (0.5 * (A_a + A_b)) * grad_kernel end +@inline transport_velocity!(dv, system, volume_term, grad_kernel, particle) = dv + @inline function transport_velocity!(dv, system::FluidSystem, volume_term, grad_kernel, particle) transport_velocity!(dv, system, system.transport_velocity, volume_term, grad_kernel, particle) end -@inline transport_velocity!(dv, system, volume_term, grad_kernel, particle) = dv - @inline transport_velocity!(dv, system, ::Nothing, volume_term, grad_kernel, particle) = dv @inline function transport_velocity!(dv, system, ::TransportVelocityAdami, volume_term, grad_kernel, particle) (; transport_velocity) = system (; background_pressure) = transport_velocity - n_dims = ndims(system) + NDIMS = ndims(system) - for dim in 1:n_dims - dv[n_dims + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] + for dim in 1:NDIMS + dv[NDIMS + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] end return dv From bc10ea38fbfb5d97d401a571c345318b189d2ff6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 3 Jan 2024 11:59:37 +0100 Subject: [PATCH 18/62] calculate volume term on the fly --- .../fluid/entropically_damped_sph/rhs.jl | 7 ++-- src/schemes/fluid/transport_velocity.jl | 38 ++++++++++++------- 2 files changed, 29 insertions(+), 16 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 00f8e1200..22a22e71c 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -46,8 +46,8 @@ function interact!(dv, v_particle_system, u_particle_system, # Add convection term when using `TransportVelocityAdami` dv_convection = momentum_convection(particle_system, neighbor_system, v_particle_system, v_neighbor_system, - rho_a, rho_b, particle, neighbor, grad_kernel, - volume_term) + 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] @@ -60,7 +60,8 @@ function interact!(dv, v_particle_system, u_particle_system, particle, pos_diff, distance, sound_speed, volume_term, m_b, p_a, p_b, rho_a, rho_b) - transport_velocity!(dv, particle_system, volume_term, grad_kernel, particle) + transport_velocity!(dv, particle_system, rho_a, rho_b, m_a, m_b, + grad_kernel, particle) end return dv diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index cb14083ab..0a291d624 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -72,28 +72,32 @@ end @inline function momentum_convection(system, neighbor_system, v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) + m_a, m_b, particle, neighbor, grad_kernel) return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) end @inline function momentum_convection(system::FluidSystem, neighbor_system::FluidSystem, v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) + 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, - particle, neighbor, grad_kernel, volume_term) + m_a, m_b, particle, neighbor, grad_kernel) end @inline function momentum_convection(system, neighbor_system, ::Nothing, v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) + 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, v_particle_system, v_neighbor_system, rho_a, rho_b, - particle, neighbor, grad_kernel, volume_term) + m_a, m_b, particle, neighbor, grad_kernel) + volume_a = m_a / rho_a + volume_b = m_b / rho_b + volume_term = (volume_a^2 + volume_b^2) / m_a + momentum_velocity_a = current_velocity(v_particle_system, system, particle) advection_velocity_a = advection_velocity(v_particle_system, system, particle) @@ -106,22 +110,30 @@ end return volume_term * (0.5 * (A_a + A_b)) * grad_kernel end -@inline transport_velocity!(dv, system, volume_term, grad_kernel, particle) = dv +@inline transport_velocity!(dv, system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) = dv -@inline function transport_velocity!(dv, system::FluidSystem, volume_term, - grad_kernel, particle) - transport_velocity!(dv, system, system.transport_velocity, volume_term, grad_kernel, - particle) +@inline function transport_velocity!(dv, system::FluidSystem, + rho_a, rho_b, m_a, m_b, grad_kernel, particle) + transport_velocity!(dv, system, system.transport_velocity, rho_a, rho_b, m_a, m_b, + grad_kernel, particle) end -@inline transport_velocity!(dv, system, ::Nothing, volume_term, grad_kernel, particle) = dv +@inline function transport_velocity!(dv, system, ::Nothing, + rho_a, rho_b, m_a, m_b, grad_kernel, particle) + return dv +end -@inline function transport_velocity!(dv, system, ::TransportVelocityAdami, volume_term, - grad_kernel, particle) +@inline function transport_velocity!(dv, system, ::TransportVelocityAdami, + rho_a, rho_b, m_a, m_b, grad_kernel, particle) (; transport_velocity) = system (; background_pressure) = transport_velocity + NDIMS = ndims(system) + volume_a = m_a / rho_a + volume_b = m_b / rho_b + volume_term = (volume_a^2 + volume_b^2) / m_a + for dim in 1:NDIMS dv[NDIMS + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] end From a3ffb4090045e1ed7417898e4611f5fbfa7d1144 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 4 Jan 2024 16:16:59 +0100 Subject: [PATCH 19/62] adapt example files --- examples/fluid/lid_driven_cavity_2d.jl | 94 +++++++++---------- .../fluid/periodic_array_of_cylinders_2d.jl | 66 +++++++------ examples/fluid/rectangular_tank_edac_2d.jl | 5 - 3 files changed, 82 insertions(+), 83 deletions(-) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index 72cf3bc0d..f50480f57 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -7,45 +7,52 @@ using TrixiParticles using OrdinaryDiffEq - -gravity = 0.0 -reynolds_number = 100.0 +# ========================================================================================== +# ==== Resolution particle_spacing = 0.02 +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 3 +spacing_ratio = 1 + # ========================================================================================== -# ==== Fluid +# ==== Experiment Setup +tspan = (0.0, 5.0) +reynolds_number = 100.0 -water_width = 1.0 -water_height = 1.0 -water_density = 1.0 +cavity_size = (1.0, 1.0) -tank_width = 1.0 -tank_height = 1.0 +fluid_density = 1.0 velocity_lid = 1.0 sound_speed = 10 * velocity_lid -pressure = sound_speed^2 * water_density +pressure = sound_speed^2 * fluid_density -smoothing_length = 1.0 * particle_spacing -smoothing_kernel = SchoenbergQuinticSplineKernel{2}() +viscosity = ViscosityAdami(; nu=velocity_lid / reynolds_number) + +cavity = RectangularTank(particle_spacing, cavity_size, cavity_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), pressure=pressure) -nu = velocity_lid / reynolds_number -viscosity = ViscosityAdami(nu) +lid_position = 0.0 - particle_spacing * boundary_layers +lid_length = cavity.n_particles_per_dimension[1] + 2boundary_layers -tank = RectangularTank(particle_spacing, (water_width, water_height), - (tank_width, tank_height), water_density, - n_layers=3, spacing_ratio=1, faces=(true, true, true, false), - pressure=pressure) +lid = RectangularShape(particle_spacing, (lid_length, 3), + (lid_position, cavity_size[2]), fluid_density) # ========================================================================================== -# ==== Lid +# ==== Fluid -lid_position = 0.0 - particle_spacing * 4 -lid_length = tank.n_particles_per_dimension[1] + 8 +smoothing_length = 1.0 * particle_spacing +smoothing_kernel = SchoenbergQuinticSplineKernel{2}() -lid = RectangularShape(particle_spacing, (lid_length, 3), - (lid_position, water_height), water_density) +fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + transport_velocity=TransportVelocityAdami(pressure)) + +# ========================================================================================== +# ==== Boundary f_y(t) = 0.0 f_x(t) = velocity_lid * t @@ -54,51 +61,36 @@ is_moving(t) = true movement = BoundaryMovement((f_x, f_y), is_moving) -# ========================================================================================== -# ==== Boundary models - -boundary_model_tank = BoundaryModelDummyParticles(tank.boundary.density, - tank.boundary.mass, - AdamiPressureExtrapolation(), - viscosity=viscosity, - smoothing_kernel, smoothing_length) +boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density, + cavity.boundary.mass, + AdamiPressureExtrapolation(), + viscosity=viscosity, + smoothing_kernel, smoothing_length) boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass, AdamiPressureExtrapolation(), viscosity=viscosity, smoothing_kernel, smoothing_length) -# ========================================================================================== -# ==== Systems - -fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, - smoothing_length, - sound_speed, - viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure), - acceleration=(0.0, gravity)) - -boundary_system_tank = BoundarySPHSystem(tank.boundary, boundary_model_tank) +boundary_system_cavity = BoundarySPHSystem(cavity.boundary, boundary_model_cavity) boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=movement) # ========================================================================================== # ==== Simulation - -tspan = (0.0, 5.0) - -semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_lid, +bnd_thickness = boundary_layers * particle_spacing +semi = Semidiscretization(fluid_system, boundary_system_cavity, boundary_system_lid, neighborhood_search=GridNeighborhoodSearch, - periodic_box_min_corner=[0.0 - 4 * particle_spacing, -0.24], - periodic_box_max_corner=[1.0 + 4 * particle_spacing, - 1.0 + 4 * particle_spacing]) + periodic_box_min_corner=[-bnd_thickness, -bnd_thickness], + periodic_box_max_corner=cavity_size .+ + [bnd_thickness, bnd_thickness]) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02) -callbacks = CallbackSet(info_callback, saving_callback) +callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep()) # Use a Runge-Kutta method with automatic (error based) time step size control. # Enable threading of the RK method for better performance on multiple threads. @@ -108,7 +100,7 @@ callbacks = CallbackSet(info_callback, saving_callback) # Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces # become extremely large when fluid particles are very close to boundary particles, # and the time integration method interprets this as an instability. -sol = solve(ode, RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!), +sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration) dtmax=1e-2, # Limit stepsize to prevent crashing diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 6a49fa57b..9ab999e68 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -8,21 +8,28 @@ using TrixiParticles using OrdinaryDiffEq -acceleration_x = 2.5e-4 +# ========================================================================================== +# ==== Resolution +n_particles_x = 144 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 3 +spacing_ratio = 1 # ========================================================================================== -# ==== Fluid +# ==== Experiment Setup +tspan = (0.0, 5.0) -n_particles_x = 144 +acceleration_x = 2.5e-4 +# Boundary geometry and initial fluid particle positions cylinder_radius = 0.02 -box_length = 6cylinder_radius -box_height = 4cylinder_radius +tank_size = (6cylinder_radius, 4cylinder_radius) +fluid_size = tank_size +initial_velocity = (1.2e-4, 0.0) fluid_density = 1000.0 -nu = 0.1 / fluid_density - -velocity_ref = 1.2e-4 +nu = 0.1 / fluid_density # viscosity parameter # Adami uses c = 0.1*sqrt(acceleration_x*cylinder_radius) but from the original setup # from M. Ellero and N. A. Adams (https://doi.org/10.1002/nme.3088) is using c = 0.02 @@ -30,44 +37,49 @@ sound_speed = 0.02 pressure = sound_speed^2 * fluid_density -particle_spacing = box_length / n_particles_x -smoothing_length = 1.2 * particle_spacing -smoothing_kernel = SchoenbergQuarticSplineKernel{2}() +particle_spacing = tank_size[1] / n_particles_x -box = RectangularTank(particle_spacing, (box_length, box_height), (box_length, box_height), - fluid_density, n_layers=3, spacing_ratio=1.0, pressure=pressure, - faces=(false, false, true, true)) +box = RectangularTank(particle_spacing, fluid_size, tank_size, + fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, + pressure=pressure, faces=(false, false, true, true)) -cylinder = SphereShape(particle_spacing, cylinder_radius, (box_length / 2, box_height / 2), +cylinder = SphereShape(particle_spacing, cylinder_radius, tank_size ./ 2, fluid_density, sphere_type=RoundSphere()) fluid = setdiff(box.fluid, cylinder) boundary = union(cylinder, box.boundary) -boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass, - AdamiPressureExtrapolation(), - viscosity=ViscosityAdami(nu), - smoothing_kernel, smoothing_length) - +# ========================================================================================== +# ==== Fluid +smoothing_length = 1.2 * particle_spacing +smoothing_kernel = SchoenbergQuarticSplineKernel{2}() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, - sound_speed, viscosity=ViscosityAdami(nu), + sound_speed, viscosity=ViscosityAdami(; nu), transport_velocity=TransportVelocityAdami(pressure), acceleration=(acceleration_x, 0.0)) +# ========================================================================================== +# ==== Boundary +boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass, + AdamiPressureExtrapolation(), + viscosity=ViscosityAdami(; nu), + smoothing_kernel, smoothing_length) + boundary_system = BoundarySPHSystem(boundary, boundary_model) +# ========================================================================================== +# ==== Simulation semi = Semidiscretization(fluid_system, boundary_system, neighborhood_search=GridNeighborhoodSearch, - periodic_box_min_corner=[0.0, -box_length * 2], - periodic_box_max_corner=[box_length, box_length * 2]) + periodic_box_min_corner=[0.0, -tank_size[2]], + periodic_box_max_corner=[tank_size[1], 2 * tank_size[2]]) -tspan = (0.0, 5.0) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) -saving_callback = SolutionSavingCallback(dt=0.02) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") -callbacks = CallbackSet(info_callback, saving_callback) +callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep()) # Use a Runge-Kutta method with automatic (error based) time step size control. # Enable threading of the RK method for better performance on multiple threads. @@ -77,7 +89,7 @@ callbacks = CallbackSet(info_callback, saving_callback) # Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces # become extremely large when fluid particles are very close to boundary particles, # and the time integration method interprets this as an instability. -sol = solve(ode, RDPK3SpFSAL49((step_limiter!)=TrixiParticles.update_transport_velocity!), +sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) dtmax=1e-2, # Limit stepsize to prevent crashing diff --git a/examples/fluid/rectangular_tank_edac_2d.jl b/examples/fluid/rectangular_tank_edac_2d.jl index c0e23b176..de6a5ad59 100644 --- a/examples/fluid/rectangular_tank_edac_2d.jl +++ b/examples/fluid/rectangular_tank_edac_2d.jl @@ -19,10 +19,7 @@ initial_fluid_size = (2.0, 0.9) tank_size = (2.0, 1.0) fluid_density = 1000.0 -atmospheric_pressure = 100000.0 sound_speed = 10 * sqrt(gravity * initial_fluid_size[2]) -state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure, - background_pressure=atmospheric_pressure) tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, @@ -33,7 +30,6 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fl smoothing_length = 1.2 * fluid_particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() -fluid_density_calculator = ContinuityDensity() alpha = 0.02 viscosity = ViscosityAdami(nu=alpha * smoothing_length * sound_speed / 8) @@ -45,7 +41,6 @@ fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, smoothi # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - state_equation=state_equation, boundary_density_calculator, smoothing_kernel, smoothing_length) From 2c45e13ef3a3a43c376a6fe96219e15ab568cffb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 5 Jan 2024 14:39:19 +0100 Subject: [PATCH 20/62] time dependent initial velocity function --- examples/fluid/taylor_green_vortex_2d.jl | 2 -- src/schemes/fluid/fluid.jl | 3 +-- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 188346a48..a9eb92e68 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -28,8 +28,6 @@ b = -8pi^2 / reynolds_number p(pos) = -U^2 * exp(2 * b * 0) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 p(pos, t) = -U^2 * exp(2 * b * t) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 -v_x(pos) = -U * exp(b * 0) * cos(2pi * pos[1]) * sin(2pi * pos[2]) -v_y(pos) = U * exp(b * 0) * sin(2pi * pos[1]) * cos(2pi * pos[2]) v_x(pos, t) = -U * exp(b * t) * cos(2pi * pos[1]) * sin(2pi * pos[2]) v_y(pos, t) = U * exp(b * t) * sin(2pi * pos[1]) * cos(2pi * pos[2]) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 93762dc32..1328f739d 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -45,9 +45,8 @@ end @inline function initial_velocity(system, particle, init_velocity_function) position = initial_coords(system, particle) - v_init = SVector(ntuple(i -> init_velocity_function[i](position), Val(ndims(system)))) - return v_init + return SVector(ntuple(i -> init_velocity_function[i](position, 0), Val(ndims(system)))) end @inline viscosity_model(system::FluidSystem) = system.viscosity From d117b8ed30ecd4596359ecb4eb6f69ac7232177f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 5 Jan 2024 15:32:27 +0100 Subject: [PATCH 21/62] time dependent pressure function --- examples/fluid/taylor_green_vortex_2d.jl | 1 - src/schemes/fluid/entropically_damped_sph/system.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index a9eb92e68..37bdb4370 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -25,7 +25,6 @@ sound_speed = 10U b = -8pi^2 / reynolds_number -p(pos) = -U^2 * exp(2 * b * 0) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 p(pos, t) = -U^2 * exp(2 * b * t) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 v_x(pos, t) = -U * exp(b * t) * cos(2pi * pos[1]) * sin(2pi * pos[2]) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 40ba090dd..27f32501a 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -225,5 +225,5 @@ end @inline function initial_pressure(system, particle, initial_pressure_function) particle_position = initial_coords(system, particle) - return initial_pressure_function(particle_position) + return initial_pressure_function(particle_position, 0) end From 6684cc44aef7b6c5ec79d32b21e0c82a7871bdc2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 27 Jan 2024 17:45:54 +0100 Subject: [PATCH 22/62] remove velocity function --- examples/fluid/taylor_green_vortex_2d.jl | 26 +++++---- src/general/initial_condition.jl | 2 +- .../fluid/entropically_damped_sph/system.jl | 57 ++++++------------- src/schemes/fluid/fluid.jl | 17 +----- src/schemes/fluid/transport_velocity.jl | 7 ++- test/systems/edac_system.jl | 10 ++-- 6 files changed, 42 insertions(+), 77 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 37bdb4370..8288402bd 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -26,9 +26,18 @@ sound_speed = 10U b = -8pi^2 / reynolds_number p(pos, t) = -U^2 * exp(2 * b * t) * (cos(4pi * pos[1]) + cos(4pi * pos[2])) / 4 +p0(pos) = p(pos, 0.0) -v_x(pos, t) = -U * exp(b * t) * cos(2pi * pos[1]) * sin(2pi * pos[2]) -v_y(pos, t) = U * exp(b * t) * sin(2pi * pos[1]) * cos(2pi * pos[2]) +function velocity_function(pos, t) + x = pos[1] + y = pos[2] + + vel = U * exp(b * t) * [-cos(2pi * x) * sin(2pi * y), sin(2pi * x) * cos(2pi * y)] + + return SVector{2}(vel) +end + +v0(pos) = velocity_function(pos, 0.0) n_particles_xy = round(Int, box_length / particle_spacing) @@ -42,8 +51,7 @@ smoothing_length = 1.0 * particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0), - fluid_density, pressure=background_pressure, - init_velocity=(0.0, 0.0)) + density=fluid_density, pressure=p0, velocity=v0) # Add small random displacement to the particles to avoid stagnant streamlines. #seed!(42); @@ -51,8 +59,7 @@ fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0. # size(fluid.coordinates)) fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, - sound_speed, initial_pressure_function=p, - initial_velocity_function=(v_x, v_y), + sound_speed, transport_velocity=TransportVelocityAdami(background_pressure), viscosity=ViscosityAdami(; nu)) @@ -69,8 +76,6 @@ ode = semidiscretize(semi, tspan) dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu)) function compute_L1v_error(v, u, t, system) - (; initial_velocity_function) = system - v_analytical_avg = 0.0 L1v = 0.0 @@ -78,8 +83,7 @@ function compute_L1v_error(v, u, t, system) position = TrixiParticles.current_coords(u, system, particle) v_mag = norm(TrixiParticles.current_velocity(v, system, particle)) - v_analytical = norm(SVector(ntuple(i -> initial_velocity_function[i](position, t), - Val(ndims(system))))) + v_analytical = norm(velocity_function(position, t)) v_analytical_avg += abs(v_analytical) L1v += abs(v_mag - v_analytical) @@ -100,7 +104,7 @@ function compute_L1p_error(v, u, t, system) position = TrixiParticles.current_coords(u, system, particle) # compute pressure error - p_analytical = system.initial_pressure_function(position, t) + p_analytical = p(position, t) p_max_exact = max(p_max_exact, abs(p_analytical)) # p_computed - p_average diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 764b82beb..cc38cd41d 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -111,7 +111,7 @@ struct InitialCondition{ELTYPE} "for $NDIMS-dimensional `coordinates`")) end velocities_svector = velocity_fun.(coordinates_svector) - velocities = reinterpret(reshape, ELTYPE, velocities_svector) + velocities = stack(velocities_svector) end if density isa AbstractVector diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 27f32501a..b5538068b 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -2,8 +2,6 @@ EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; alpha=0.5, viscosity=NoViscosity(), - initial_pressure_function=nothing, - initial_velocity_function=nothing, transport_velocity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS)) @@ -35,28 +33,24 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). In: Computers and Fluids 179 (2019), pages 579-594. [doi: 10.1016/j.compfluid.2018.11.023](https://doi.org/10.1016/j.compfluid.2018.11.023) """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, VF, TV, C} <: +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, TV, C} <: FluidSystem{NDIMS} - initial_condition :: InitialCondition{ELTYPE} - mass :: Array{ELTYPE, 1} # [particle] - density :: Array{ELTYPE, 1} # [particle] - density_calculator :: DC - smoothing_kernel :: K - smoothing_length :: ELTYPE - sound_speed :: ELTYPE - viscosity :: V - nu_edac :: ELTYPE - initial_pressure_function :: PF - initial_velocity_function :: VF - acceleration :: SVector{NDIMS, ELTYPE} - transport_velocity :: TV - cache :: C + initial_condition :: InitialCondition{ELTYPE} + mass :: Array{ELTYPE, 1} # [particle] + density :: Array{ELTYPE, 1} # [particle] + density_calculator :: DC + smoothing_kernel :: K + smoothing_length :: ELTYPE + sound_speed :: ELTYPE + viscosity :: V + nu_edac :: ELTYPE + acceleration :: SVector{NDIMS, ELTYPE} + transport_velocity :: TV + cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; alpha=0.5, viscosity=NoViscosity(), - initial_pressure_function=nothing, - initial_velocity_function=nothing, transport_velocity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel))) @@ -83,12 +77,10 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, PF, VF, TV, cache = create_cache_edac(initial_condition, transport_velocity) new{NDIMS, ELTYPE, typeof(density_calculator), typeof(smoothing_kernel), - typeof(viscosity), typeof(initial_pressure_function), - typeof(initial_velocity_function), typeof(transport_velocity), + typeof(viscosity), typeof(transport_velocity), typeof(cache)}(initial_condition, mass, density, density_calculator, smoothing_kernel, smoothing_length, sound_speed, viscosity, - nu_edac, initial_pressure_function, initial_velocity_function, - acceleration_, transport_velocity, cache) + nu_edac, acceleration_, transport_velocity, cache) end end @@ -113,8 +105,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "initial pressure function", system.initial_pressure_function) - summary_line(io, "initial velocity function", system.initial_velocity_function) summary_line(io, "acceleration", system.acceleration) summary_footer(io) end @@ -180,7 +170,7 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode system_coords = current_coordinates(u, system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - neighborhood_search = neighborhood_searches(system, neighbor_system, semi) + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) # Loop over all pairs of particles and neighbors within the kernel cutoff. for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, @@ -201,7 +191,7 @@ end function write_v0!(v0, density_calculator, system::EntropicallyDampedSPHSystem) for particle in eachparticle(system) - v0[end, particle] = initial_pressure(system, particle) + v0[end, particle] = system.initial_condition.pressure[particle] end return v0 @@ -214,16 +204,3 @@ function restart_with!(system::EntropicallyDampedSPHSystem, v, u) system.initial_condition.pressure[particle] = v[end, particle] end end - -@inline function initial_pressure(system, particle) - initial_pressure(system, particle, system.initial_pressure_function) -end - -@inline function initial_pressure(system, particle, ::Nothing) - return system.initial_condition.pressure[particle] -end - -@inline function initial_pressure(system, particle, initial_pressure_function) - particle_position = initial_coords(system, particle) - return initial_pressure_function(particle_position, 0) -end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index c2549613f..2cd159350 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -20,9 +20,8 @@ end function write_v0!(v0, system::FluidSystem, ::Nothing) for particle in eachparticle(system) # Write particle velocities - v_init = initial_velocity(system, particle) for dim in 1:ndims(system) - v0[dim, particle] = v_init[dim] + v0[dim, particle] = system.initial_condition.velocity[dim, particle] end end @@ -31,20 +30,6 @@ function write_v0!(v0, system::FluidSystem, ::Nothing) return v0 end -@inline function initial_velocity(system, particle) - initial_velocity(system, particle, system.initial_velocity_function) -end - -@inline function initial_velocity(system, particle, ::Nothing) - return extract_svector(system.initial_condition.velocity, system, particle) -end - -@inline function initial_velocity(system, particle, init_velocity_function) - position = initial_coords(system, particle) - - return SVector(ntuple(i -> init_velocity_function[i](position, 0), Val(ndims(system)))) -end - @inline viscosity_model(system::FluidSystem) = system.viscosity include("pressure_acceleration.jl") diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 0a291d624..6828106ff 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -30,12 +30,13 @@ end end function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami) + (; initial_condition) = system + for particle in eachparticle(system) # Write particle velocities - v_init = initial_velocity(system, particle) for dim in 1:ndims(system) - v0[dim, particle] = v_init[dim] - v0[ndims(system) + dim, particle] = v_init[dim] + v0[dim, particle] = initial_condition.velocity[dim, particle] + v0[ndims(system) + dim, particle] = initial_condition.velocity[dim, particle] end end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 9fc99fbfe..c135fa8ad 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -33,8 +33,6 @@ @test system.smoothing_kernel == smoothing_kernel @test system.smoothing_length == smoothing_length @test system.viscosity isa TrixiParticles.NoViscosity - @test system.initial_pressure_function isa Nothing - @test system.initial_velocity_function isa Nothing @test system.transport_velocity isa Nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -90,8 +88,6 @@ @test system.smoothing_kernel == smoothing_kernel @test system.smoothing_length == smoothing_length @test system.viscosity isa TrixiParticles.NoViscosity - @test system.initial_pressure_function isa Nothing - @test system.initial_velocity_function isa Nothing @test system.transport_velocity isa Nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -196,9 +192,11 @@ @test v0 == vcat(velocity, pressure') + initial_condition = InitialCondition(; coordinates, velocity, mass, density, + pressure=pressure_function) + system = EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, - smoothing_length, sound_speed, - initial_pressure_function=pressure_function) + smoothing_length, sound_speed) v0 = zeros(TrixiParticles.v_nvariables(system), TrixiParticles.n_moving_particles(system)) From 2cd470ac4895360e7ae7725f89ac4b1fe895ba32 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 27 Jan 2024 18:55:29 +0100 Subject: [PATCH 23/62] multi dimensional functions --- src/general/initial_condition.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 764b82beb..79d544e10 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -111,7 +111,7 @@ struct InitialCondition{ELTYPE} "for $NDIMS-dimensional `coordinates`")) end velocities_svector = velocity_fun.(coordinates_svector) - velocities = reinterpret(reshape, ELTYPE, velocities_svector) + velocities = stack(velocities_svector) end if density isa AbstractVector @@ -135,7 +135,10 @@ struct InitialCondition{ELTYPE} pressures = pressure else pressure_fun = wrap_function(pressure, Val(NDIMS)) - pressures = pressure_fun.(coordinates_svector) + pressures = stack((pressure_fun.(coordinates_svector))) + if !isa(pressures, Vector) + pressures = pressures[findfirst(p -> abs(p) > eps(), pressures)[1], :] + end end if mass isa AbstractVector From 3207c9a2b74c34938c115b8bb93743dcfb47b8d3 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 27 Jan 2024 19:12:07 +0100 Subject: [PATCH 24/62] apply formatter --- src/general/initial_condition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 79d544e10..e921ab4dc 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -136,7 +136,7 @@ struct InitialCondition{ELTYPE} else pressure_fun = wrap_function(pressure, Val(NDIMS)) pressures = stack((pressure_fun.(coordinates_svector))) - if !isa(pressures, Vector) + if !isa(pressures, Vector) pressures = pressures[findfirst(p -> abs(p) > eps(), pressures)[1], :] end end From 5cb9586e6f7a071dad28e129a6d9246d3a8a7568 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sun, 28 Jan 2024 20:02:50 +0100 Subject: [PATCH 25/62] fix for open boundaries --- src/schemes/fluid/fluid.jl | 22 ++++++++++++++++++++++ src/schemes/fluid/transport_velocity.jl | 19 ------------------- 2 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 2cd159350..ebbb6b53e 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -37,3 +37,25 @@ include("viscosity.jl") include("transport_velocity.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") + +@inline function add_velocity!(du, v, particle, + system::Union{EntropicallyDampedSPHSystem, + WeaklyCompressibleSPHSystem}) + add_velocity!(du, v, particle, system, system.transport_velocity) +end + +@inline function momentum_convection(system, neighbor_system, + 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::Union{EntropicallyDampedSPHSystem, + WeaklyCompressibleSPHSystem}, + 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) +end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 6828106ff..01f397a7a 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -45,10 +45,6 @@ function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami) return v0 end -@inline function add_velocity!(du, v, particle, system::FluidSystem) - add_velocity!(du, v, particle, system, system.transport_velocity) -end - # Add momentum velocity. @inline function add_velocity!(du, v, particle, system, ::Nothing) for i in 1:ndims(system) @@ -71,21 +67,6 @@ end return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) end -@inline function momentum_convection(system, neighbor_system, - 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::FluidSystem, - neighbor_system::FluidSystem, - 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) -end - @inline function momentum_convection(system, neighbor_system, ::Nothing, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) From 38902d9df68cd671b56128e0b88b8565ef520a0f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 22 Feb 2024 11:05:59 +0100 Subject: [PATCH 26/62] add `UpdateCallback` --- src/TrixiParticles.jl | 2 +- src/callbacks/callbacks.jl | 1 + src/callbacks/update.jl | 139 ++++++++++++++++++++++++++++++++++++ test/callbacks/callbacks.jl | 1 + test/callbacks/update.jl | 65 +++++++++++++++++ 5 files changed, 207 insertions(+), 1 deletion(-) create mode 100644 src/callbacks/update.jl create mode 100644 test/callbacks/update.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 3aa241248..d8fb0fe31 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -44,7 +44,7 @@ export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, BoundarySPHSystem export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback + PostprocessCallback, StepsizeCallback, UpdateCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 7cc580038..62b6c106c 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -31,3 +31,4 @@ include("solution_saving.jl") include("density_reinit.jl") include("post_process.jl") include("stepsize.jl") +include("update.jl") diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl new file mode 100644 index 000000000..6bb1b7736 --- /dev/null +++ b/src/callbacks/update.jl @@ -0,0 +1,139 @@ +struct UpdateCallback{I} + interval :: I + update :: Bool +end + +""" + UpdateCallback(; update=true, interval::Integer, dt=0.0) + +Callback to update quantities either at the end of every `interval` integration step or at +regular intervals at `dt` in terms of integration time. + +# Keywords +- `update`: Callback is only applied when `true` (default) +- `interval`: Update quantities at the end of every `interval` time steps (default `inverval=1`) +- `dt`: Update quantities in regular intervals of `dt` in terms of integration time +""" +function UpdateCallback(; update=true, 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 update && interval == -1 + interval = 1 + end + + update_callback! = UpdateCallback(interval, update) + + if dt > 0 && update + # Add a `tstop` every `dt`, and save the final solution. + return PeriodicCallback(update_callback!, dt, + initialize=initial_update!, + save_positions=(false, false)) + else + # The first one is the condition, the second the affect! + return DiscreteCallback(update_callback!, update_callback!, + initialize=initial_update!, + save_positions=(false, false)) + end +end + +# initialize +function initial_update!(cb, u, t, integrator) + # The `UpdateCallback` is either `cb.affect!` (with `DiscreteCallback`) + # or `cb.affect!.affect!` (with `PeriodicCallback`). + # Let recursive dispatch handle this. + + initial_update!(cb.affect!, u, t, integrator) +end + +initial_update!(cb::UpdateCallback, u, t, integrator) = cb.update && cb(integrator) + +# condition +function (update_callback!::UpdateCallback)(u, t, integrator) + (; interval, update) = update_callback! + + # With error-based step size control, some steps can be rejected. Thus, + # `integrator.iter >= integrator.stats.naccept` + # (total #steps) (#accepted steps) + # We need to check the number of accepted steps since callbacks are not + # activated after a rejected step. + return update && (integrator.stats.naccept % interval == 0) +end + +# affect +function (update_callback!::UpdateCallback)(integrator) + t = integrator.t + semi = integrator.p + v_ode, u_ode = integrator.u.x + + # Update quantities that are stored in the systems. These quantities (e.g. pressure) + # still have the values from the last stage of the previous step if not updated here. + update_systems_and_nhs(v_ode, u_ode, semi, t) + + # Other updates might be added here later (e.g. Transport Velocity Formulation). + # @trixi_timeit timer() "update open boundary" foreach_system(semi) do system + # update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) + # end + # + # @trixi_timeit timer() "update TVF" foreach_system(semi) do system + # update_transport_velocity_eachstep!(system, v_ode, u_ode, semi, t) + # end + + + # Tell OrdinaryDiffEq that u has been modified + u_modified!(integrator, true) + + return integrator +end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:UpdateCallback}) + @nospecialize cb # reduce precompilation time + print(io, "UpdateCallback(interval=", (cb.affect!.update ? cb.affect!.interval : "-"), + ")") +end + +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:UpdateCallback}}) + @nospecialize cb # reduce precompilation time + print(io, "UpdateCallback(dt=", cb.affect!.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:UpdateCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect! + setup = [ + "interval" => update_cb.update ? update_cb.interval : "-", + "update" => update_cb.update ? "yes" : "no", + ] + summary_box(io, "UpdateCallback", setup) + end +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:UpdateCallback}}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect!.affect! + setup = [ + "dt" => update_cb.interval, + "update" => update_cb.update ? "yes" : "no", + ] + summary_box(io, "UpdateCallback", setup) + end +end diff --git a/test/callbacks/callbacks.jl b/test/callbacks/callbacks.jl index 3076b9274..3acaa93b5 100644 --- a/test/callbacks/callbacks.jl +++ b/test/callbacks/callbacks.jl @@ -2,4 +2,5 @@ include("info.jl") include("stepsize.jl") include("postprocess.jl") + include("update.jl") end diff --git a/test/callbacks/update.jl b/test/callbacks/update.jl new file mode 100644 index 000000000..5515b94a6 --- /dev/null +++ b/test/callbacks/update.jl @@ -0,0 +1,65 @@ +@testset verbose=true "UpdateCallback" begin + @testset verbose=true "show" begin + # Default + callback0 = UpdateCallback() + + show_compact = "UpdateCallback(interval=1)" + @test repr(callback0) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ UpdateCallback │ + │ ══════════════ │ + │ interval: ……………………………………………………… 1 │ + │ update: …………………………………………………………… yes │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback0) == show_box + + callback1 = UpdateCallback(interval=11) + + show_compact = "UpdateCallback(interval=11)" + @test repr(callback1) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ UpdateCallback │ + │ ══════════════ │ + │ interval: ……………………………………………………… 11 │ + │ update: …………………………………………………………… yes │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback1) == show_box + + callback2 = UpdateCallback(dt=1.2) + + show_compact = "UpdateCallback(dt=1.2)" + @test repr(callback2) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ UpdateCallback │ + │ ══════════════ │ + │ dt: ……………………………………………………………………… 1.2 │ + │ update: …………………………………………………………… yes │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback2) == show_box + + callback3 = UpdateCallback(update=false) + + show_compact = "UpdateCallback(interval=-)" + @test repr(callback3) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ UpdateCallback │ + │ ══════════════ │ + │ interval: ……………………………………………………… - │ + │ update: …………………………………………………………… no │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback3) == show_box + end + + @testset "Illegal Input" begin + error_str = "Setting both interval and dt is not supported!" + @test_throws ArgumentError(error_str) UpdateCallback(dt=0.1, interval=1) + end +end From cc3477cd3968505ddd016af48b8c3cd2c66bcf47 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 22 Feb 2024 11:07:45 +0100 Subject: [PATCH 27/62] fix typo --- src/callbacks/update.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 6bb1b7736..875fe63e2 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -11,7 +11,7 @@ regular intervals at `dt` in terms of integration time. # Keywords - `update`: Callback is only applied when `true` (default) -- `interval`: Update quantities at the end of every `interval` time steps (default `inverval=1`) +- `interval`: Update quantities at the end of every `interval` time steps (default `interval=1`) - `dt`: Update quantities in regular intervals of `dt` in terms of integration time """ function UpdateCallback(; update=true, interval::Integer=-1, dt=0.0) From 65823a4fb59240c00b8e26def7d6f446112fce7d Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 22 Feb 2024 11:10:43 +0100 Subject: [PATCH 28/62] prepare for merge `update-callback` --- examples/fluid/lid_driven_cavity_2d.jl | 2 +- .../fluid/periodic_array_of_cylinders_2d.jl | 2 +- examples/fluid/taylor_green_vortex_2d.jl | 2 +- src/callbacks/callbacks.jl | 1 - src/callbacks/update.jl | 44 ------------------- 5 files changed, 3 insertions(+), 48 deletions(-) delete mode 100644 src/callbacks/update.jl diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index f50480f57..28505ec18 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -90,7 +90,7 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02) -callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep()) +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) # Use a Runge-Kutta method with automatic (error based) time step size control. # Enable threading of the RK method for better performance on multiple threads. diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 9ab999e68..c2523b817 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -79,7 +79,7 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") -callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep()) +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) # Use a Runge-Kutta method with automatic (error based) time step size control. # Enable threading of the RK method for better performance on multiple threads. diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 49f4d20be..f002ca20d 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -142,7 +142,7 @@ saving_callback = SolutionSavingCallback(dt=0.02, p_avg=average_pressure, t=time_vector) -callbacks = CallbackSet(info_callback, saving_callback, UpdateEachTimeStep()) +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) # Use a Runge-Kutta method with automatic (error based) time step size control. # Enable threading of the RK method for better performance on multiple threads. diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 39ad4b564..7cc580038 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -29,6 +29,5 @@ end include("info.jl") include("solution_saving.jl") include("density_reinit.jl") -include("update.jl") include("post_process.jl") include("stepsize.jl") diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl deleted file mode 100644 index 9bfcb067c..000000000 --- a/src/callbacks/update.jl +++ /dev/null @@ -1,44 +0,0 @@ -# TODO: naming -struct UpdateEachDt end - -function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:UpdateEachDt}) - @nospecialize cb # reduce precompilation time - print(io, "UpdateEachTimeStep()") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, <:UpdateEachDt}) - @nospecialize cb # reduce precompilation time - if get(io, :compact, false) - show(io, cb) - else - summary_box(io, "UpdateEachTimeStep") - end -end - -function UpdateEachTimeStep() - update_each_dt = UpdateEachDt() - return DiscreteCallback(update_each_dt, update_each_dt, initialize=initial_update!, - save_positions=(false, false)) -end - -# condition -(update_each_dt::UpdateEachDt)(u, t, integrator) = true - -# affect -function (update_each_dt::UpdateEachDt)(integrator) - semi = integrator.p - v_ode, u_ode = integrator.u.x - - foreach_system(semi) do system - update_transport_velocity!(system, v_ode, semi) - end - - # Tell OrdinaryDiffEq that u has been modified - u_modified!(integrator, true) - - return integrator -end - -# initialize -initial_update!(cb, u, t, integrator) = cb.affect!(integrator) From 04fdc54b1498ebde8caf2ecdb14ac5d9081b50f3 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 5 Mar 2024 19:11:32 +0100 Subject: [PATCH 29/62] fix bug --- .../fluid/weakly_compressible_sph/system.jl | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 8e98e3a13..85bfc2bb5 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -53,6 +53,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF + transport_velocity :: Nothing # TODO source_terms :: ST cache :: C @@ -96,17 +97,13 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, create_cache_wcsph(correction, initial_condition.density, NDIMS, n_particles)..., cache...) - return new{NDIMS, ELTYPE, typeof(density_calculator), - typeof(state_equation), typeof(smoothing_kernel), - typeof(viscosity), typeof(density_diffusion), - typeof(correction), typeof(pressure_acceleration), - typeof(source_terms), typeof(cache)}(initial_condition, mass, pressure, - density_calculator, state_equation, - smoothing_kernel, smoothing_length, - acceleration_, viscosity, - density_diffusion, correction, - pressure_acceleration, - source_terms, cache) + return new{NDIMS, ELTYPE, typeof(density_calculator), typeof(state_equation), + typeof(smoothing_kernel), typeof(viscosity), typeof(density_diffusion), + typeof(correction), typeof(pressure_acceleration), typeof(source_terms), + typeof(cache)}(initial_condition, mass, pressure, density_calculator, + state_equation, smoothing_kernel, smoothing_length, + acceleration_, viscosity, density_diffusion, correction, + pressure_acceleration, nothing, source_terms, cache) end end From 2ae70a3f908719673ebd598f7e26277230c17928 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 5 Mar 2024 19:27:28 +0100 Subject: [PATCH 30/62] fix example --- examples/fluid/lid_driven_cavity_2d.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index 28505ec18..db39a6c2c 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -39,7 +39,7 @@ lid_position = 0.0 - particle_spacing * boundary_layers lid_length = cavity.n_particles_per_dimension[1] + 2boundary_layers lid = RectangularShape(particle_spacing, (lid_length, 3), - (lid_position, cavity_size[2]), fluid_density) + (lid_position, cavity_size[2]), density=fluid_density) # ========================================================================================== # ==== Fluid @@ -54,12 +54,11 @@ fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel, smoot # ========================================================================================== # ==== Boundary -f_y(t) = 0.0 -f_x(t) = velocity_lid * t +movement_function(t) = SVector(velocity_lid * t, 0.0) is_moving(t) = true -movement = BoundaryMovement((f_x, f_y), is_moving) +movement = BoundaryMovement(movement_function, is_moving) boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density, cavity.boundary.mass, From eee3d0710d3431d13a3ac3b1c5e2fa885ea6e385 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 07:57:21 +0100 Subject: [PATCH 31/62] fix tests --- examples/fluid/periodic_array_of_cylinders_2d.jl | 2 +- test/examples/examples.jl | 9 ++++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index c2523b817..88e915ef3 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -20,7 +20,7 @@ spacing_ratio = 1 # ==== Experiment Setup tspan = (0.0, 5.0) -acceleration_x = 2.5e-4 +const acceleration_x = 2.5e-4 # Boundary geometry and initial fluid particle positions cylinder_radius = 0.02 diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 0e38faf50..429a03b56 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -121,27 +121,30 @@ end @trixi_testset "fluid/lid_driven_cavity_2d.jl" begin - @test_nowarn trixi_include(@__MODULE__, + @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "lid_driven_cavity_2d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 end @trixi_testset "fluid/taylor_green_vortex_2d.jl" begin - @test_nowarn trixi_include(@__MODULE__, + @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "taylor_green_vortex_2d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 end @trixi_testset "fluid/periodic_array_of_cylinders_2d.jl" begin - @test_nowarn trixi_include(@__MODULE__, + @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "periodic_array_of_cylinders_2d.jl"), tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 end include("dam_break_2d_corrections.jl") From cf288218d545fb1e56aa2dff0dca2d16202e1ae2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 09:32:31 +0100 Subject: [PATCH 32/62] fix update bug --- src/callbacks/update.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 5bbf6fe9d..1d61170da 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -81,9 +81,10 @@ function (update_callback!::UpdateCallback)(integrator) # update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) # end # - # @trixi_timeit timer() "update TVF" foreach_system(semi) do system - # update_transport_velocity_eachstep!(system, v_ode, u_ode, semi, t) - # end + + @trixi_timeit timer() "update TVF" foreach_system(semi) do system + update_transport_velocity!(system, v_ode, semi) + end # Tell OrdinaryDiffEq that u has been modified u_modified!(integrator, true) From 4d3917425ab2b8ef7942e37900c17e24790cd04b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 09:35:17 +0100 Subject: [PATCH 33/62] apply formatter --- test/examples/examples.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 429a03b56..83c283c36 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -122,27 +122,27 @@ @trixi_testset "fluid/lid_driven_cavity_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "lid_driven_cavity_2d.jl"), - tspan=(0.0, 0.1)) + joinpath(examples_dir(), "fluid", + "lid_driven_cavity_2d.jl"), + tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @trixi_testset "fluid/taylor_green_vortex_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "taylor_green_vortex_2d.jl"), - tspan=(0.0, 0.1)) + joinpath(examples_dir(), "fluid", + "taylor_green_vortex_2d.jl"), + tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @trixi_testset "fluid/periodic_array_of_cylinders_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "periodic_array_of_cylinders_2d.jl"), - tspan=(0.0, 0.1)) + joinpath(examples_dir(), "fluid", + "periodic_array_of_cylinders_2d.jl"), + tspan=(0.0, 0.1)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end From 9c81ecfa1fcff98fd8bacb4dd128859b9ba24902 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 10:28:13 +0100 Subject: [PATCH 34/62] fix test --- examples/fluid/lid_driven_cavity_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index db39a6c2c..bddcf64f8 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -24,7 +24,7 @@ cavity_size = (1.0, 1.0) fluid_density = 1.0 -velocity_lid = 1.0 +const velocity_lid = 1.0 sound_speed = 10 * velocity_lid pressure = sound_speed^2 * fluid_density From 09870658da4309d4c480a51a476d112dde417bf1 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 15:45:11 +0100 Subject: [PATCH 35/62] minor changes --- examples/fluid/taylor_green_vortex_2d.jl | 12 ++++++------ src/schemes/fluid/entropically_damped_sph/system.jl | 2 ++ src/visualization/write2vtk.jl | 3 +++ test/systems/edac_system.jl | 1 + 4 files changed, 12 insertions(+), 6 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index f002ca20d..c5f0a23fb 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -127,19 +127,19 @@ function compute_L1p_error(v, u, t, system) return L1p /= p_max_exact end -function average_pressure(v, u, t, system) - return [TrixiParticles.average_pressure(system, i) for i in axes(v, 2)] -end +# The pressure plotted in the paper is the difference of the local pressure minus +# the average of the pressure of all particles. +function diff_p_loc_p_avg(v, u, t, system) + p_avg_tot = avg_pressure(v, u, t, system) -function time_vector(v, u, t, system) - return t + return v[end, :] .- p_avg_tot end info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02, L1v=compute_L1v_error, L1p=compute_L1p_error, - p_avg=average_pressure, + diff_p_loc_p_avg=diff_p_loc_p_avg, t=time_vector) callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index bda26c2b1..90770389f 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -124,6 +124,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) + summary_line(io, "tansport velocity formulation", + system.transport_velocity |> typeof |> nameof) summary_line(io, "acceleration", system.acceleration) summary_footer(io) end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index bb7fc101a..b2b211859 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -202,6 +202,9 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) else vtk["solver"] = "EDAC" vtk["sound_speed"] = system.sound_speed + vtk["background_pressure_TVF"] = system.transport_velocity isa Nothing ? + "-" : + system.transport_velocity.background_pressure end end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 479a99e18..da4409fde 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -139,6 +139,7 @@ │ viscosity: …………………………………………………… Nothing │ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ + │ tansport velocity formulation: Nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box From 1ce8e09e263ae741215cee98f4f1edcd519fda49 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 16:47:53 +0100 Subject: [PATCH 36/62] add tests --- test/schemes/fluid/fluid.jl | 1 + test/schemes/fluid/transport_velocity.jl | 35 ++++++++++++++++++++++++ test/systems/edac_system.jl | 24 ++++++++++++++++ 3 files changed, 60 insertions(+) create mode 100644 test/schemes/fluid/transport_velocity.jl diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index bf7c3720c..937bfde9a 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -1,3 +1,4 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") +include("transport_velocity.jl") diff --git a/test/schemes/fluid/transport_velocity.jl b/test/schemes/fluid/transport_velocity.jl new file mode 100644 index 000000000..f53cd813a --- /dev/null +++ b/test/schemes/fluid/transport_velocity.jl @@ -0,0 +1,35 @@ +@trixi_testset "Transport Velocity Formulation" begin + particle_spacing = 0.1 + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 1.2particle_spacing + + fluid = rectangular_patch(particle_spacing, (3, 3), seed=1) + + v0_tvf = zeros(5, nparticles(fluid)) + + system_tvf = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, + transport_velocity=TransportVelocityAdami(0.0), + smoothing_length, 0.0) + system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, 0.0) + + @testset "Number of Variables" begin + @test TrixiParticles.v_nvariables(system_tvf) == 5 + @test TrixiParticles.v_nvariables(system) == 3 + end + + @testset "write_v0!" begin + TrixiParticles.write_v0!(v0_tvf, system_tvf) + + @test vcat(fluid.velocity, fluid.velocity, fluid.pressure') ≈ v0_tvf + end + + @testset "Update" begin + semi = Semidiscretization(system_tvf) + fill!(v0_tvf, 1.5) + v0_tvf[1:2, :] .= 2.5 + + TrixiParticles.update_transport_velocity!(system_tvf, vec(v0_tvf), semi) + + @test fill(2.5, (4, nparticles(system_tvf))) ≈ v0_tvf[1:4, :] + end +end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index da4409fde..2b34ff2d9 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -204,4 +204,28 @@ @test v0 == vcat(velocity, [0.8, 1.0]') end + + @trixi_testset "Average Pressure" begin + particle_spacing = 0.1 + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 1.6particle_spacing + + fluid = rectangular_patch(particle_spacing, (3, 3), seed=1) + + system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, + transport_velocity=TransportVelocityAdami(0.0), + smoothing_length, 0.0) + semi = Semidiscretization(system) + + u_ode = vec(fluid.coordinates) + v_ode = vec(vcat(fluid.velocity, fluid.velocity, fluid.pressure')) + + TrixiParticles.update_average_pressure!(system, system.transport_velocity, v_ode, + u_ode, semi) + + @test all(i -> system.cache.neighbor_counter[i] == nparticles(system), + nparticles(system)) + @test all(i -> system.cache.pressure_average[i] ≈ -50.968532955185964, + nparticles(system)) + end end From 8315ac93b3f06fabad85980f2e0209e6d2ca4881 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 21:19:41 +0100 Subject: [PATCH 37/62] add docs --- docs/src/systems/entropically_damped_sph.md | 60 +++++++++++++++++++ .../fluid/entropically_damped_sph/rhs.jl | 6 +- .../fluid/entropically_damped_sph/system.jl | 4 ++ src/schemes/fluid/transport_velocity.jl | 17 +++++- 4 files changed, 85 insertions(+), 2 deletions(-) diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 499c5aafe..5a0257b00 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -54,3 +54,63 @@ Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")] - Jonathan R. Clausen. "Entropically damped form of artificial compressibility for explicit simulation of incompressible flow". In: American Physical Society 87 (2013), page 13309. [doi: 10.1103/PhysRevE.87.013309](http://dx.doi.org/10.1103/PhysRevE.87.013309) + +## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) +Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow. +Adami et al (2013) modified the advection velocity and added an extra term in the momentum equation to avoid this problems. +The authors introduced the so called Transport Velocity Formulation (TVF) for WCSPH. Ramachandran (2019) et al applied the TVF +also for the [EDAC](@ref edac) scheme. + +The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by +```math +\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a +``` +and is obtained at every time-step ``\Delta t`` from +```math +\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right) +``` + +where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field. + +The discretized form of the last term is +```math + -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab} +``` +Note that although ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, +which means that there is a non-vanishing contribution only when particles are disordered. +That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions. +Suggested is a background pressure which is on the order of the reference pressure but can be chosen arbitrarily large when time-step criterion is adjusted. + +The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is +```math +\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A} +``` + where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)`` is a consequence of the modified + advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. + +The discretized form of the momentum equation for a particle ``a`` reads as +```math +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right] +``` +where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively and ``v_{ab}= v_a -v_b`` is the relative velocity. +Here, ``\tilde{p}_{ab}`` is the density-weighted pressure +```math +\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, +``` +with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. + +The convection tensor is calculated as ``\bm{A} = \rho v (\tilde{v}-v)^T``. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")] +``` + +### References +- S. Adami, X. Y. Hu, N. A. Adams. + "A transport-velocity formulation for smoothed particle hydrodynamics". + In: Journal of Computational Physics 241, (2013), pages 292--307. + [doi: 10.1016/j.jcp.2013.01.043](http://dx.doi.org/10.1016/j.jcp.2013.01.043) +- Prabhu Ramachandran. "Entropically damped artificial compressibility for SPH". + In: Computers and Fluids 179 (2019), pages 579--594. + [doi: 10.1016/j.compfluid.2018.11.023](https://doi.org/10.1016/j.compfluid.2018.11.023) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index b4fddf9a4..af7fa59a6 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -21,7 +21,11 @@ function interact!(dv, v_particle_system, u_particle_system, p_a = particle_pressure(v_particle_system, particle_system, particle) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) - p_avg = average_pressure(particle_system, particle) # Only used with EDAC-TVF + + # This technique is for a more robust `pressure_acceleration` but only with TVF. + # It results only in signifcant improvement for EDAC and not for WCSPH. + # See Ramachandran (2019) p. 582 + p_avg = average_pressure(particle_system, particle) m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 90770389f..21467093a 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -29,6 +29,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) +- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). Default is no TVF. - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure)` (which are the quantities of a single particle), returning a `Tuple` @@ -185,6 +186,9 @@ function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi) return system end +# This technique is for a more robust `pressure_acceleration` but only with TVF. +# It results only in signifcant improvement for EDAC and not for WCSPH. +# See Ramachandran (2019) p. 582 function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode, semi) (; cache) = system (; pressure_average, neighbor_counter) = cache diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 66908b15f..aa1332e03 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -1,6 +1,21 @@ +""" + TransportVelocityAdami(background_pressure::Real) + +Transport Velocity Formulation (TVF) to supress pairing and tensile instability. +See [TVF](@ref transport_velocity_formulation) for more details on the method. + +# Arguments +- `background_pressure`: Background pressure. Suggested is a background pressure which is + on the order of the reference pressure. + +!!! note + There is no need for an artificial viscosity to suppress tensile instability when using `TransportVelocityAdami`. + Thus, it is highly recommended to use [`ViscosityAdami`](@ref) as viscosity model, + since [`ArtificialViscosityMonaghan`](@ref) leads to bad results. +""" struct TransportVelocityAdami{ELTYPE} background_pressure::ELTYPE - function TransportVelocityAdami(background_pressure) + function TransportVelocityAdami(background_pressure::Real) new{typeof(background_pressure)}(background_pressure) end end From a8d21a8b0076d31187475a91200261d5e4cc7ea3 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 6 Mar 2024 21:22:15 +0100 Subject: [PATCH 38/62] fix tpos --- src/schemes/fluid/entropically_damped_sph/rhs.jl | 2 +- src/schemes/fluid/entropically_damped_sph/system.jl | 2 +- src/schemes/fluid/transport_velocity.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index af7fa59a6..4f395140c 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -23,7 +23,7 @@ function interact!(dv, v_particle_system, u_particle_system, p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) # This technique is for a more robust `pressure_acceleration` but only with TVF. - # It results only in signifcant improvement for EDAC and not for WCSPH. + # It results only in significant improvement for EDAC and not for WCSPH. # See Ramachandran (2019) p. 582 p_avg = average_pressure(particle_system, particle) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 21467093a..52625c53f 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -187,7 +187,7 @@ function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi) end # This technique is for a more robust `pressure_acceleration` but only with TVF. -# It results only in signifcant improvement for EDAC and not for WCSPH. +# It results only in significant improvement for EDAC and not for WCSPH. # See Ramachandran (2019) p. 582 function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode, semi) (; cache) = system diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index aa1332e03..3385c340d 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -1,7 +1,7 @@ """ TransportVelocityAdami(background_pressure::Real) -Transport Velocity Formulation (TVF) to supress pairing and tensile instability. +Transport Velocity Formulation (TVF) to suppress pairing and tensile instability. See [TVF](@ref transport_velocity_formulation) for more details on the method. # Arguments From 7b8b1d5a82c920a686ddab2816ec9f88f25829e4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 8 Mar 2024 13:35:12 +0100 Subject: [PATCH 39/62] add configuration check --- src/general/semidiscretization.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 222b0d906..d3b3e8172 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -670,3 +670,12 @@ function check_configuration(system::TotalLagrangianSPHSystem, systems) "`ContinuityDensity` is not yet supported for a `TotalLagrangianSPHSystem`")) end end + +function check_configuration(system::FluidSystem, systems) + (; viscosity) = system + + if viscosity isa ArtificialViscosityMonaghan && + system.transport_velocity isa TransportVelocityAdami + throw(ArgumentError("Please use `ViscosityAdami` when simulating with `TransportVelocityAdami`")) + end +end From 4891f175dbe750d99c39f275b37e609536b5d983 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 25 Jun 2024 21:25:17 +0200 Subject: [PATCH 40/62] add setter for tvf --- src/schemes/boundary/open_boundary/system.jl | 3 ++- src/schemes/fluid/transport_velocity.jl | 21 ++++++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 8c580cf44..254b031ac 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -435,7 +435,8 @@ end v_new[dim, particle_new] = v_old[dim, particle_old] end - # TODO: Only when using TVF: set tvf + # Only when using TVF + set_transport_velocity!(system_new, particle_new, particle_old, v_new, v_old) return system_new end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 9ae26ca87..db43bd48b 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -177,3 +177,24 @@ function update_final!(system::FluidSystem, transport_velocity, return system end + +# WARNING! +# These functions are intended to be used internally to set the transport velocity +# of newly activated particles in a callback. +# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` +# outside of callbacks. +function set_transport_velocity!(system::FluidSystem, particle, particle_old, v, v_old) + set_transport_velocity!(system, particle, particle_old, v, v_old, + system.transport_velocity) +end + +set_transport_velocity!(system, particle, particle_old, v, v_old) = system + +set_transport_velocity!(system, particle, particle_old, v, v_old, ::Nothing) = system + +function set_transport_velocity!(system, particle, particle_old, v, v_old, + ::TransportVelocityAdami) + for i in 1:ndims(system) + v[ndims(system) + i, particle] = v_old[i, particle_old] + end +end From 16e7ce57edabc27e47c1ad29fee3bbeb8045fba0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 3 Jul 2024 16:12:53 +0200 Subject: [PATCH 41/62] fix callback used check --- src/schemes/fluid/transport_velocity.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index db43bd48b..46e9ad08a 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -142,7 +142,7 @@ end reset_callback_flag!(system, ::Nothing) = system -function reset_callback_flag!(system::FluidSystem, transport_velocity) +function reset_callback_flag!(system::FluidSystem, ::TransportVelocityAdami) system.cache.update_callback_used[] = false return system @@ -165,11 +165,11 @@ function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t; end function update_final!(system::FluidSystem, ::Nothing, - v, u, v_ode, u_ode, semi, t) + v, u, v_ode, u_ode, semi, t; update_from_callback=false) return system end -function update_final!(system::FluidSystem, transport_velocity, +function update_final!(system::FluidSystem, ::TransportVelocityAdami, v, u, v_ode, u_ode, semi, t; update_from_callback=false) if !update_from_callback && !(system.cache.update_callback_used[]) throw(ArgumentError("`UpdateCallback` is required when using `TransportVelocityAdami`")) From 5ab14ee147b3e5215fa0e1c58362ed43917fe8b5 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 6 Aug 2024 10:23:32 +0200 Subject: [PATCH 42/62] adapt examples --- examples/fluid/lid_driven_cavity_2d.jl | 32 ++++++++----------- .../fluid/periodic_array_of_cylinders_2d.jl | 24 +++++--------- 2 files changed, 21 insertions(+), 35 deletions(-) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index bddcf64f8..fd9064080 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -7,13 +7,13 @@ using TrixiParticles using OrdinaryDiffEq + # ========================================================================================== # ==== Resolution particle_spacing = 0.02 -# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled boundary_layers = 3 -spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup @@ -24,15 +24,15 @@ cavity_size = (1.0, 1.0) fluid_density = 1.0 -const velocity_lid = 1.0 -sound_speed = 10 * velocity_lid +const VELOCITY_LID = 1.0 +sound_speed = 10 * VELOCITY_LID pressure = sound_speed^2 * fluid_density -viscosity = ViscosityAdami(; nu=velocity_lid / reynolds_number) +viscosity = ViscosityAdami(; nu=VELOCITY_LID / reynolds_number) cavity = RectangularTank(particle_spacing, cavity_size, cavity_size, fluid_density, - n_layers=boundary_layers, spacing_ratio=spacing_ratio, + n_layers=boundary_layers, faces=(true, true, true, false), pressure=pressure) lid_position = 0.0 - particle_spacing * boundary_layers @@ -48,13 +48,14 @@ smoothing_length = 1.0 * particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel, smoothing_length, + density_calculator=ContinuityDensity(), sound_speed, viscosity=viscosity, transport_velocity=TransportVelocityAdami(pressure)) # ========================================================================================== # ==== Boundary -movement_function(t) = SVector(velocity_lid * t, 0.0) +movement_function(t) = SVector(VELOCITY_LID * t, 0.0) is_moving(t) = true @@ -78,11 +79,11 @@ boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=moveme # ========================================================================================== # ==== Simulation bnd_thickness = boundary_layers * particle_spacing +periodic_box = PeriodicBox(min_corner=[-bnd_thickness, -bnd_thickness], + max_corner=cavity_size .+ [bnd_thickness, bnd_thickness]) + semi = Semidiscretization(fluid_system, boundary_system_cavity, boundary_system_lid, - neighborhood_search=GridNeighborhoodSearch, - periodic_box_min_corner=[-bnd_thickness, -bnd_thickness], - periodic_box_max_corner=cavity_size .+ - [bnd_thickness, bnd_thickness]) + neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box)) ode = semidiscretize(semi, tspan) @@ -91,14 +92,7 @@ info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02) callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) -# Use a Runge-Kutta method with automatic (error based) time step size control. -# Enable threading of the RK method for better performance on multiple threads. -# Limiting of the maximum stepsize is necessary to prevent crashing. -# When particles are approaching a wall in a uniform way, they can be advanced -# with large time steps. Close to the wall, the stepsize has to be reduced drastically. -# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces -# become extremely large when fluid particles are very close to boundary particles, -# and the time integration method interprets this as an instability. +# Use a Runge-Kutta method with automatic (error based) time step size control sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration) diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 88e915ef3..6b70f4724 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -12,9 +12,8 @@ using OrdinaryDiffEq # ==== Resolution n_particles_x = 144 -# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled boundary_layers = 3 -spacing_ratio = 1 # ========================================================================================== # ==== Experiment Setup @@ -31,8 +30,8 @@ initial_velocity = (1.2e-4, 0.0) fluid_density = 1000.0 nu = 0.1 / fluid_density # viscosity parameter -# Adami uses c = 0.1*sqrt(acceleration_x*cylinder_radius) but from the original setup -# from M. Ellero and N. A. Adams (https://doi.org/10.1002/nme.3088) is using c = 0.02 +# Adami uses `c = 0.1 * sqrt(acceleration_x * cylinder_radius)`` but the original setup +# from M. Ellero and N. A. Adams (https://doi.org/10.1002/nme.3088) uses `c = 0.02` sound_speed = 0.02 pressure = sound_speed^2 * fluid_density @@ -40,7 +39,7 @@ pressure = sound_speed^2 * fluid_density particle_spacing = tank_size[1] / n_particles_x box = RectangularTank(particle_spacing, fluid_size, tank_size, - fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, + fluid_density, n_layers=boundary_layers, pressure=pressure, faces=(false, false, true, true)) cylinder = SphereShape(particle_spacing, cylinder_radius, tank_size ./ 2, @@ -69,10 +68,10 @@ boundary_system = BoundarySPHSystem(boundary, boundary_model) # ========================================================================================== # ==== Simulation +periodic_box = PeriodicBox(min_corner=[0.0, -tank_size[2]], + max_corner=[tank_size[1], 2 * tank_size[2]]) semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch, - periodic_box_min_corner=[0.0, -tank_size[2]], - periodic_box_max_corner=[tank_size[1], 2 * tank_size[2]]) + neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box)) ode = semidiscretize(semi, tspan) @@ -81,14 +80,7 @@ saving_callback = SolutionSavingCallback(dt=0.02, prefix="") callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) -# Use a Runge-Kutta method with automatic (error based) time step size control. -# Enable threading of the RK method for better performance on multiple threads. -# Limiting of the maximum stepsize is necessary to prevent crashing. -# When particles are approaching a wall in a uniform way, they can be advanced -# with large time steps. Close to the wall, the stepsize has to be reduced drastically. -# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces -# become extremely large when fluid particles are very close to boundary particles, -# and the time integration method interprets this as an instability. +# Use a Runge-Kutta method with automatic (error based) time step size control sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) From c491b1a5f27e5d471fd426f59913e0c6fb1fae43 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 12 Aug 2024 11:10:13 +0200 Subject: [PATCH 43/62] minor changes --- docs/src/systems/entropically_damped_sph.md | 29 ++++++++++++------- .../fluid/periodic_array_of_cylinders_2d.jl | 2 +- examples/fluid/taylor_green_vortex_2d.jl | 8 ++--- src/TrixiParticles.jl | 2 -- src/callbacks/update.jl | 2 +- src/schemes/boundary/open_boundary/system.jl | 3 -- src/schemes/fluid/transport_velocity.jl | 23 +-------------- 7 files changed, 26 insertions(+), 43 deletions(-) diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 5a0257b00..5734c2214 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -57,49 +57,58 @@ Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")] ## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow. -Adami et al (2013) modified the advection velocity and added an extra term in the momentum equation to avoid this problems. -The authors introduced the so called Transport Velocity Formulation (TVF) for WCSPH. Ramachandran (2019) et al applied the TVF +To address these problems, Adami et al. (2013) modified the advection velocity and added an extra term to the momentum equation. +The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. Ramachandran et al. (2019) applied the TVF also for the [EDAC](@ref edac) scheme. The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by + ```math \frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a ``` + and is obtained at every time-step ``\Delta t`` from + ```math -\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right) +\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), ``` where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field. The discretized form of the last term is + ```math - -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab} + -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, ``` + +where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. Note that although ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, which means that there is a non-vanishing contribution only when particles are disordered. That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions. Suggested is a background pressure which is on the order of the reference pressure but can be chosen arbitrarily large when time-step criterion is adjusted. The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is + ```math -\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A} +\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, ``` - where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)`` is a consequence of the modified + + where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. The discretized form of the momentum equation for a particle ``a`` reads as + ```math -\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right] +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. ``` -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively and ``v_{ab}= v_a -v_b`` is the relative velocity. + Here, ``\tilde{p}_{ab}`` is the density-weighted pressure + ```math \tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, ``` -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. -The convection tensor is calculated as ``\bm{A} = \rho v (\tilde{v}-v)^T``. +with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and is given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. ```@autodocs Modules = [TrixiParticles] diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 6b70f4724..2734a8529 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -23,7 +23,7 @@ const acceleration_x = 2.5e-4 # Boundary geometry and initial fluid particle positions cylinder_radius = 0.02 -tank_size = (6cylinder_radius, 4cylinder_radius) +tank_size = (6 * cylinder_radius, 4 * cylinder_radius) fluid_size = tank_size initial_velocity = (1.2e-4, 0.0) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 10ac35f29..964245502 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -25,7 +25,7 @@ sound_speed = 10U b = -8pi^2 / reynolds_number -# Pressure function +# Taylor Green Vortex Pressure Function function pressure_function(pos, t) x = pos[1] y = pos[2] @@ -35,7 +35,7 @@ end initial_pressure_function(pos) = pressure_function(pos, 0.0) -# Velocity function +# Taylor Green Vortex Velocity Function function velocity_function(pos, t) x = pos[1] y = pos[2] @@ -89,8 +89,8 @@ function compute_L1v_error(v, u, t, system) for particle in TrixiParticles.eachparticle(system) position = TrixiParticles.current_coords(u, system, particle) - v_mag = norm(TrixiParticles.current_velocity(v, system, particle)) - v_analytical = norm(velocity_function(position, t)) + v_mag = TrixiParticles.norm(TrixiParticles.current_velocity(v, system, particle)) + v_analytical = TrixiParticles.norm(velocity_function(position, t)) v_analytical_avg += abs(v_analytical) L1v += abs(v_mag - v_analytical) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 7eb2ad693..27b080215 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -22,8 +22,6 @@ using RecipesBase: RecipesBase, @series using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem @reexport using StaticArrays: SVector -@reexport using Random: seed! -@reexport using LinearAlgebra: norm using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer! diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index ce9ec77b7..821d36e33 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -80,7 +80,7 @@ function (update_callback!::UpdateCallback)(integrator) # still have the values from the last stage of the previous step if not updated here. update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) - # Other updates might be added here later (e.g. Transport Velocity Formulation). + # Update open boundaries first, since particles might be activated or deactivated @trixi_timeit timer() "update open boundary" foreach_system(semi) do system update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 294ab1a7f..f8e5f9f16 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -294,9 +294,6 @@ end v_new[dim, particle_new] = v_old[dim, particle_old] end - # Only when using TVF - set_transport_velocity!(system_new, particle_new, particle_old, v_new, v_old) - return system_new end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 46e9ad08a..8080f4702 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -2,7 +2,7 @@ TransportVelocityAdami(background_pressure::Real) Transport Velocity Formulation (TVF) to suppress pairing and tensile instability. -See [TVF](@ref transport_velocity_formulation) for more details on the method. +See [TVF](@ref transport_velocity_formulation) for more details of the method. # Arguments - `background_pressure`: Background pressure. Suggested is a background pressure which is @@ -177,24 +177,3 @@ function update_final!(system::FluidSystem, ::TransportVelocityAdami, return system end - -# WARNING! -# These functions are intended to be used internally to set the transport velocity -# of newly activated particles in a callback. -# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` -# outside of callbacks. -function set_transport_velocity!(system::FluidSystem, particle, particle_old, v, v_old) - set_transport_velocity!(system, particle, particle_old, v, v_old, - system.transport_velocity) -end - -set_transport_velocity!(system, particle, particle_old, v, v_old) = system - -set_transport_velocity!(system, particle, particle_old, v, v_old, ::Nothing) = system - -function set_transport_velocity!(system, particle, particle_old, v, v_old, - ::TransportVelocityAdami) - for i in 1:ndims(system) - v[ndims(system) + i, particle] = v_old[i, particle_old] - end -end From 597b80fb51d7c3774f1c211d73d74b0281f05230 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 12 Aug 2024 11:49:12 +0200 Subject: [PATCH 44/62] clean up --- examples/fluid/taylor_green_vortex_2d.jl | 70 ++---------------------- 1 file changed, 6 insertions(+), 64 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 964245502..22c62edd3 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -80,75 +80,17 @@ semi = Semidiscretization(fluid_system, ode = semidiscretize(semi, tspan) -dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu)) - -function compute_L1v_error(v, u, t, system) - v_analytical_avg = 0.0 - L1v = 0.0 - - for particle in TrixiParticles.eachparticle(system) - position = TrixiParticles.current_coords(u, system, particle) - - v_mag = TrixiParticles.norm(TrixiParticles.current_velocity(v, system, particle)) - v_analytical = TrixiParticles.norm(velocity_function(position, t)) - - v_analytical_avg += abs(v_analytical) - L1v += abs(v_mag - v_analytical) - end - v_analytical_avg /= TrixiParticles.nparticles(system) - - L1v /= TrixiParticles.nparticles(system) - - return L1v /= v_analytical_avg -end - -function compute_L1p_error(v, u, t, system) - p_max_exact = 0.0 - - L1p = 0.0 - - for particle in TrixiParticles.eachparticle(system) - position = TrixiParticles.current_coords(u, system, particle) - - # compute pressure error - p_analytical = pressure_function(position, t) - p_max_exact = max(p_max_exact, abs(p_analytical)) +info_callback = InfoCallback(interval=100) - # p_computed - p_average - p_computed = TrixiParticles.particle_pressure(v, system, particle) - - TrixiParticles.average_pressure(system, particle) - L1p += abs(p_computed - p_analytical) - end +saving_callback = SolutionSavingCallback(dt=0.02) - L1p /= TrixiParticles.nparticles(system) +pp_callback = nothing - return L1p /= p_max_exact -end +callbacks = CallbackSet(info_callback, saving_callback, pp_callback, UpdateCallback()) -# The pressure plotted in the paper is the difference of the local pressure minus -# the average of the pressure of all particles. -function diff_p_loc_p_avg(v, u, t, system) - p_avg_tot = avg_pressure(v, u, t, system) - - return v[end, :] .- p_avg_tot -end +dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu)) -info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02, - L1v=compute_L1v_error, - L1p=compute_L1p_error, - diff_p_loc_p_avg=diff_p_loc_p_avg) - -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) - -# Use a Runge-Kutta method with automatic (error based) time step size control. -# Enable threading of the RK method for better performance on multiple threads. -# Limiting of the maximum stepsize is necessary to prevent crashing. -# When particles are approaching a wall in a uniform way, they can be advanced -# with large time steps. Close to the wall, the stepsize has to be reduced drastically. -# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces -# become extremely large when fluid particles are very close to boundary particles, -# and the time integration method interprets this as an instability. +# Use a Runge-Kutta method with automatic (error based) time step size control sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) From 24f946d599f3ed7a65295f300d3268c4bfdeb4ff Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 12 Aug 2024 14:13:14 +0200 Subject: [PATCH 45/62] add tgv validation example --- .../validation_taylor_green_vortex_2d.jl | 88 +++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl new file mode 100644 index 000000000..070e1265b --- /dev/null +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -0,0 +1,88 @@ +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacings = [0.02, 0.015, 0.01, 0.005] + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 5.0) +reynolds_number = 100.0 + +box_length = 1.0 + +function compute_L1v_error(v, u, t, system) + v_analytical_avg = 0.0 + L1v = 0.0 + + for particle in TrixiParticles.eachparticle(system) + position = TrixiParticles.current_coords(u, system, particle) + + v_mag = TrixiParticles.norm(TrixiParticles.current_velocity(v, system, particle)) + v_analytical = TrixiParticles.norm(velocity_function(position, t)) + + v_analytical_avg += abs(v_analytical) + L1v += abs(v_mag - v_analytical) + end + v_analytical_avg /= nparticles(system) + + L1v /= nparticles(system) + + return L1v /= v_analytical_avg +end + +function compute_L1p_error(v, u, t, system) + p_max_exact = 0.0 + + L1p = 0.0 + + for particle in TrixiParticles.eachparticle(system) + position = TrixiParticles.current_coords(u, system, particle) + + # compute pressure error + p_analytical = pressure_function(position, t) + p_max_exact = max(p_max_exact, abs(p_analytical)) + + # p_computed - p_average + p_computed = TrixiParticles.particle_pressure(v, system, particle) - + TrixiParticles.average_pressure(system, particle) + L1p += abs(p_computed - p_analytical) + end + + L1p /= nparticles(system) + + return L1p /= p_max_exact +end + +# The pressure plotted in the paper is the difference of the local pressure minus +# the average of the pressure of all particles. +function diff_p_loc_p_avg(v, u, t, system) + p_avg_tot = avg_pressure(v, u, t, system) + + return v[end, :] .- p_avg_tot +end + +for particle_spacing in particle_spacings + n_particles_xy = round(Int, box_length / particle_spacing) + + output_directory = joinpath("out_tgv", + "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)") + saving_callback = SolutionSavingCallback(dt=0.02, + output_directory=output_directory, + p_avg=diff_p_loc_p_avg) + + + pp_callback = PostprocessCallback(; dt=0.02, + L1v=compute_L1v_error, + L1p=compute_L1p_error, + output_directory=output_directory, + filename="errors", + write_csv=true, write_file_interval=1) + + # Import variables into scope + trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "taylor_green_vortex_2d.jl"), + particle_spacing=particle_spacing, reynolds_number=reynolds_number, + box_length=box_length, + saving_callback=saving_callback, pp_callback=pp_callback) +end From 2dd205ea78bdac035ab55e477fdb015014b8b048 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 12 Aug 2024 15:47:36 +0200 Subject: [PATCH 46/62] add ldc validation --- examples/fluid/lid_driven_cavity_2d.jl | 7 ++- .../validation_lid_driven_cavity_2d.jl | 44 +++++++++++++++++++ .../validation_taylor_green_vortex_2d.jl | 8 +--- 3 files changed, 51 insertions(+), 8 deletions(-) create mode 100644 validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index fd9064080..ee06d644b 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -13,7 +13,7 @@ using OrdinaryDiffEq particle_spacing = 0.02 # Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 3 +boundary_layers = 4 # ========================================================================================== # ==== Experiment Setup @@ -90,7 +90,10 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02) -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +pp_callback = nothing + +callbacks = CallbackSet(info_callback, saving_callback, pp_callback, UpdateCallback()) # Use a Runge-Kutta method with automatic (error based) time step size control sol = solve(ode, RDPK3SpFSAL49(), diff --git a/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl new file mode 100644 index 000000000..4bdfb2ede --- /dev/null +++ b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl @@ -0,0 +1,44 @@ +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacings = [0.02, 0.01, 0.005] + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 0.001) +reynolds_numbers = [100.0, 1000.0, 10_000.0] + +for particle_spacing in particle_spacings, reynolds_number in reynolds_numbers + n_particles_xy = round(Int, 1.0 / particle_spacing) + + Re = Int(reynolds_number) + + output_directory = joinpath("out_ldc", + "validation_run_lid_driven_cavity_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)_Re_$Re") + + saving_callback = SolutionSavingCallback(dt=0.02, output_directory=output_directory) + + # Import variables into scope + trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "lid_driven_cavity_2d.jl"), + saving_callback=saving_callback, tspan=tspan, + particle_spacing=particle_spacing, reynolds_number=reynolds_number) + + values_y = interpolate_line([0.5, 0.0], [0.5, 1.0], n_particles_xy, semi, fluid_system, + sol; endpoint=true, cut_off_bnd=true) + + values_x = interpolate_line([0.0, 0.5], [1.0, 0.5], n_particles_xy, semi, fluid_system, + sol; endpoint=true, cut_off_bnd=true) + + vy_y = stack(values_y.velocity)[2, :] + vy_x = stack(values_y.velocity)[1, :] + + vx_y = stack(values_x.velocity)[2, :] + vx_x = stack(values_x.velocity)[1, :] + + df = TrixiParticles.DataFrame(pos=collect(LinRange(0.0, 1.0, n_particles_xy)), + vy_y=vy_y, vy_x=vy_x, vx_y=vx_y, vx_x=vx_x) + + TrixiParticles.CSV.write(output_directory * "/interpolated_velocities.csv", df) +end diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 070e1265b..705a650fd 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -9,8 +9,6 @@ particle_spacings = [0.02, 0.015, 0.01, 0.005] tspan = (0.0, 5.0) reynolds_number = 100.0 -box_length = 1.0 - function compute_L1v_error(v, u, t, system) v_analytical_avg = 0.0 L1v = 0.0 @@ -63,7 +61,7 @@ function diff_p_loc_p_avg(v, u, t, system) end for particle_spacing in particle_spacings - n_particles_xy = round(Int, box_length / particle_spacing) + n_particles_xy = round(Int, 1.0 / particle_spacing) output_directory = joinpath("out_tgv", "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)") @@ -71,7 +69,6 @@ for particle_spacing in particle_spacings output_directory=output_directory, p_avg=diff_p_loc_p_avg) - pp_callback = PostprocessCallback(; dt=0.02, L1v=compute_L1v_error, L1p=compute_L1p_error, @@ -83,6 +80,5 @@ for particle_spacing in particle_spacings trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "taylor_green_vortex_2d.jl"), particle_spacing=particle_spacing, reynolds_number=reynolds_number, - box_length=box_length, - saving_callback=saving_callback, pp_callback=pp_callback) + tspan=tspan, saving_callback=saving_callback, pp_callback=pp_callback) end From 288e2d7977cd3b86c37dc32f8d60db2f14959d57 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 12 Aug 2024 17:01:14 +0200 Subject: [PATCH 47/62] modify validation --- .../validation_lid_driven_cavity_2d.jl | 81 +++++++++++++++---- 1 file changed, 65 insertions(+), 16 deletions(-) diff --git a/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl index 4bdfb2ede..376f730f5 100644 --- a/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl +++ b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl @@ -6,39 +6,88 @@ particle_spacings = [0.02, 0.01, 0.005] # ========================================================================================== # ==== Experiment Setup -tspan = (0.0, 0.001) +tspan = (0.0, 25.0) reynolds_numbers = [100.0, 1000.0, 10_000.0] +interpolated_velocity(v, u, t, system) = nothing + +function interpolated_velocity(v, u, t, system::TrixiParticles.FluidSystem) + if t < 24.8 + return nothing + end + + n_particles_xy = round(Int, 1.0 / system.initial_condition.particle_spacing) + + values_y = interpolate_line([0.5, 0.0], [0.5, 1.0], n_particles_xy, semi, system, + v, u; endpoint=true, cut_off_bnd=true) + + values_x = interpolate_line([0.0, 0.5], [1.0, 0.5], n_particles_xy, semi, system, + v, u; endpoint=true, cut_off_bnd=true) + + vy_y = stack(values_y.velocity)[2, :] + vy_x = stack(values_y.velocity)[1, :] + + vx_y = stack(values_x.velocity)[2, :] + vx_x = stack(values_x.velocity)[1, :] + + file = joinpath(output_directory, "interpolated_velocities.csv") + + if isfile(file) + data = TrixiParticles.CSV.read(file, TrixiParticles.DataFrame) + vy_y_ = (data.vy_y .+ vy_y) + vy_x_ = (data.vy_x .+ vy_x) + vx_y_ = (data.vx_y .+ vx_y) + vx_x_ = (data.vx_x .+ vx_x) + + n_evaluations = first(data.counter) + 1 + + df = TrixiParticles.DataFrame(pos=data.pos, counter=n_evaluations, + vy_y=vy_y_, vy_x=vy_x_, vx_y=vx_y_, vx_x=vx_x_) + + TrixiParticles.CSV.write(output_directory * "/interpolated_velocities.csv", df) + else + df = TrixiParticles.DataFrame(pos=collect(LinRange(0.0, 1.0, n_particles_xy)), + counter=1, vy_y=vy_y, vy_x=vy_x, vx_y=vx_y, vx_x=vx_x) + + TrixiParticles.CSV.write(output_directory * "/interpolated_velocities.csv", df) + end + + return nothing +end + for particle_spacing in particle_spacings, reynolds_number in reynolds_numbers n_particles_xy = round(Int, 1.0 / particle_spacing) Re = Int(reynolds_number) - output_directory = joinpath("out_ldc", - "validation_run_lid_driven_cavity_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)_Re_$Re") + global output_directory = joinpath("out_ldc", + "validation_run_lid_driven_cavity_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)_Re_$Re") - saving_callback = SolutionSavingCallback(dt=0.02, output_directory=output_directory) + saving_callback = SolutionSavingCallback(dt=0.5, output_directory=output_directory) + + pp_callback = PostprocessCallback(; dt=0.02, + interpolated_velocity=interpolated_velocity, + filename="interpolated_velocities", + write_file_interval=0) # Import variables into scope trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "lid_driven_cavity_2d.jl"), - saving_callback=saving_callback, tspan=tspan, + saving_callback=saving_callback, tspan=tspan, pp_callback=pp_callback, particle_spacing=particle_spacing, reynolds_number=reynolds_number) - values_y = interpolate_line([0.5, 0.0], [0.5, 1.0], n_particles_xy, semi, fluid_system, - sol; endpoint=true, cut_off_bnd=true) - - values_x = interpolate_line([0.0, 0.5], [1.0, 0.5], n_particles_xy, semi, fluid_system, - sol; endpoint=true, cut_off_bnd=true) + file = joinpath(output_directory, "interpolated_velocities.csv") - vy_y = stack(values_y.velocity)[2, :] - vy_x = stack(values_y.velocity)[1, :] + data = TrixiParticles.CSV.read(file, TrixiParticles.DataFrame) - vx_y = stack(values_x.velocity)[2, :] - vx_x = stack(values_x.velocity)[1, :] + n_evaluations = first(data.counter) - df = TrixiParticles.DataFrame(pos=collect(LinRange(0.0, 1.0, n_particles_xy)), - vy_y=vy_y, vy_x=vy_x, vx_y=vx_y, vx_x=vx_x) + df = TrixiParticles.DataFrame(pos=data.pos, + avg_vx_x=data.vx_x ./ n_evaluations, + avg_vx_y=data.vx_y ./ n_evaluations, + avg_vy_x=data.vy_x ./ n_evaluations, + avg_vy_y=data.vy_y ./ n_evaluations, + counter=n_evaluations) TrixiParticles.CSV.write(output_directory * "/interpolated_velocities.csv", df) end From 46a999eb596a1bfe7044a8a1e9600973ae2cd0aa Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 12 Aug 2024 19:17:30 +0200 Subject: [PATCH 48/62] increase `maxiters` --- examples/fluid/lid_driven_cavity_2d.jl | 1 + .../lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index ee06d644b..1a9aeb98a 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -100,4 +100,5 @@ sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration) dtmax=1e-2, # Limit stepsize to prevent crashing + maxiters=Int(1e7), save_everystep=false, callback=callbacks); diff --git a/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl index 376f730f5..4d7cb6b23 100644 --- a/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl +++ b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl @@ -63,7 +63,7 @@ for particle_spacing in particle_spacings, reynolds_number in reynolds_numbers global output_directory = joinpath("out_ldc", "validation_run_lid_driven_cavity_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)_Re_$Re") - saving_callback = SolutionSavingCallback(dt=0.5, output_directory=output_directory) + saving_callback = SolutionSavingCallback(dt=0.1, output_directory=output_directory) pp_callback = PostprocessCallback(; dt=0.02, interpolated_velocity=interpolated_velocity, From 1dfc4a32da72277211c961244ba0ac5396109859 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 Aug 2024 08:51:56 +0200 Subject: [PATCH 49/62] remove double velocity initialization --- src/schemes/fluid/transport_velocity.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 8080f4702..07505bb29 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -50,7 +50,6 @@ function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami) for particle in eachparticle(system) # Write particle velocities for dim in 1:ndims(system) - v0[dim, particle] = initial_condition.velocity[dim, particle] v0[ndims(system) + dim, particle] = initial_condition.velocity[dim, particle] end end From 93c944ce31b92211916196c6555eb6dfa201426b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 Aug 2024 09:30:22 +0200 Subject: [PATCH 50/62] add random displacement --- examples/fluid/taylor_green_vortex_2d.jl | 6 +----- src/TrixiParticles.jl | 1 + src/setups/rectangular_shape.jl | 10 ++++++++++ .../validation_taylor_green_vortex_2d.jl | 2 +- 4 files changed, 13 insertions(+), 6 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 22c62edd3..c11649aeb 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -59,14 +59,10 @@ smoothing_length = 1.0 * particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0), + coordinates_perturbation=0.2, # To avoid stagnant streamlines density=fluid_density, pressure=initial_pressure_function, velocity=initial_velocity_function) -# Add small random displacement to the particles to avoid stagnant streamlines. -#seed!(42); -#fluid.coordinates .+= rand((-particle_spacing / 5):1e-5:(particle_spacing / 5), -# size(fluid.coordinates)) - fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, transport_velocity=TransportVelocityAdami(background_pressure), diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 27b080215..6f4739cb0 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -19,6 +19,7 @@ using MuladdMacro: @muladd using Polyester: Polyester, @batch 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 @reexport using StaticArrays: SVector diff --git a/src/setups/rectangular_shape.jl b/src/setups/rectangular_shape.jl index d520bbfa5..db592ecc5 100644 --- a/src/setups/rectangular_shape.jl +++ b/src/setups/rectangular_shape.jl @@ -43,6 +43,8 @@ Rectangular shape filled with particles. Returns an [`InitialCondition`](@ref). - `tlsph`: With the [`TotalLagrangianSPHSystem`](@ref), particles need to be placed on the boundary of the shape and not one particle radius away, as for fluids. When `tlsph=true`, particles will be placed on the boundary of the shape. +- `coordinates_perturbation`: Add a small random displacement to the particle positions, + where the amplitude is `coordinates_perturbation * particle_spacing`. # Examples ```jldoctest; output = false, setup = :(particle_spacing = 0.1) @@ -64,6 +66,7 @@ InitialCondition{Float64}(0.1, [1.05 1.15 … 1.35 1.45; 2.05 2.05 … 2.35 2.35 """ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coordinates; velocity=zeros(length(n_particles_per_dimension)), + coordinates_perturbation=nothing, mass=nothing, density=nothing, pressure=0.0, acceleration=nothing, state_equation=nothing, tlsph=false, loop_order=nothing) @@ -89,6 +92,13 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord min_coordinates, tlsph=tlsph, loop_order=loop_order) + if !isnothing(coordinates_perturbation) + seed!(1) + amplitude = coordinates_perturbation * particle_spacing + coordinates .+= rand((-amplitude):(particle_spacing * 1e-3):(amplitude), + NDIMS, n_particles) + end + # Allow zero acceleration with state equation, but interpret `nothing` acceleration # with state equation as a likely mistake. if acceleration isa AbstractVector || acceleration isa Tuple diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 705a650fd..337fd9cb0 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -2,7 +2,7 @@ using TrixiParticles # ========================================================================================== # ==== Resolution -particle_spacings = [0.02, 0.015, 0.01, 0.005] +particle_spacings = [0.02, 0.01, 0.005] # ========================================================================================== # ==== Experiment Setup From b6e734602e164f8c64d922041656fc17f1e4c47e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 19 Sep 2024 16:26:20 +0200 Subject: [PATCH 51/62] implement suggestions --- docs/src/systems/entropically_damped_sph.md | 4 ++-- examples/fluid/lid_driven_cavity_2d.jl | 6 +++--- src/schemes/fluid/entropically_damped_sph/rhs.jl | 1 + src/schemes/fluid/entropically_damped_sph/system.jl | 2 +- src/schemes/fluid/fluid.jl | 2 +- 5 files changed, 8 insertions(+), 7 deletions(-) diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 5734c2214..d7bb711fa 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -82,10 +82,10 @@ The discretized form of the last term is ``` where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. -Note that although ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, +Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, which means that there is a non-vanishing contribution only when particles are disordered. That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions. -Suggested is a background pressure which is on the order of the reference pressure but can be chosen arbitrarily large when time-step criterion is adjusted. +Suggested is a background pressure which is in the order of the reference pressure but can be chosen arbitrarily large when the time-step criterion is adjusted. The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index 1a9aeb98a..77b409588 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -55,11 +55,11 @@ fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel, smoot # ========================================================================================== # ==== Boundary -movement_function(t) = SVector(VELOCITY_LID * t, 0.0) +lid_movement_function(t) = SVector(VELOCITY_LID * t, 0.0) is_moving(t) = true -movement = BoundaryMovement(movement_function, is_moving) +lid_movement = BoundaryMovement(lid_movement_function, is_moving) boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density, cavity.boundary.mass, @@ -74,7 +74,7 @@ boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass, boundary_system_cavity = BoundarySPHSystem(cavity.boundary, boundary_model_cavity) -boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=movement) +boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=lid_movement) # ========================================================================================== # ==== Simulation diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 4ac10709d..3644c3630 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -24,6 +24,7 @@ function interact!(dv, v_particle_system, u_particle_system, # This technique is for a more robust `pressure_acceleration` but only with TVF. # It results only in significant improvement for EDAC and not for WCSPH. # See Ramachandran (2019) p. 582 + # Note that the return value is zero when not using EDAC with TVF. p_avg = average_pressure(particle_system, particle) m_a = hydrodynamic_mass(particle_system, particle) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index e1320d82e..d93ae6982 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -215,7 +215,7 @@ end # This technique is for a more robust `pressure_acceleration` but only with TVF. # It results only in significant improvement for EDAC and not for WCSPH. -# See Ramachandran (2019) p. 582 +# See Ramachandran (2019) p. 582. function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode, semi) (; cache) = system (; pressure_average, neighbor_counter) = cache diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 995d87856..e0f07796c 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -92,7 +92,7 @@ end @inline function momentum_convection(system, neighbor_system, 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)))) + return zero(grad_kernel) end @inline function momentum_convection(system, From 196cebfe81c21f66365b944fb373598128bde842 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 4 Oct 2024 16:06:41 +0200 Subject: [PATCH 52/62] update `NEWS.md` --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index 8b1c87cde..e60c63e9f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,11 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Version 0.2.3 + +### Highlights +Transport Velocity Formulation (TVF) based on the work of Ramachandran et al. "Entropically damped artificial compressibility for SPH" (2019) was added. + ## Version 0.2.2 ### Highlights From 66241b2f954f593ecbd18118497525a307ba13cd Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 4 Oct 2024 16:32:25 +0200 Subject: [PATCH 53/62] remove check for viscosity --- src/general/semidiscretization.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index d7ab3d6b4..600d573d5 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -854,12 +854,3 @@ function check_configuration(system::TotalLagrangianSPHSystem, systems) "`ContinuityDensity` is not yet supported for a `TotalLagrangianSPHSystem`")) end end - -function check_configuration(system::FluidSystem, systems) - (; viscosity) = system - - if viscosity isa ArtificialViscosityMonaghan && - system.transport_velocity isa TransportVelocityAdami - throw(ArgumentError("Please use `ViscosityAdami` when simulating with `TransportVelocityAdami`")) - end -end From 4b05674d7fd47d74045dd980485c30f830c86e85 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 4 Oct 2024 16:32:31 +0200 Subject: [PATCH 54/62] fix tests --- test/systems/edac_system.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 2b34ff2d9..dd6209f79 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -217,6 +217,8 @@ smoothing_length, 0.0) semi = Semidiscretization(system) + TrixiParticles.initialize_neighborhood_searches!(semi) + u_ode = vec(fluid.coordinates) v_ode = vec(vcat(fluid.velocity, fluid.velocity, fluid.pressure')) From e22702dcece6113ed1515567e63890e1b089927a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 4 Oct 2024 16:43:46 +0200 Subject: [PATCH 55/62] add short tilde description in docs --- docs/src/systems/entropically_damped_sph.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index d7bb711fa..2f285df71 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -74,6 +74,7 @@ and is obtained at every time-step ``\Delta t`` from ``` where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field. +The tilde in the second term of the right hand side indicates that the material derivative has an advection part. The discretized form of the last term is From a98e2781635cd418169ae3e1e2f10cb989e1335a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 7 Oct 2024 16:30:40 +0200 Subject: [PATCH 56/62] remove discarded example from tests --- test/examples/examples.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 1441f76e5..d4bbf5985 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -235,15 +235,6 @@ @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/dam_break_2d_surface_tension.jl" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "dam_break_2d_surface_tension.jl"), - tspan=(0.0, 0.1)) - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - @trixi_testset "fluid/taylor_green_vortex_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", From f1a1d48cb6ec7d265cb2e24f6917b5ba4e2e2893 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 7 Oct 2024 16:33:36 +0200 Subject: [PATCH 57/62] implement suggestions --- docs/src/systems/entropically_damped_sph.md | 2 +- src/schemes/fluid/entropically_damped_sph/system.jl | 8 +++----- src/schemes/fluid/transport_velocity.jl | 7 ++----- 3 files changed, 6 insertions(+), 11 deletions(-) diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 2f285df71..e81de5cad 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -109,7 +109,7 @@ Here, ``\tilde{p}_{ab}`` is the density-weighted pressure \tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, ``` -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and is given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. +with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and are given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. ```@autodocs Modules = [TrixiParticles] diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index d93ae6982..614190934 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -246,11 +246,9 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode end end - for particle in eachparticle(system) - if neighbor_counter[particle] > 0 - pressure_average[particle] /= neighbor_counter[particle] - end - end + pressure_average ./= neighbor_counter + + return system end function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::SummationDensity) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 07505bb29..26a3688b6 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -13,11 +13,8 @@ See [TVF](@ref transport_velocity_formulation) for more details of the method. Thus, it is highly recommended to use [`ViscosityAdami`](@ref) as viscosity model, since [`ArtificialViscosityMonaghan`](@ref) leads to bad results. """ -struct TransportVelocityAdami{ELTYPE} - background_pressure::ELTYPE - function TransportVelocityAdami(background_pressure::Real) - new{typeof(background_pressure)}(background_pressure) - end +struct TransportVelocityAdami{T<:Real} + background_pressure::T end # Calculate `v_nvariables` appropriately From d963af0564a4e87349461fb2d94366ac11d87d1f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 16 Oct 2024 11:03:11 +0200 Subject: [PATCH 58/62] apply formatter --- src/schemes/fluid/transport_velocity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 26a3688b6..bc8fdf005 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -13,7 +13,7 @@ See [TVF](@ref transport_velocity_formulation) for more details of the method. Thus, it is highly recommended to use [`ViscosityAdami`](@ref) as viscosity model, since [`ArtificialViscosityMonaghan`](@ref) leads to bad results. """ -struct TransportVelocityAdami{T<:Real} +struct TransportVelocityAdami{T <: Real} background_pressure::T end From 0589a04b04c718b7b3f2b87d5d529dae93d797c7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 23 Oct 2024 10:36:34 +0200 Subject: [PATCH 59/62] implement suggestions --- examples/fluid/periodic_array_of_cylinders_2d.jl | 5 ++--- examples/fluid/taylor_green_vortex_2d.jl | 5 ++--- test/schemes/fluid/transport_velocity.jl | 6 +++--- test/systems/edac_system.jl | 2 +- .../validation_taylor_green_vortex_2d.jl | 16 ++++++++-------- 5 files changed, 16 insertions(+), 18 deletions(-) diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 2734a8529..a8ecb035c 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -1,6 +1,6 @@ # Channel flow through periodic array of cylinders # -# S. Adami et al +# S. Adami et al. # "A transport-velocity formulation for smoothed particle hydrodynamics". # In: Journal of Computational Physics, Volume 241 (2013), pages 292-307. # https://doi.org/10.1016/j.jcp.2013.01.043 @@ -19,13 +19,12 @@ boundary_layers = 3 # ==== Experiment Setup tspan = (0.0, 5.0) -const acceleration_x = 2.5e-4 +acceleration_x = 2.5e-4 # Boundary geometry and initial fluid particle positions cylinder_radius = 0.02 tank_size = (6 * cylinder_radius, 4 * cylinder_radius) fluid_size = tank_size -initial_velocity = (1.2e-4, 0.0) fluid_density = 1000.0 nu = 0.1 / fluid_density # viscosity parameter diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index c11649aeb..551863ccb 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -59,7 +59,7 @@ smoothing_length = 1.0 * particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0), - coordinates_perturbation=0.2, # To avoid stagnant streamlines + coordinates_perturbation=0.2, # To avoid stagnant streamlines when not using TVF. density=fluid_density, pressure=initial_pressure_function, velocity=initial_velocity_function) @@ -90,5 +90,4 @@ dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=dt_max,#1e-2, # Limit stepsize to prevent crashing - save_everystep=false, callback=callbacks); + dtmax=dt_max, save_everystep=false, callback=callbacks); diff --git a/test/schemes/fluid/transport_velocity.jl b/test/schemes/fluid/transport_velocity.jl index f53cd813a..bfa71e125 100644 --- a/test/schemes/fluid/transport_velocity.jl +++ b/test/schemes/fluid/transport_velocity.jl @@ -3,7 +3,7 @@ smoothing_kernel = SchoenbergCubicSplineKernel{2}() smoothing_length = 1.2particle_spacing - fluid = rectangular_patch(particle_spacing, (3, 3), seed=1) + fluid = rectangular_patch(particle_spacing, (3, 3)) v0_tvf = zeros(5, nparticles(fluid)) @@ -20,7 +20,7 @@ @testset "write_v0!" begin TrixiParticles.write_v0!(v0_tvf, system_tvf) - @test vcat(fluid.velocity, fluid.velocity, fluid.pressure') ≈ v0_tvf + @test isapprox(vcat(fluid.velocity, fluid.velocity, fluid.pressure'), v0_tvf) end @testset "Update" begin @@ -30,6 +30,6 @@ TrixiParticles.update_transport_velocity!(system_tvf, vec(v0_tvf), semi) - @test fill(2.5, (4, nparticles(system_tvf))) ≈ v0_tvf[1:4, :] + @test isapprox(fill(2.5, (4, nparticles(system_tvf))), v0_tvf[1:4, :]) end end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index dd6209f79..9409e4cf2 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -227,7 +227,7 @@ @test all(i -> system.cache.neighbor_counter[i] == nparticles(system), nparticles(system)) - @test all(i -> system.cache.pressure_average[i] ≈ -50.968532955185964, + @test all(isapprox(i -> system.cache.pressure_average[i], -50.968532955185964), nparticles(system)) end end diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 337fd9cb0..cd9d18159 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -9,9 +9,9 @@ particle_spacings = [0.02, 0.01, 0.005] tspan = (0.0, 5.0) reynolds_number = 100.0 -function compute_L1v_error(v, u, t, system) +function compute_l1v_error(v, u, t, system) v_analytical_avg = 0.0 - L1v = 0.0 + v_avg = 0.0 for particle in TrixiParticles.eachparticle(system) position = TrixiParticles.current_coords(u, system, particle) @@ -20,16 +20,16 @@ function compute_L1v_error(v, u, t, system) v_analytical = TrixiParticles.norm(velocity_function(position, t)) v_analytical_avg += abs(v_analytical) - L1v += abs(v_mag - v_analytical) + v_avg += abs(v_mag - v_analytical) end v_analytical_avg /= nparticles(system) - L1v /= nparticles(system) + v_avg /= nparticles(system) - return L1v /= v_analytical_avg + return v_avg /= v_analytical_avg end -function compute_L1p_error(v, u, t, system) +function compute_l1p_error(v, u, t, system) p_max_exact = 0.0 L1p = 0.0 @@ -70,8 +70,8 @@ for particle_spacing in particle_spacings p_avg=diff_p_loc_p_avg) pp_callback = PostprocessCallback(; dt=0.02, - L1v=compute_L1v_error, - L1p=compute_L1p_error, + L1v=compute_l1v_error, + L1p=compute_l1p_error, output_directory=output_directory, filename="errors", write_csv=true, write_file_interval=1) From 605bb20c4edc65f104cd1a4c042db944009d55f6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 29 Oct 2024 18:09:11 +0100 Subject: [PATCH 60/62] add comment --- src/schemes/fluid/entropically_damped_sph/system.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 614190934..eb2aa1c1d 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -246,6 +246,8 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode end end + # We do not need to check for zero division here, as `neighbor_counter = 1` + # for zero neighbors. That is, the `particle` itself is also taken into account. pressure_average ./= neighbor_counter return system From d776e62eeddd3a50a58d8d97cb8474a39f09fb9a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 29 Oct 2024 18:17:08 +0100 Subject: [PATCH 61/62] fix test --- test/systems/edac_system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 9409e4cf2..ea188fed9 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -227,7 +227,7 @@ @test all(i -> system.cache.neighbor_counter[i] == nparticles(system), nparticles(system)) - @test all(isapprox(i -> system.cache.pressure_average[i], -50.968532955185964), + @test all(i -> isapprox(system.cache.pressure_average[i], -50.968532955185964), nparticles(system)) end end From b5e0e981cdd4879b293276a215195322ae64fc8d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 29 Oct 2024 19:29:35 +0100 Subject: [PATCH 62/62] Set to Version 0.2.3 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e00f645f3..956cdbad5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TrixiParticles" uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a" authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"] -version = "0.2.3-dev" +version = "0.2.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"