diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index f69f5b78d1..7e4045a3e9 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -369,7 +369,7 @@ function InitialCondition(sol::ODESolution, system, semi; use_final_velocity=fal v_ode, u_ode = sol.u[end].x u = wrap_u(u_ode, system, semi) - v = wrap_u(v_ode, system, semi) + v = wrap_v(v_ode, system, semi) # Check if particles come too close especially when the surface exhibits large curvature too_close = find_too_close_particles(u, min_particle_distance) diff --git a/test/general/initial_condition.jl b/test/general/initial_condition.jl index 0c2949adef..39e0edceec 100644 --- a/test/general/initial_condition.jl +++ b/test/general/initial_condition.jl @@ -155,6 +155,48 @@ @test_throws ArgumentError(error_str) apply_angular_velocity(ic_1d, 2.0) end + @testset "InitialCondition from solution uses final velocity" begin + struct FinalVelocityMock{IC, M} <: TrixiParticles.AbstractStructureSystem{2} + initial_condition::IC + mass::M + end + + Base.eltype(::FinalVelocityMock) = Float64 + TrixiParticles.v_nvariables(::FinalVelocityMock) = 3 + TrixiParticles.compact_support(::FinalVelocityMock, neighbor) = 1.0 + function TrixiParticles.write_u0!(u0, system::FinalVelocityMock) + u0 .= system.initial_condition.coordinates + return u0 + end + function TrixiParticles.write_v0!(v0, system::FinalVelocityMock) + v0[1:2, :] .= system.initial_condition.velocity + v0[3, :] .= system.initial_condition.density + return v0 + end + + coordinates = [0.0 1.0 2.0 + 0.0 0.0 0.0] + initial_condition = InitialCondition(; coordinates, + velocity=zeros(2, 3), + density=ones(3), mass=ones(3)) + system = FinalVelocityMock(initial_condition, initial_condition.mass) + semi = Semidiscretization(system; neighborhood_search=nothing) + ode = semidiscretize(semi, (0.0, 1.0)) + + final_state = copy(ode.u0) + expected_velocity = [1.0 2.0 3.0 + 4.0 5.0 6.0] + v = TrixiParticles.wrap_v(final_state.x[1], system, semi) + v[1:2, :] .= expected_velocity + v[3, :] .= initial_condition.density + + sol = TrixiParticles.SciMLBase.build_solution(ode, nothing, [0.0, 1.0], + [ode.u0, final_state]) + restart_ic = InitialCondition(sol, system, semi; use_final_velocity=true) + + @test restart_ic.velocity == expected_velocity + end + @testset "Automatic Mass Calculation" begin particle_spacing = 0.13 coordinates = [88.3 10.4 5.2 48.3 58.9;