Change API of DensityDiffusionAntuono

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)
 

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,

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,

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,

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,

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,

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...)
 

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

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))

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
 

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