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 diff --git a/Project.toml b/Project.toml index c632a1994..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" @@ -21,6 +21,7 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 86df0782f..e1b287b07 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -13,6 +13,20 @@ @Article{Adami2012 } +@Article{Adami2013, + author = {S. Adami and X.Y. Hu and N.A. Adams}, + journal = {Journal of Computational Physics}, + title = {A transport-velocity formulation for smoothed particle hydrodynamics}, + year = {2013}, + month = {may}, + pages = {292--307}, + volume = {241}, + doi = {10.1016/j.jcp.2013.01.043}, + groups = {enhancement}, + publisher = {Elsevier {BV}}, +} + + @Article{Akinci2012, author = {Akinci, Nadir and Ihmsen, Markus and Akinci, Gizem and Solenthaler, Barbara and Teschner, Matthias}, journal = {ACM Transactions on Graphics}, diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index d035fb6fc..0db8d275e 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -45,3 +45,64 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). Modules = [TrixiParticles] 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. +To address these problems, [Adami (2013)](@cite Adami2013) 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 (2019)](@cite Ramachandran2019) 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 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 + +```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}, +``` + +where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. +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 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 + +```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)^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]. +``` + +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. ``\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] +Pages = [joinpath("schemes", "fluid", "transport_velocity.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..77b409588 --- /dev/null +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -0,0 +1,104 @@ +# 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 + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.02 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 4 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 5.0) +reynolds_number = 100.0 + +cavity_size = (1.0, 1.0) + +fluid_density = 1.0 + +const VELOCITY_LID = 1.0 +sound_speed = 10 * VELOCITY_LID + +pressure = sound_speed^2 * fluid_density + +viscosity = ViscosityAdami(; nu=VELOCITY_LID / reynolds_number) + +cavity = RectangularTank(particle_spacing, cavity_size, cavity_size, fluid_density, + n_layers=boundary_layers, + faces=(true, true, true, false), pressure=pressure) + +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]), density=fluid_density) + +# ========================================================================================== +# ==== Fluid + +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 + +lid_movement_function(t) = SVector(VELOCITY_LID * t, 0.0) + +is_moving(t) = true + +lid_movement = BoundaryMovement(lid_movement_function, is_moving) + +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) + +boundary_system_cavity = BoundarySPHSystem(cavity.boundary, boundary_model_cavity) + +boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=lid_movement) + +# ========================================================================================== +# ==== 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{2}(; periodic_box)) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) + +saving_callback = SolutionSavingCallback(dt=0.02) + +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(), + 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/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl new file mode 100644 index 000000000..a8ecb035c --- /dev/null +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -0,0 +1,87 @@ +# 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 + +# ========================================================================================== +# ==== Resolution +n_particles_x = 144 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 5.0) + +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 + +fluid_density = 1000.0 +nu = 0.1 / fluid_density # viscosity parameter + +# 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 + +particle_spacing = tank_size[1] / n_particles_x + +box = RectangularTank(particle_spacing, fluid_size, tank_size, + fluid_density, n_layers=boundary_layers, + pressure=pressure, faces=(false, false, true, true)) + +cylinder = SphereShape(particle_spacing, cylinder_radius, tank_size ./ 2, + fluid_density, sphere_type=RoundSphere()) + +fluid = setdiff(box.fluid, cylinder) +boundary = union(cylinder, box.boundary) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 1.2 * particle_spacing +smoothing_kernel = SchoenbergQuarticSplineKernel{2}() +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + 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 +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{2}(; periodic_box)) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +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 +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 + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl new file mode 100644 index 000000000..551863ccb --- /dev/null +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -0,0 +1,93 @@ +# 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 + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.02 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 5.0) +reynolds_number = 100.0 + +box_length = 1.0 + +U = 1.0 # m/s +fluid_density = 1.0 +sound_speed = 10U + +b = -8pi^2 / reynolds_number + +# Taylor Green Vortex Pressure Function +function pressure_function(pos, t) + x = pos[1] + y = pos[2] + + return -U^2 * exp(2 * b * t) * (cos(4pi * x) + cos(4pi * y)) / 4 +end + +initial_pressure_function(pos) = pressure_function(pos, 0.0) + +# Taylor Green Vortex Velocity Function +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 + +initial_velocity_function(pos) = velocity_function(pos, 0.0) + +n_particles_xy = round(Int, box_length / particle_spacing) + +# ========================================================================================== +# ==== Fluid +nu = U * box_length / reynolds_number + +background_pressure = sound_speed^2 * fluid_density + +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 when not using TVF. + density=fluid_density, pressure=initial_pressure_function, + velocity=initial_velocity_function) + +fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + sound_speed, + transport_velocity=TransportVelocityAdami(background_pressure), + viscosity=ViscosityAdami(; nu)) + +# ========================================================================================== +# ==== Simulation +periodic_box = PeriodicBox(min_corner=[0.0, 0.0], max_corner=[box_length, box_length]) +semi = Semidiscretization(fluid_system, + neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box)) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) + +saving_callback = SolutionSavingCallback(dt=0.02) + +pp_callback = nothing + +callbacks = CallbackSet(info_callback, saving_callback, pp_callback, UpdateCallback()) + +dt_max = min(smoothing_length / 4 * (sound_speed + U), smoothing_length^2 / (8 * nu)) + +# 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) + dtmax=dt_max, save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d6e0c2996..bcdce98ac 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 @@ -59,7 +60,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback export ContinuityDensity, SummationDensity -export PenaltyForceGanzenmueller +export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 3595d5c36..821d36e33 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -80,11 +80,15 @@ 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 + @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) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 56e7746fc..9d2c79b00 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -292,8 +292,6 @@ end v_new[dim, particle_new] = v_old[dim, particle_old] end - # TODO: Only when using TVF: set tvf - return system_new end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index be8d689f5..3644c3630 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -21,22 +21,35 @@ 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) + # 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) m_b = hydrodynamic_mass(neighbor_system, neighbor) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) dv_pressure = pressure_acceleration(particle_system, neighbor_system, neighbor, - m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) + m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, + rho_b, pos_diff, distance, grad_kernel, + correction) dv_viscosity_ = dv_viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + # Add convection term when using `TransportVelocityAdami` + dv_convection = momentum_convection(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - @@ -46,6 +59,9 @@ function interact!(dv, v_particle_system, u_particle_system, particle, pos_diff, distance, sound_speed, m_a, m_b, p_a, p_b, rho_a, rho_b) + transport_velocity!(dv, particle_system, rho_a, rho_b, m_a, m_b, + grad_kernel, particle) + continuity_equation!(dv, density_calculator, v_diff, particle, m_b, rho_a, rho_b, particle_system, grad_kernel) end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index cf73f5876..eb2aa1c1d 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -3,6 +3,7 @@ smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), + transport_velocity=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), buffer_size=nothing, source_terms=nothing) @@ -28,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. - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `source_terms`: Additional source terms for this system. Has to be either `nothing` @@ -41,7 +43,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more The keyword argument `acceleration` should be used instead for gravity-like source terms. """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, PF, ST, B, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] @@ -54,6 +56,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, acceleration :: SVector{NDIMS, ELTYPE} correction :: Nothing pressure_acceleration_formulation :: PF + transport_velocity :: TV source_terms :: ST buffer :: B cache :: C @@ -62,6 +65,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), + transport_velocity=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), @@ -93,14 +97,15 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = create_cache_density(initial_condition, density_calculator) + cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), - typeof(density_calculator), typeof(smoothing_kernel), - typeof(viscosity), typeof(pressure_acceleration), typeof(source_terms), + typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), + typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms), typeof(buffer), typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, - acceleration_, nothing, pressure_acceleration, source_terms, + smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, + nothing, pressure_acceleration, transport_velocity, source_terms, buffer, cache) end end @@ -134,21 +139,33 @@ 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 end +create_cache_edac(initial_condition, ::Nothing) = (;) + +function create_cache_edac(initial_condition, ::TransportVelocityAdami) + pressure_average = copy(initial_condition.pressure) + neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) + update_callback_used = Ref(false) + + return (; pressure_average, neighbor_counter, update_callback_used) +end + @inline function v_nvariables(system::EntropicallyDampedSPHSystem) return v_nvariables(system, system.density_calculator) end @inline function v_nvariables(system::EntropicallyDampedSPHSystem, density_calculator) - return ndims(system) + 1 + return ndims(system) * factor_tvf(system) + 1 end @inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::ContinuityDensity) - return ndims(system) + 2 + return ndims(system) * factor_tvf(system) + 2 end @inline function particle_density(v, ::ContinuityDensity, @@ -174,12 +191,69 @@ end @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed +@inline average_pressure(system, particle) = zero(eltype(system)) + +@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) = zero(eltype(system)) + function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) compute_density!(system, u, u_ode, semi, system.density_calculator) + update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi) +end + +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 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 + (; pressure_average, neighbor_counter) = cache + + set_zero!(pressure_average) + set_zero!(neighbor_counter) + + u = wrap_u(u_ode, system, semi) + + # Use all other systems for the average pressure + @trixi_timeit timer() "compute average pressure" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, + pos_diff, distance + pressure_average[particle] += particle_pressure(v_neighbor_system, + neighbor_system, + neighbor) + neighbor_counter[particle] += 1 + 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 end -function write_v0!(v0, system::EntropicallyDampedSPHSystem, density_calculator) +function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::SummationDensity) for particle in eachparticle(system) v0[end, particle] = system.initial_condition.pressure[particle] end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 9e2289a8c..e0f07796c 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -31,11 +31,12 @@ function write_v0!(v0, system::FluidSystem) copyto!(v0, indices, system.initial_condition.velocity, indices) write_v0!(v0, system, system.density_calculator) + write_v0!(v0, system, system.transport_velocity) return v0 end -write_v0!(v0, system, density_calculator) = v0 +write_v0!(v0, system::FluidSystem, _) = v0 # To account for boundary effects in the viscosity term of the RHS, use the viscosity model # of the neighboring particle systems. @@ -77,6 +78,29 @@ end include("pressure_acceleration.jl") include("viscosity.jl") +include("transport_velocity.jl") include("surface_tension.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 zero(grad_kernel) +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 new file mode 100644 index 000000000..bc8fdf005 --- /dev/null +++ b/src/schemes/fluid/transport_velocity.jl @@ -0,0 +1,175 @@ +""" + TransportVelocityAdami(background_pressure::Real) + +Transport Velocity Formulation (TVF) to suppress pairing and tensile instability. +See [TVF](@ref transport_velocity_formulation) for more details of 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{T <: Real} + background_pressure::T +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) + update_transport_velocity!(system, v_ode, semi, system.transport_velocity) +end + +@inline update_transport_velocity!(system, v_ode, semi, ::Nothing) = system + +@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] + end + end + + return system +end + +function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami) + (; initial_condition) = system + + for particle in eachparticle(system) + # Write particle velocities + for dim in 1:ndims(system) + v0[ndims(system) + dim, particle] = initial_condition.velocity[dim, particle] + end + end + + return v0 +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, particle) + return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) +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) + 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, + 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) + + 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 transport_velocity!(dv, system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) = dv + +@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 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, + 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 + + return dv +end + +function reset_callback_flag!(system::FluidSystem) + reset_callback_flag!(system, system.transport_velocity) +end + +reset_callback_flag!(system, ::Nothing) = system + +function reset_callback_flag!(system::FluidSystem, ::TransportVelocityAdami) + system.cache.update_callback_used[] = false + + return system +end + +function update_callback_used!(system::FluidSystem) + update_callback_used!(system, system.transport_velocity) +end + +update_callback_used!(system, ::Nothing) = system + +function update_callback_used!(system, transport_velocity) + system.cache.update_callback_used[] = true +end + +function update_final!(system::FluidSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + update_final!(system, system.transport_velocity, + v, u, v_ode, u_ode, semi, t; update_from_callback) +end + +function update_final!(system::FluidSystem, ::Nothing, + v, u, v_ode, u_ode, semi, t; update_from_callback=false) + return system +end + +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`")) + end + + return system +end diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 482261fc1..eb6c06c84 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -58,6 +58,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF + transport_velocity :: Nothing # TODO source_terms :: ST surface_tension :: SRFT buffer :: B @@ -121,7 +122,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, smoothing_kernel, smoothing_length, acceleration_, viscosity, density_diffusion, correction, - pressure_acceleration, + pressure_acceleration, nothing, source_terms, surface_tension, buffer, cache) end 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/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 63d7628f4..ef5f68cc6 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -270,6 +270,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/examples/examples.jl b/test/examples/examples.jl index c962b5912..d4bbf5985 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -226,6 +226,24 @@ @test count_rhs_allocations(sol, semi) == 0 end + @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)) + @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)) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/sphere_surface_tension_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", @@ -234,6 +252,15 @@ @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)) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/sphere_surface_tension_3d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index 639fa6429..cfb75a5c7 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -1,5 +1,6 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") +include("transport_velocity.jl") include("surface_tension.jl") include("viscosity.jl") diff --git a/test/schemes/fluid/transport_velocity.jl b/test/schemes/fluid/transport_velocity.jl new file mode 100644 index 000000000..bfa71e125 --- /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)) + + 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 isapprox(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 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 33b06384d..ea188fed9 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -32,6 +32,7 @@ @test system.mass == mass @test system.smoothing_kernel == smoothing_kernel @test system.smoothing_length == smoothing_length + @test system.transport_velocity isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -86,6 +87,7 @@ @test system.mass == setup.mass @test system.smoothing_kernel == smoothing_kernel @test system.smoothing_length == smoothing_length + @test system.transport_velocity isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -137,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 @@ -201,4 +204,30 @@ @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) + + TrixiParticles.initialize_neighborhood_searches!(semi) + + 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 -> isapprox(system.cache.pressure_average[i], -50.968532955185964), + nparticles(system)) + end end 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..4d7cb6b23 --- /dev/null +++ b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl @@ -0,0 +1,93 @@ +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacings = [0.02, 0.01, 0.005] + +# ========================================================================================== +# ==== Experiment Setup +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) + + 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.1, 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, pp_callback=pp_callback, + particle_spacing=particle_spacing, reynolds_number=reynolds_number) + + file = joinpath(output_directory, "interpolated_velocities.csv") + + data = TrixiParticles.CSV.read(file, TrixiParticles.DataFrame) + + n_evaluations = first(data.counter) + + 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 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..cd9d18159 --- /dev/null +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -0,0 +1,84 @@ +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacings = [0.02, 0.01, 0.005] + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 5.0) +reynolds_number = 100.0 + +function compute_l1v_error(v, u, t, system) + v_analytical_avg = 0.0 + v_avg = 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) + v_avg += abs(v_mag - v_analytical) + end + v_analytical_avg /= nparticles(system) + + v_avg /= nparticles(system) + + return v_avg /= 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, 1.0 / 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, + tspan=tspan, saving_callback=saving_callback, pp_callback=pp_callback) +end