diff --git a/NEWS.md b/NEWS.md index 14a12e0635..fe05cfcf50 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### API Changes +- `DensityDiffusionAntuono` now only takes only one kwarg `delta` (#1142). - Return type of `vtk2trixi` changed to `NamedTuple` including an optional `:initial_condition` field if `create_initial_condition=true` is passed. (#959) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 56303b1877..1f36fac5f1 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -63,7 +63,7 @@ viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) # free surface in long-running simulations, but is significantly faster than the model # by Antuono. This simulation is short enough to use the faster model. density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) -# density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) +# density_diffusion = DensityDiffusionAntuono(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 778d404d69..1e952316d7 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -57,7 +57,7 @@ fluid_density_calculator = ContinuityDensity() nu = 0.005 alpha = 8 * nu / (fluid_smoothing_length * sound_speed) viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) -density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) +density_diffusion = DensityDiffusionAntuono(delta=0.1) sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, fluid_smoothing_length, diff --git a/examples/fluid/oscillating_drop_2d.jl b/examples/fluid/oscillating_drop_2d.jl index fb96050f42..ebc169ce23 100644 --- a/examples/fluid/oscillating_drop_2d.jl +++ b/examples/fluid/oscillating_drop_2d.jl @@ -55,7 +55,7 @@ smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) -density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1) +density_diffusion = DensityDiffusionAntuono(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index b2dad5bca6..7b4d14be7e 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -63,7 +63,7 @@ smoothing_kernel = WendlandC2Kernel{2}() state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=1, clip_negative_pressure=false) -density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1) +density_diffusion = DensityDiffusionAntuono(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(), state_equation, smoothing_kernel, smoothing_length, density_diffusion=density_diffusion, diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl index 22a132fd3e..35a7c4fe76 100644 --- a/examples/fsi/hydrostatic_water_column_2d.jl +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -103,7 +103,7 @@ if use_edac else fluid_density_calculator = ContinuityDensity() density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) - # density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) + # density_diffusion = DensityDiffusionAntuono(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length_fluid, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 75d3ba6074..f1c1fe7caf 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -178,14 +178,8 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit # as it was already verified in `allocate_buffer` that the density array is constant. density_rest = first(initial_condition.density) - if density_diffusion isa DensityDiffusionAntuono - density_diffusion_ = DensityDiffusionAntuono(initial_condition; - delta=density_diffusion.delta) - else - density_diffusion_ = density_diffusion - end - - cache = (; density_diffusion=density_diffusion_, + cache = (; density_diffusion=density_diffusion, + create_cache_density_diffusion(initial_condition, density_diffusion)..., pressure_boundary=pressure_boundary, density_rest=density_rest, cache...) diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index 4fabd5317a..10aee9d66c 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -88,7 +88,7 @@ end end @doc raw""" - DensityDiffusionAntuono(initial_condition; delta) + DensityDiffusionAntuono(; delta) The commonly used density diffusion terms by [Antuono (2010)](@cite Antuono2010), also referred to as δ-SPH. The density diffusion term by [Molteni (2009)](@cite Molteni2009) is extended by a second @@ -115,33 +115,31 @@ where ``d`` is the number of dimensions. See [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) for an overview and comparison of implemented density diffusion terms. """ -struct DensityDiffusionAntuono{NDIMS, ELTYPE, ARRAY2D, ARRAY3D} <: AbstractDensityDiffusion - delta :: ELTYPE - correction_matrix :: ARRAY3D # Array{ELTYPE, 3}: [i, j, particle] - normalized_density_gradient :: ARRAY2D # Array{ELTYPE, 2}: [i, particle] - - function DensityDiffusionAntuono(delta, correction_matrix, normalized_density_gradient) - new{size(correction_matrix, 1), typeof(delta), - typeof(normalized_density_gradient), - typeof(correction_matrix)}(delta, correction_matrix, - normalized_density_gradient) +struct DensityDiffusionAntuono{ELTYPE} <: AbstractDensityDiffusion + delta::ELTYPE + + function DensityDiffusionAntuono(; delta) + new{typeof(delta)}(delta) end end -function DensityDiffusionAntuono(initial_condition; delta) +create_cache_density_diffusion(initial_condition, density_diffusion) = (;) + +function create_cache_density_diffusion(initial_condition, + ::DensityDiffusionAntuono) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) - correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, - nparticles(initial_condition)) + n_particles = nparticles(initial_condition) - normalized_density_gradient = Array{ELTYPE, 2}(undef, NDIMS, - nparticles(initial_condition)) + density_diffusion_correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, + n_particles) + density_diffusion_normalized_density_gradient = Array{ELTYPE, 2}(undef, NDIMS, + n_particles) - return DensityDiffusionAntuono(delta, correction_matrix, normalized_density_gradient) + return (; density_diffusion_correction_matrix, + density_diffusion_normalized_density_gradient) end -@inline Base.ndims(::DensityDiffusionAntuono{NDIMS}) where {NDIMS} = NDIMS - function Base.show(io::IO, density_diffusion::DensityDiffusionAntuono) @nospecialize density_diffusion # reduce precompilation time @@ -154,22 +152,19 @@ function allocate_buffer(initial_condition, density_diffusion, buffer) return allocate_buffer(initial_condition, buffer), density_diffusion end -function allocate_buffer(ic, dd::DensityDiffusionAntuono, buffer::SystemBuffer) - initial_condition = allocate_buffer(ic, buffer) - return initial_condition, DensityDiffusionAntuono(initial_condition; delta=dd.delta) -end - @propagate_inbounds function density_diffusion_psi(density_diffusion::DensityDiffusionAntuono, rho_a, rho_b, pos_diff, distance, system, particle, neighbor) - (; normalized_density_gradient) = density_diffusion + (; density_diffusion_normalized_density_gradient) = system.cache # First term by Molteni & Colagrossi result = 2 * (rho_a - rho_b) # Second correction term - normalized_gradient_a = extract_svector(normalized_density_gradient, system, particle) - normalized_gradient_b = extract_svector(normalized_density_gradient, system, neighbor) + normalized_gradient_a = extract_svector(density_diffusion_normalized_density_gradient, + system, particle) + normalized_gradient_b = extract_svector(density_diffusion_normalized_density_gradient, + system, neighbor) result -= dot(normalized_gradient_a + normalized_gradient_b, pos_diff) # Since this is one of the most performance critical functions, using fast divisions @@ -179,13 +174,14 @@ end end function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi) - (; normalized_density_gradient) = density_diffusion + density_diffusion_correction_matrix = system.cache.density_diffusion_correction_matrix + normalized_density_gradient = system.cache.density_diffusion_normalized_density_gradient # Compute correction matrix density_fun = @propagate_inbounds(particle->current_density(v, system, particle)) system_coords = current_coordinates(u, system) - compute_gradient_correction_matrix!(density_diffusion.correction_matrix, + compute_gradient_correction_matrix!(density_diffusion_correction_matrix, system, system_coords, density_fun, semi) # Compute normalized density gradient @@ -198,7 +194,8 @@ function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi) rho_b = @inbounds current_density(v, system, neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - L = @inbounds correction_matrix(density_diffusion, particle) + L = @inbounds extract_smatrix(density_diffusion_correction_matrix, system, + particle) m_b = @inbounds hydrodynamic_mass(system, neighbor) volume_b = m_b / rho_b diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 7a9b290965..b00845b904 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -152,6 +152,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, stat n_particles)..., create_cache_refinement(initial_condition, particle_refinement, smoothing_length)..., + create_cache_density_diffusion(initial_condition, density_diffusion)..., create_cache_shifting(initial_condition, shifting_technique)..., # Per-system color tag for colorfield surface-normal logic and VTK output. color=Int(color_value)) diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl index 251c5465e8..2b8f361205 100644 --- a/test/io/read_vtk.jl +++ b/test/io/read_vtk.jl @@ -106,13 +106,13 @@ end @testset verbose=true "Exact Field Matches" begin - trixi2vtk(expected_ic; filename="tmp_initial_condition_exact_field_match", + trixi2vtk(expected_data; filename="tmp_initial_condition_exact_field_match", output_directory=tmp_dir, - center_of_mass_velocity=fill(42.0, size(expected_ic.velocity))) + center_of_mass_velocity=fill(42.0, size(expected_data.velocity))) file = joinpath(tmp_dir, "tmp_initial_condition_exact_field_match.vtu") test_ic = vtk2trixi(file) - @test isapprox(expected_ic.velocity, test_ic.velocity, rtol=1e-5) + @test isapprox(expected_data.velocity, test_ic.velocity, rtol=1e-5) end end diff --git a/test/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/test/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index 9e2c82a0f6..7dafc0377e 100644 --- a/test/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/test/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -2,13 +2,7 @@ @testset verbose=true "DensityDiffusionAntuono" begin # Use `@trixi_testset` to isolate the mock functions in a separate namespace @trixi_testset "show" begin - # Mock initial condition - initial_condition = Val(:initial_condition) - Base.ndims(::Val{:initial_condition}) = 2 - Base.eltype(::Val{:initial_condition}) = Float64 - TrixiParticles.nparticles(::Val{:initial_condition}) = 15 - - density_diffusion = DensityDiffusionAntuono(delta=0.1, initial_condition) + density_diffusion = DensityDiffusionAntuono(delta=0.1) @test repr(density_diffusion) == "DensityDiffusionAntuono(0.1)" end