trixi-framework · LasNikas · Apr 10, 2025
diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl
@@ -77,15 +77,17 @@ if wcsph
                                        exponent=1)
     fluid_system = WeaklyCompressibleSPHSystem(fluid, density_calculator,
                                                state_equation, smoothing_kernel,
+                                               density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),
                                                pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure,
                                                smoothing_length,
-                                               viscosity=ViscosityAdami(; nu),
-                                               transport_velocity=TransportVelocityAdami(background_pressure))
+                                               transport_velocity=TransportVelocityAdami(background_pressure),
+                                               viscosity=ViscosityAdami(; nu))
 else
     density_calculator = SummationDensity()
     fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
                                                sound_speed,
                                                density_calculator=density_calculator,
+                                            #    density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),
                                                transport_velocity=TransportVelocityAdami(background_pressure),
                                                viscosity=ViscosityAdami(; nu))
 end

diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -79,8 +79,10 @@ function interact!(dv, v_particle_system, u_particle_system,
         transport_velocity!(dv, particle_system, neighbor_system, particle, neighbor,
                             m_a, m_b, grad_kernel)
 
-        continuity_equation!(dv, density_calculator, v_diff, particle, m_b, rho_a, rho_b,
-                             particle_system, grad_kernel)
+        continuity_equation!(dv, density_calculator, v_particle_system,
+                             v_neighbor_system, particle, neighbor,
+                             pos_diff, distance, m_b, rho_a, rho_b,
+                             particle_system, neighbor_system, grad_kernel)
     end
 
     return dv
@@ -125,20 +127,3 @@ end
 
     return dv
 end
-
-# We need a separate method for EDAC since the density is stored in `v[end-1,:]`.
-@inline function continuity_equation!(dv, density_calculator::ContinuityDensity,
-                                      vdiff, particle, m_b, rho_a, rho_b,
-                                      particle_system::EntropicallyDampedSPHSystem,
-                                      grad_kernel)
-    dv[end - 1, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel)
-
-    return dv
-end
-
-@inline function continuity_equation!(dv, density_calculator,
-                                      vdiff, particle, m_b, rho_a, rho_b,
-                                      particle_system::EntropicallyDampedSPHSystem,
-                                      grad_kernel)
-    return dv
-end
diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl
@@ -3,7 +3,7 @@
                                 smoothing_length, sound_speed;
                                 pressure_acceleration=inter_particle_averaged_pressure,
                                 density_calculator=SummationDensity(),
-                                transport_velocity=nothing,
+                                transport_velocity=nothing, density_diffusion=nothing,
                                 alpha=0.5, viscosity=nothing,
                                 acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing,
                                 surface_normal_method=nothing, buffer_size=nothing,
@@ -31,6 +31,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                                 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.
+- `density_diffusion`:          Density diffusion terms for this system. See [`DensityDiffusion`](@ref).
 - `buffer_size`:                Number of buffer particles.
                                 This is needed when simulating with [`OpenBoundarySPHSystem`](@ref).
 - `correction`:                 Correction method used for this system. (default: no correction, see [Corrections](@ref corrections))
@@ -52,7 +53,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
 - `color_value`:                The value used to initialize the color of particles in the system.
 
 """
-struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
+struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, DD,
                                    PF, ST, SRFT, SRFN, B, PR, C} <: FluidSystem{NDIMS}
     initial_condition                 :: IC
     mass                              :: M # Vector{ELTYPE}: [particle]
@@ -65,6 +66,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
     correction                        :: Nothing
     pressure_acceleration_formulation :: PF
     transport_velocity                :: TV
+    density_diffusion                 :: DD
     source_terms                      :: ST
     surface_tension                   :: SRFT
     surface_normal_method             :: SRFN
@@ -79,7 +81,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
                                      smoothing_length, sound_speed;
                                      pressure_acceleration=inter_particle_averaged_pressure,
                                      density_calculator=SummationDensity(),
-                                     transport_velocity=nothing,
+                                     transport_velocity=nothing, density_diffusion=nothing,
                                      alpha=0.5, viscosity=nothing,
                                      acceleration=ntuple(_ -> 0.0,
                                                          ndims(smoothing_kernel)),
@@ -152,16 +154,17 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
     EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass),
                                 typeof(density_calculator), typeof(smoothing_kernel),
                                 typeof(viscosity), typeof(transport_velocity),
+                                typeof(density_diffusion),
                                 typeof(pressure_acceleration), typeof(source_terms),
                                 typeof(surface_tension), typeof(surface_normal_method),
                                 typeof(buffer), Nothing,
                                 typeof(cache)}(initial_condition, mass, density_calculator,
                                                smoothing_kernel, sound_speed, viscosity,
                                                nu_edac, acceleration_, correction,
                                                pressure_acceleration, transport_velocity,
-                                               source_terms, surface_tension,
-                                               surface_normal_method, buffer,
-                                               particle_refinement, cache)
+                                               density_diffusion, source_terms,
+                                               surface_tension, surface_normal_method,
+                                               buffer, particle_refinement, cache)
 end
 
 function Base.show(io::IO, system::EntropicallyDampedSPHSystem)
@@ -200,6 +203,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst
         summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof)
         summary_line(io, "tansport velocity formulation",
                      system.transport_velocity |> typeof |> nameof)
+        summary_line(io, "density diffusion", system.density_diffusion)
         summary_line(io, "acceleration", system.acceleration)
         summary_line(io, "surface tension", system.surface_tension)
         summary_line(io, "surface normal method", system.surface_normal_method)
@@ -233,6 +237,9 @@ end
     return ndims(system) * factor_tvf(system) + 2
 end
 
+# To find the correct index of the density in the integration array
+@inline density_index(system::EntropicallyDampedSPHSystem) = v_nvariables(system) - 1
+
 system_correction(system::EntropicallyDampedSPHSystem) = system.correction
 
 @inline function particle_density(v, ::ContinuityDensity,

diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
@@ -206,8 +206,7 @@ end
     distance < sqrt(eps(typeof(distance))) && return
 
     (; delta) = density_diffusion
-    (; state_equation) = particle_system
-    (; sound_speed) = state_equation
+    sound_speed = system_sound_speed(particle_system)
 
     volume_b = m_b / rho_b
 
@@ -217,8 +216,9 @@ end
 
     smoothing_length_avg = (smoothing_length(particle_system, particle) +
                             smoothing_length(particle_system, neighbor)) / 2
-    dv[end, particle] += delta * smoothing_length_avg * sound_speed *
-                         density_diffusion_term
+
+    dv[density_index(particle_system), particle] += delta * smoothing_length_avg *
+                                                    sound_speed * density_diffusion_term
 end
 
 # Density diffusion `nothing` or interaction other than fluid-fluid

diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -126,14 +126,15 @@ end
                                                   v_particle_system, v_neighbor_system,
                                                   particle, neighbor, pos_diff, distance,
                                                   m_b, rho_a, rho_b,
-                                                  particle_system::WeaklyCompressibleSPHSystem,
+                                                  particle_system::FluidSystem,
                                                   neighbor_system, grad_kernel)
     (; density_diffusion) = particle_system
 
     vdiff = current_velocity(v_particle_system, particle_system, particle) -
             current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
-    dv[end, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel)
+    dv[density_index(particle_system), particle] += rho_a / rho_b *
+                                                    m_b * dot(vdiff, grad_kernel)
 
     # Artificial density diffusion should only be applied to system(s) representing a fluid
     # with the same physical properties i.e. density and viscosity.

diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -237,6 +237,9 @@ end
     return ndims(system) * factor_tvf(system) + 1
 end
 
+# To find the correct index of the density in the integration array
+@inline density_index(system::WeaklyCompressibleSPHSystem) = v_nvariables(system)
+
 system_correction(system::WeaklyCompressibleSPHSystem) = system.correction
 
 @propagate_inbounds function particle_pressure(v, system::WeaklyCompressibleSPHSystem,

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
@@ -2,16 +2,13 @@ using TrixiParticles
 
 # ==========================================================================================
 # ==== Resolution
-particle_spacings = [0.02, 0.01, 0.005]
+particle_spacings = [0.01]#, 0.005]
 
 # ==========================================================================================
 # ==== Experiment Setup
 tspan = (0.0, 5.0)
 reynolds_number = 100.0
 
-density_calculators = [ContinuityDensity(), SummationDensity()]
-perturb_coordinates = [false, true]
-
 function compute_l1v_error(v, u, t, system)
     v_analytical_avg = 0.0
     v_avg = 0.0
@@ -63,8 +60,8 @@ function diff_p_loc_p_avg(v, u, t, system)
     return v[end, :] .- p_avg_tot
 end
 
-for density_calculator in density_calculators, perturbation in perturb_coordinates,
-    particle_spacing in particle_spacings, wcsph in [false, true]
+for perturbation in [false], particle_spacing in particle_spacings,
+    density_calculator in [ContinuityDensity()], wcsph in [false]
 
     n_particles_xy = round(Int, 1.0 / particle_spacing)
 
@@ -94,3 +91,51 @@ for density_calculator in density_calculators, perturbation in perturb_coordinat
                   particle_spacing=particle_spacing, reynolds_number=reynolds_number,
                   tspan=tspan, saving_callback=saving_callback, pp_callback=pp_callback)
 end
+
+L_1_v_wcsph_conti = [Vector{Float64}() for _ in eachindex(particle_spacings)]
+L_1_v_wcsph_sum = [Vector{Float64}() for _ in eachindex(particle_spacings)]
+L_1_v_edac_conti = [Vector{Float64}() for _ in eachindex(particle_spacings)]
+L_1_v_edac_sum = [Vector{Float64}() for _ in eachindex(particle_spacings)]
+for perturbation in [false], (i, particle_spacing) in enumerate(particle_spacings[1:2]),
+    density_calculator in [ContinuityDensity(), SummationDensity()], wcsph in [false, true]
+
+    n_particles_xy = round(Int, 1.0 / particle_spacing)
+
+    name_density_calculator = density_calculator isa SummationDensity ?
+                              "summation_density" : "continuity_density"
+    name_perturbation = perturbation ? "_perturbed" : ""
+
+    input_directory = joinpath("out_tgv", name_density_calculator * name_perturbation,
+                               wcsph ? "wcsph" : "edac",
+                               "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)")
+    data = TrixiParticles.CSV.read(joinpath(input_directory, "errors.csv"),
+                                   TrixiParticles.DataFrame)
+
+    if wcsph && density_calculator isa ContinuityDensity
+        append!(L_1_v_wcsph_conti[i], copy(Vector(data.var"L1v_fluid_1")))
+    elseif wcsph && density_calculator isa SummationDensity
+        append!(L_1_v_wcsph_sum[i], copy(Vector(data.var"L1v_fluid_1")))
+    elseif !wcsph && density_calculator isa ContinuityDensity
+        append!(L_1_v_edac_conti[i], copy(Vector(data.var"L1v_fluid_1")))
+    elseif !wcsph && density_calculator isa SummationDensity
+        append!(L_1_v_edac_sum[i], copy(Vector(data.var"L1v_fluid_1")))
+    end
+end
+
+using Plots
+
+line_colors = cgrad(:coolwarm, length(particle_spacings), categorical=true)
+
+labels1 = "Δx" .* ["0.02" "0.01"] .* " ContinuityDensity"
+labels2 = "Δx" .* ["0.02" "0.01"] .* " SummationDensity"
+
+p1 = plot(L_1_v_wcsph_conti, label=labels1, linewidth=2, title="WCSPH",
+          palette=line_colors.colors, xlabel="Time", ylabel="L_1", legend=:topright)
+plot!(p1, L_1_v_wcsph_sum, label=labels2, linestyle=:dash, linewidth=2,
+      palette=line_colors.colors, xlabel="t/sec", ylabel="L_1", legend=:topright)
+
+p2 = plot(L_1_v_edac_conti, label=labels1, linewidth=2, title="EDAC",
+          palette=line_colors.colors,
+          xlabel="t/sec", ylabel="L_1", legend=:topright)
+plot!(p2, L_1_v_edac_sum, label=labels2, linestyle=:dash, linewidth=2,
+      palette=line_colors.colors, xlabel="t/sec", ylabel="L_1", legend=:topright)