diff --git a/Project.toml b/Project.toml index 85abbabaca..f00a02ee74 100644 --- a/Project.toml +++ b/Project.toml @@ -34,13 +34,13 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [weakdeps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" [extensions] -TrixiParticlesOrdinaryDiffEqExt = ["OrdinaryDiffEq", "OrdinaryDiffEqCore"] TrixiParticlesCUDAExt = "CUDA" +TrixiParticlesOrdinaryDiffEqExt = ["OrdinaryDiffEq", "OrdinaryDiffEqCore"] [compat] Accessors = "0.1.43" diff --git a/src/general/abstract_system.jl b/src/general/abstract_system.jl index fde3b92a7b..2b2b71f8f1 100644 --- a/src/general/abstract_system.jl +++ b/src/general/abstract_system.jl @@ -68,17 +68,18 @@ end end # Return `A[:, :, i]` as an `SMatrix`. -@inline function extract_smatrix(A, system, particle) - @boundscheck checkbounds(A, ndims(system), ndims(system), particle) +@propagate_inbounds function extract_smatrix(A, system, particle) + return extract_smatrix(A, Val(ndims(system)), particle) +end + +@inline function extract_smatrix(A, ::Val{NDIMS}, particle) where {NDIMS} + @boundscheck checkbounds(A, NDIMS, NDIMS, particle) # Extract the matrix elements for this particle as a tuple to pass to SMatrix - return SMatrix{ndims(system), - ndims(system)}(ntuple(@inline(i->@inbounds A[mod(i - 1, - ndims(system)) + 1, - div(i - 1, - ndims(system)) + 1, - particle]), - Val(ndims(system)^2))) + return SMatrix{NDIMS, NDIMS}(ntuple(@inline(i->@inbounds A[mod(i - 1, NDIMS) + 1, + div(i - 1, NDIMS) + 1, + particle]), + Val(NDIMS^2))) end # Specifically get the current coordinates of a particle for all system types. diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 163bec348e..b2953a6a20 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -11,6 +11,12 @@ function PointNeighbors.foreach_point_neighbor(f, system, neighbor_system, points, parallelization_backend) end +@propagate_inbounds function foreach_neighbor(f, system_coords, neighbor_coords, + neighborhood_search, backend, particle) + PointNeighbors.foreach_neighbor(f, system_coords, neighbor_coords, + neighborhood_search, particle) +end + # === Compact support selection === # -- Generic @inline function compact_support(system, neighbor) diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index a176be0327..3a546f4013 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -179,6 +179,11 @@ end return kernel(smoothing_kernel, distance, smoothing_length) end +@inline function smoothing_kernel_unsafe(system::WallBoundarySystem, distance, particle) + (; smoothing_kernel, smoothing_length) = system.boundary_model + return kernel_unsafe(smoothing_kernel, distance, smoothing_length) +end + @inline function smoothing_length(system::WallBoundarySystem, particle) return smoothing_length(system.boundary_model, particle) end diff --git a/src/schemes/structure/total_lagrangian_sph/penalty_force.jl b/src/schemes/structure/total_lagrangian_sph/penalty_force.jl index 025f1ab660..5dad18257e 100644 --- a/src/schemes/structure/total_lagrangian_sph/penalty_force.jl +++ b/src/schemes/structure/total_lagrangian_sph/penalty_force.jl @@ -14,38 +14,53 @@ struct PenaltyForceGanzenmueller{ELTYPE} end end -@inline function dv_penalty_force(penalty_force::Nothing, - particle, neighbor, initial_pos_diff, initial_distance, - current_pos_diff, current_distance, - system, m_a, m_b, rho_a, rho_b) - return zero(initial_pos_diff) +@inline function dv_penalty_force!(dv_particle, penalty_force::Nothing, + particle, neighbor, initial_pos_diff, initial_distance, + current_pos_diff, current_distance, + system, m_a, m_b, rho_a, rho_b, F_a, F_b) + return dv_particle end -@propagate_inbounds function dv_penalty_force(penalty_force::PenaltyForceGanzenmueller, - particle, neighbor, initial_pos_diff, - initial_distance, - current_pos_diff, current_distance, - system, m_a, m_b, rho_a, rho_b) - volume_a = m_a / rho_a - volume_b = m_b / rho_b +@propagate_inbounds function dv_penalty_force!(dv_particle, + penalty_force::PenaltyForceGanzenmueller, + particle, neighbor, initial_pos_diff, + initial_distance, + current_pos_diff, current_distance, + system, m_a, m_b, rho_a, rho_b, F_a, F_b) + (; alpha) = penalty_force - kernel_weight = smoothing_kernel(system, initial_distance, particle) + # Since this is one of the most performance critical functions, using fast divisions + # here gives a significant speedup on GPUs. + # See the docs page "Development" for more details on `div_fast`. + volume_a = div_fast(m_a, rho_a) + volume_b = div_fast(m_b, rho_b) - F_a = deformation_gradient(system, particle) - F_b = deformation_gradient(system, neighbor) + # This function is called after a compact support check, so we can use the unsafe + # kernel function, which does not check the distance again. + kernel_weight = smoothing_kernel_unsafe(system, initial_distance, particle) - inv_current_distance = 1 / current_distance + E_a = young_modulus(system, particle) + E_b = young_modulus(system, neighbor) - # Use the symmetry of epsilon to simplify computations - eps_sum = (F_a + F_b) * initial_pos_diff - 2 * current_pos_diff - delta_sum = dot(eps_sum, current_pos_diff) * inv_current_distance + # Note that this is actually ϵ^b_ab = -ϵ^b_ba in the paper, so we later compute + # δ^b_ab instead of δ^b_ba, but δ^b_ab = δ^b_ba because of antisymmetry of x_ab and ϵ_ab. + eps_a = F_a * initial_pos_diff - current_pos_diff + eps_b = F_b * initial_pos_diff - current_pos_diff - E = young_modulus(system, particle) + # This is (E_a * delta_a + E_b * delta_b) * current_distance. + # Pulling the division by `current_distance` out allows us to do one division by + # `current_distance^2` instead. + delta_sum = E_a * dot(eps_a, current_pos_diff) + E_b * dot(eps_b, current_pos_diff) - f = (penalty_force.alpha / 2) * volume_a * volume_b * - kernel_weight / initial_distance^2 * E * delta_sum * current_pos_diff * - inv_current_distance + # The division contains all scalar factors, which are then multiplied by + # the vector `current_pos_diff` at the end. + # We already divide by `m_a` to obtain an acceleration. + # Since this is one of the most performance critical functions, using fast divisions + # here gives a significant speedup on GPUs. + # See the docs page "Development" for more details on `div_fast`. + dv_particle[] += div_fast((alpha / 2) * volume_a * volume_b * kernel_weight * delta_sum, + initial_distance^2 * current_distance^2 * m_a) * + current_pos_diff - # Divide force by mass to obtain acceleration - return f / m_a + return dv_particle end diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index 874f1c92f7..fa000ef38f 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -19,6 +19,8 @@ end # Everything here is done in the initial coordinates system_coords = initial_coordinates(system) + neighborhood_search = get_neighborhood_search(system, semi) + backend = semi.parallelization_backend # For `distance == 0`, the analytical gradient is zero, but the unsafe gradient # and the density diffusion divide by zero. @@ -29,48 +31,61 @@ end h = initial_smoothing_length(system) almostzero = sqrt(eps(h^2)) - # Loop over all pairs of particles and neighbors within the kernel cutoff. - # For structure-structure interaction, this has to happen in the initial coordinates. - foreach_point_neighbor(system, system, system_coords, system_coords, semi; - points=each_integrated_particle(system)) do particle, neighbor, - initial_pos_diff, - initial_distance - # Skip neighbors with the same position because the kernel gradient is zero. - # Note that `return` only exits the closure, i.e., skips the current neighbor. - skip_zero_distance(system) && initial_distance < almostzero && return - - # Now that we know that `distance` is not zero, we can safely call the unsafe - # version of the kernel gradient to avoid redundant zero checks. - grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, - initial_distance, particle) - - rho_a = @inbounds system.material_density[particle] - rho_b = @inbounds system.material_density[neighbor] - + @threaded semi for particle in each_integrated_particle(system) + # We are looping over the particles of `system`, so it is guaranteed + # that `particle` is in bounds of `system`. m_a = @inbounds system.mass[particle] - m_b = @inbounds system.mass[neighbor] - + rho_a = @inbounds system.material_density[particle] # PK1 / rho^2 pk1_rho2_a = @inbounds pk1_rho2(system, particle) - pk1_rho2_b = @inbounds pk1_rho2(system, neighbor) - - current_pos_diff_ = @inbounds current_coords(system, particle) - - current_coords(system, neighbor) - # On GPUs, convert `Float64` coordinates to `Float32` after computing the difference - current_pos_diff = convert.(eltype(system), current_pos_diff_) - current_distance = norm(current_pos_diff) - - dv_stress = m_b * (pk1_rho2_a + pk1_rho2_b) * grad_kernel - - dv_penalty_force_ = @inbounds dv_penalty_force(penalty_force, particle, neighbor, - initial_pos_diff, initial_distance, - current_pos_diff, current_distance, - system, m_a, m_b, rho_a, rho_b) - - dv_particle = Ref(dv_stress + dv_penalty_force_) - @inbounds dv_viscosity_tlsph!(dv_particle, system, v_system, particle, neighbor, - current_pos_diff, current_distance, - m_a, m_b, rho_a, rho_b, grad_kernel) + current_coords_a = @inbounds current_coords(system, particle) + F_a = @inbounds deformation_gradient(system, particle) + + # Accumulate the RHS contributions over all neighbors before writing to `dv` + # to reduce the number of memory writes. + # Note that we need a `Ref` in order to be able to update these variables + # inside the closure in the `foreach_neighbor` loop. + dv_particle = Ref(zero(current_coords_a)) + + # Loop over all neighbors within the kernel cutoff + @inbounds foreach_neighbor(system_coords, system_coords, + neighborhood_search, backend, + particle) do particle, neighbor, + initial_pos_diff, initial_distance + # Skip neighbors with the same position because the kernel gradient is zero. + # Note that `return` only exits the closure, i.e., skips the current neighbor. + skip_zero_distance(system) && initial_distance < almostzero && return + + # Now that we know that `distance` is not zero, we can safely call the unsafe + # version of the kernel gradient to avoid redundant zero checks. + grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, + initial_distance, particle) + + rho_b = @inbounds system.material_density[neighbor] + m_b = @inbounds system.mass[neighbor] + # PK1 / rho^2 + pk1_rho2_b = @inbounds pk1_rho2(system, neighbor) + current_coords_b = @inbounds current_coords(system, neighbor) + + # The compiler is smart enough to optimize this away if no penalty force is used + F_b = @inbounds deformation_gradient(system, neighbor) + + current_pos_diff_ = current_coords_a - current_coords_b + # On GPUs, convert `Float64` coordinates to `Float32` after computing the difference + current_pos_diff = convert.(eltype(system), current_pos_diff_) + current_distance = norm(current_pos_diff) + + dv_particle[] += m_b * (pk1_rho2_a + pk1_rho2_b) * grad_kernel + + @inbounds dv_penalty_force!(dv_particle, penalty_force, particle, neighbor, + initial_pos_diff, initial_distance, + current_pos_diff, current_distance, + system, m_a, m_b, rho_a, rho_b, F_a, F_b) + + @inbounds dv_viscosity_tlsph!(dv_particle, system, v_system, particle, neighbor, + current_pos_diff, current_distance, + m_a, m_b, rho_a, rho_b, F_a, grad_kernel) + end for i in 1:ndims(system) @inbounds dv[i, particle] += dv_particle[][i] diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 81237caf14..dcb777d84e 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -489,30 +489,49 @@ end # Loop over all pairs of particles and neighbors within the kernel cutoff initial_coords = initial_coordinates(system) - foreach_point_neighbor(system, system, initial_coords, initial_coords, - semi) do particle, neighbor, initial_pos_diff, initial_distance - # Skip neighbors with the same position because the kernel gradient is zero. - # Note that `return` only exits the closure, i.e., skips the current neighbor. - skip_zero_distance(system) && initial_distance < almostzero && return + neighborhood_search = get_neighborhood_search(system, semi) + backend = semi.parallelization_backend - # Now that we know that `distance` is not zero, we can safely call the unsafe - # version of the kernel gradient to avoid redundant zero checks. - grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, - initial_distance, particle) - - volume = @inbounds mass[neighbor] / material_density[neighbor] - pos_diff_ = @inbounds current_coords(system, particle) - - current_coords(system, neighbor) - # On GPUs, convert `Float64` coordinates to `Float32` after computing the difference - pos_diff = convert.(eltype(system), pos_diff_) - - # Multiply by L_{0a} - L = @inbounds correction_matrix(system, particle) - - result = volume * pos_diff * grad_kernel' * L' + # The deformation gradient is computed for all particles, including the clamped ones + @threaded semi for particle in eachparticle(system) + # We are looping over the particles of `system`, so it is guaranteed + # that `particle` is in bounds of `system`. + current_coords_a = @inbounds current_coords(system, particle) + L_a = @inbounds correction_matrix(system, particle) + + # Accumulate the contributions over all neighbors before writing + # to `deformation_grad` to reduce the number of memory writes. + # Note that we need a `Ref` in order to be able to update these variables + # inside the closure in the `foreach_neighbor` loop. + result = Ref(zero(L_a)) + + # Loop over all neighbors within the kernel cutoff + @inbounds foreach_neighbor(initial_coords, initial_coords, + neighborhood_search, backend, + particle) do particle, neighbor, + initial_pos_diff, initial_distance + # Skip neighbors with the same position because the kernel gradient is zero. + # Note that `return` only exits the closure, i.e., skips the current neighbor. + skip_zero_distance(system) && initial_distance < almostzero && return + + # Now that we know that `distance` is not zero, we can safely call the unsafe + # version of the kernel gradient to avoid redundant zero checks. + grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff, + initial_distance, particle) + + volume = @inbounds mass[neighbor] / material_density[neighbor] + current_coords_b = @inbounds current_coords(system, neighbor) + pos_diff_ = current_coords_a - current_coords_b + # On GPUs, convert `Float64` coordinates to `Float32` after computing the difference + pos_diff = convert.(eltype(system), pos_diff_) + + # The tensor product pos_diff ⊗ (L_{0a} * ∇W) is equivalent to multiplication + # by the transpose: pos_diff * (L_{0a} * ∇W)ᵀ = pos_diff * ∇Wᵀ * L_{0a}ᵀ. + result[] -= volume * pos_diff * grad_kernel' * L_a' + end for j in 1:ndims(system), i in 1:ndims(system) - @inbounds deformation_grad[i, j, particle] -= result[i, j] + @inbounds deformation_grad[i, j, particle] += result[][i, j] end end diff --git a/src/schemes/structure/total_lagrangian_sph/viscosity.jl b/src/schemes/structure/total_lagrangian_sph/viscosity.jl index b19aec0117..618bfa7293 100644 --- a/src/schemes/structure/total_lagrangian_sph/viscosity.jl +++ b/src/schemes/structure/total_lagrangian_sph/viscosity.jl @@ -4,28 +4,28 @@ @propagate_inbounds function dv_viscosity_tlsph!(dv_particle, system, v_system, particle, neighbor, current_pos_diff, current_distance, - m_a, m_b, rho_a, rho_b, grad_kernel) + m_a, m_b, rho_a, rho_b, F_a, grad_kernel) viscosity = system.viscosity return dv_viscosity_tlsph!(dv_particle, viscosity, system, v_system, particle, neighbor, current_pos_diff, current_distance, - m_a, m_b, rho_a, rho_b, grad_kernel) + m_a, m_b, rho_a, rho_b, F_a, grad_kernel) end @propagate_inbounds function dv_viscosity_tlsph!(dv_particle, viscosity, system, v_system, particle, neighbor, current_pos_diff, current_distance, - m_a, m_b, rho_a, rho_b, grad_kernel) + m_a, m_b, rho_a, rho_b, F_a, grad_kernel) return viscosity(dv_particle, system, v_system, particle, neighbor, current_pos_diff, current_distance, - m_a, m_b, rho_a, rho_b, grad_kernel) + m_a, m_b, rho_a, rho_b, F_a, grad_kernel) end @inline function dv_viscosity_tlsph!(dv_particle, viscosity::Nothing, system, v_system, particle, neighbor, current_pos_diff, current_distance, - m_a, m_b, rho_a, rho_b, grad_kernel) - return zero(current_pos_diff) + m_a, m_b, rho_a, rho_b, F_a, grad_kernel) + return dv_particle end # Applying the viscosity according to Lin et al. (2015): @@ -38,7 +38,8 @@ end current_pos_diff, current_distance, m_a, m_b, rho_a, - rho_b, grad_kernel) + rho_b, F_a, + grad_kernel) v_a = current_velocity(v_system, system, particle) v_b = current_velocity(v_system, system, neighbor) v_diff = v_a - v_b @@ -54,10 +55,14 @@ end # Compute bulk modulus from Young's modulus and Poisson's ratio. # See the table at the end of https://en.wikipedia.org/wiki/Lam%C3%A9_parameters E = young_modulus(system, particle) + # A fast division is slower here for some reason K = E / (ndims(system) * (1 - 2 * poisson_ratio(system, particle))) # Newton–Laplace equation - sound_speed = sqrt(K / rho_a) + # Since this is one of the most performance critical functions, using fast divisions + # here gives a significant speedup on GPUs. + # See the docs page "Development" for more details on `div_fast`. + sound_speed = sqrt(div_fast(K, rho_a)) h_a = smoothing_length(system, particle) h_b = smoothing_length(system, neighbor) @@ -65,18 +70,20 @@ end rho_mean = (rho_a + rho_b) / 2 + # Since this is one of the most performance critical functions, using fast divisions + # here gives a significant speedup on GPUs. + # See the docs page "Development" for more details on `div_fast`. (; alpha, beta, epsilon) = viscosity - mu = h * vr / (current_distance^2 + epsilon * h^2) + mu = div_fast(h * vr, (current_distance^2 + epsilon * h^2)) c = sound_speed - pi_ab = (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel + pi_ab = div_fast(alpha * c * mu + beta * mu^2, rho_mean) * grad_kernel - F = deformation_gradient(system, particle) - det_F = det(F) + det_F = det(F_a) if abs(det_F) < 1.0f-9 return dv_particle end # See eq. 26 of Lin et al. (2015) - dv_particle[] += m_b * det_F * inv(F)' * pi_ab + dv_particle[] += m_b * det_F * inv(F_a)' * pi_ab end return dv_particle diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 2351ec7024..ae5add731a 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -21,6 +21,31 @@ end end + @trixi_testset "structure/oscillating_beam_2d.jl with penalty force" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "structure", + "oscillating_beam_2d.jl"), + tspan=(0.0, 0.1), + penalty_force=PenaltyForceGanzenmueller(alpha=1.0), + sol=nothing) [ + r"\[ Info: To create the self-interaction neighborhood search.*\n" + ] + # Use a fixed time step, tuned to the maximum stable step size for this example. + # Together with the very large penalty force alpha, this test will crash with + # "Instability detected" if the penalty force is not working correctly. + callbacks = CallbackSet((@__MODULE__).callbacks, StepsizeCallback(cfl=1.6)) + sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=1.0, + save_everystep=false, callback=callbacks) + @test sol.retcode == ReturnCode.Success + if VERSION < v"1.12" + # Older Julia versions produce allocations because `get_neighborhood_search` + # is not type-stable with TLSPH. + @test count_rhs_allocations(sol, semi) < 200 + else + @test count_rhs_allocations(sol, semi) == 0 + end + end + @trixi_testset "structure/oscillating_beam_2d.jl with penalty force and viscosity" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "structure", diff --git a/test/schemes/structure/total_lagrangian_sph/rhs.jl b/test/schemes/structure/total_lagrangian_sph/rhs.jl index 2b145859df..ce498b4db2 100644 --- a/test/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/test/schemes/structure/total_lagrangian_sph/rhs.jl @@ -58,37 +58,23 @@ kernel_deriv = 1.0 #### Mocking - struct MockSystem <: TrixiParticles.AbstractSystem{2} end + struct MockSystem <: TrixiParticles.AbstractSystem{2} + current_coordinates::Any + material_density::Any + pk1_rho2::Any + mass::Any + penalty_force::Any + viscosity::Any + buffer::Any + end @inline Base.eltype(::MockSystem) = Float64 - system = MockSystem() + system = MockSystem(current_coordinates, material_density, pk1_rho2, mass, + nothing, nothing, nothing) function TrixiParticles.initial_coordinates(::MockSystem) return initial_coordinates end - # Unpack calls should return predefined values or - # another mock object of the type Val{:mock_property_name}. - function Base.getproperty(::MockSystem, f::Symbol) - if f === :current_coordinates - return current_coordinates - elseif f === :material_density - return material_density - elseif f === :pk1_rho2 - return pk1_rho2 - elseif f === :mass - return mass - elseif f === :penalty_force - return nothing - elseif f === :viscosity - return nothing - elseif f === :buffer - return nothing - end - - # For all other properties, return mock objects - return Val(Symbol("mock_" * string(f))) - end - TrixiParticles.eachparticle(::MockSystem) = eachparticle TrixiParticles.each_integrated_particle(::MockSystem) = each_integrated_particle TrixiParticles.smoothing_length(::MockSystem, _) = eps() @@ -102,6 +88,9 @@ function TrixiParticles.current_coords(system::MockSystem, particle) return TrixiParticles.current_coords(initial_coordinates, system, particle) end + function TrixiParticles.deformation_gradient(::MockSystem, particle) + return nothing + end #### Verification backends = [ diff --git a/test/systems/tlsph_system.jl b/test/systems/tlsph_system.jl index 3bc3ea2ceb..8cb261f2d2 100644 --- a/test/systems/tlsph_system.jl +++ b/test/systems/tlsph_system.jl @@ -140,8 +140,8 @@ mass = [2.0, 2.0] density = 4.0 - correction_matrix = [1 0; 0 1;;; - 1 0; 0 1] + correction_matrix = [1.0 0.0; 0.0 1.0;;; + 1.0 0.0; 0.0 1.0] # This will cause the computed gradient # to be equal to `initial_coords[particle] - initial_coords[neighbor]` @@ -149,42 +149,32 @@ #### Mocking # Mock the system - system = Val(:mock_system_tensor) - TrixiParticles.ndims(::Val{:mock_system_tensor}) = 2 - Base.ntuple(f, ::Symbol) = ntuple(f, 2) # Make `extract_svector` work - function TrixiParticles.current_coords(system::Val{:mock_system_tensor}, - particle) + struct MockSystem <: TrixiParticles.AbstractStructureSystem{2} + mass::Any + correction_matrix::Any + material_density::Any + smoothing_length::Any + end + Base.eltype(::MockSystem) = Float64 + system = MockSystem(mass, correction_matrix, Val(:mock_material_density), + 0.12) + + function TrixiParticles.current_coords(::MockSystem, particle) return TrixiParticles.extract_svector(current_coordinates[i], Val(2), particle) end - function TrixiParticles.initial_coordinates(::Val{:mock_system_tensor}) + function TrixiParticles.initial_coordinates(::MockSystem) return initial_coordinates[i] end - TrixiParticles.smoothing_length(::Val{:mock_system_tensor}, _) = 0.12 - - # All unpack calls should return another mock object - # of the type `Val{:mock_property_name}`, but we want to have some real matrices - # as properties as opposed to only mock objects. - function Base.getproperty(::Val{:mock_system_tensor}, f::Symbol) - if f === :correction_matrix - return correction_matrix - elseif f === :mass - return mass - end - - # For all other properties, return mock objects. - return Val(Symbol("mock_" * string(f))) - end - function TrixiParticles.compact_support(::Val{:mock_system_tensor}, - ::Val{:mock_system_tensor}) + function TrixiParticles.compact_support(::MockSystem, ::MockSystem) return Inf end Base.getindex(::Val{:mock_material_density}, ::Int64) = density - function TrixiParticles.smoothing_kernel_grad_unsafe(::Val{:mock_system_tensor}, + function TrixiParticles.smoothing_kernel_grad_unsafe(::MockSystem, pos_diff, distance, particle) return kernel_derivative * pos_diff / distance diff --git a/test/test_util.jl b/test/test_util.jl index 51b39f0b0d..12ead6e7a9 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -78,7 +78,13 @@ end @inline function TrixiParticles.get_neighborhood_search(system, semi::DummySemidiscretization) - return get_neighborhood_search(system, system, semi) + return TrixiParticles.get_neighborhood_search(system, system, semi) +end + +# Avoid method ambiguity +@inline function TrixiParticles.get_neighborhood_search(system::TotalLagrangianSPHSystem, + semi::DummySemidiscretization) + return TrixiParticles.get_neighborhood_search(system, system, semi) end include("count_allocations.jl")