Merge branch 'main' into more-bc

Project.toml

@@ -6,9 +6,11 @@ version = "0.1.0"
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Morton = "2a6d852e-3fac-5a38-885c-fe708af2d09e"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -21,8 +23,10 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 DiffEqCallbacks = "2.25"
+FastPow = "0.1"
 ForwardDiff = "0.10"
 Morton = "0.1"
+MuladdMacro = "0.2"
 Polyester = "0.7.5"
 Reexport = "1"
 SciMLBase = "1, 2"

src/TrixiParticles.jl

@@ -4,8 +4,10 @@ using Reexport: @reexport
 
 using Dates
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
+using FastPow: @fastpow
 using LinearAlgebra: norm, dot, I, tr
 using Morton: cartesian2morton
+using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,

src/general/density_calculators.jl

@@ -59,17 +59,15 @@ function summation_density!(system, system_index, semi, u, u_ode, density;
     # Use all other systems for the density summation
     @trixi_timeit timer() "compute density" foreach_enumerate(systems) do (neighbor_system_index,
                                                                            neighbor_system)
-        u_neighbor_system = wrap_u(u_ode, neighbor_system_index,
-                                   neighbor_system, semi)
+        u_neighbor_system = wrap_u(u_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,
+        for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords,
                               neighborhood_search,
                               particles=particles) do particle, neighbor, pos_diff, distance
             mass = hydrodynamic_mass(neighbor_system, neighbor)

src/general/initial_condition.jl

@@ -66,24 +66,28 @@ struct InitialCondition{ELTYPE}
 
     function InitialCondition(coordinates, velocities, masses, densities; pressure=0.0,
                               particle_spacing=-1.0)
+        ELTYPE = eltype(coordinates)
+
         if size(coordinates) != size(velocities)
             throw(ArgumentError("`coordinates` and `velocities` must be of the same size"))
         end
 
         if !(size(coordinates, 2) == length(masses) == length(densities))
-            throw(ArgumentError("the following must hold: " *
-                                "`size(coordinates, 2) == length(masses) == length(densities)`"))
+            throw(ArgumentError("Expected: size(coordinates, 2) == length(masses) == length(densities)\n" *
+                                "Got: size(coordinates, 2) = $(size(coordinates, 2)), " *
+                                "length(masses) = $(length(masses)), " *
+                                "length(densities) = $(length(densities))"))
         end
 
         if pressure isa Number
-            pressure = pressure * ones(length(masses))
+            pressure = pressure * ones(ELTYPE, length(masses))
         elseif length(pressure) != length(masses)
             throw(ArgumentError("`pressure` must either be a scalar or a vector of the " *
                                 "same length as `masses`"))
         end
 
-        return new{eltype(coordinates)}(particle_spacing, coordinates, velocities, masses,
-                                        densities, pressure)
+        return new{ELTYPE}(particle_spacing, coordinates, velocities, masses,
+                           densities, pressure)
     end
 end
 

src/general/smoothing_kernels.jl

@@ -163,24 +163,29 @@ For an analytic formula for higher order kernels, see (Monaghan, 1985).
 """
 struct SchoenbergQuarticSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end
 
-function kernel(kernel::SchoenbergQuarticSplineKernel, r::Real, h)
+# Note that `floating_point_number^integer_literal` is lowered to `Base.literal_pow`.
+# Currently, specializations reducing this to simple multiplications exist only up
+# to a power of three, see
+# https://github.com/JuliaLang/julia/blob/34934736fa4dcb30697ac1b23d11d5ad394d6a4d/base/intfuncs.jl#L327-L339
+# Higher powers just call a generic power implementation. Here, we accept to lose
+# some precision but gain performance by using plain multiplications instead via
+# `@fastpow`.
+@muladd @fastpow @inline function kernel(kernel::SchoenbergQuarticSplineKernel, r::Real, h)
     q = r / h
 
-    if q >= 5 / 2
-        return 0.0
-    end
+    q5_2_pow4 = (5 / 2 - q)^4
+    q3_2_pow4 = (3 / 2 - q)^4
+    q1_2_pow4 = (1 / 2 - q)^4
 
-    result = (5 / 2 - q)^4
+    # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl.
+    result = q5_2_pow4
+    result = result - 5 * (q < 3 / 2) * q3_2_pow4
+    result = result + 10 * (q < 1 / 2) * q1_2_pow4
 
-    if q < 3 / 2
-        result -= 5 * (3 / 2 - q)^4
+    # Zero out result if q >= 5/2
+    result = ifelse(q < 5 / 2, normalization_factor(kernel, h) * result, zero(result))
 
-        if q < 1 / 2
-            result += 10 * (1 / 2 - q)^4
-        end
-    end
-
-    return normalization_factor(kernel, h) * result
+    return result
 end
 
 function kernel_deriv(kernel::SchoenbergQuarticSplineKernel, r::Real, h)

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -101,6 +101,7 @@ function create_cache(n_particles, ELTYPE, ::SummationDensity)
 end
 
 function create_cache(n_particles, ELTYPE, ::ContinuityDensity)
+    # Density in this case is added to the end of 'v' and allocated by modifying 'v_nvariables'.
     return (;)
 end
 

test/general/density_calculator.jl

@@ -0,0 +1,36 @@
+using OrdinaryDiffEq
+
+# Setup a single particle and calculate its density
+@testset verbose=true "DensityCalculators" begin
+    @testset verbose=true "SummationDensity" begin
+        water_density = 1000.0
+
+        initial_condition = InitialCondition(zeros(2, 1), zeros(2, 1), [water_density],
+                                             [water_density])
+
+        smoothing_length = 1.0
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+
+        state_equation = StateEquationCole(10, 7, water_density, 100000.0,
+                                           background_pressure=100000.0)
+        viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
+
+        fluid_system = WeaklyCompressibleSPHSystem(initial_condition, SummationDensity(),
+                                                   state_equation,
+                                                   smoothing_kernel, smoothing_length,
+                                                   viscosity=viscosity)
+
+        (; cache) = fluid_system
+        (; density) = cache # Density is in the cache for SummationDensity
+
+        semi = Semidiscretization(fluid_system, neighborhood_search=GridNeighborhoodSearch,
+                                  damping_coefficient=1e-5)
+
+        tspan = (0.0, 0.01)
+        ode = semidiscretize(semi, tspan)
+        TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0)
+
+        @test density[1] ===
+              water_density * TrixiParticles.kernel(smoothing_kernel, 0.0, smoothing_length)
+    end
+end

test/general/general.jl

@@ -1,4 +1,5 @@
 # `util.jl` is added here to avoid confusion with `test_util.jl`
 include("util.jl")
 include("initial_condition.jl")
+include("density_calculator.jl")
 include("semidiscretization.jl")