From 231e2e020cf484ccb46fc28d0eee39c4490aef91 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Fri, 22 Nov 2024 15:13:20 +0100 Subject: [PATCH 1/5] Validation TVF (#640) * add `density_calculator` to example * add perturbation option * fix typo * implement suggestions --------- Co-authored-by: Sven Berger --- examples/fluid/taylor_green_vortex_2d.jl | 7 ++++++- .../validation_taylor_green_vortex_2d.jl | 15 +++++++++++++-- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 551863ccb..fc4ad9e63 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -58,13 +58,18 @@ background_pressure = sound_speed^2 * fluid_density smoothing_length = 1.0 * particle_spacing smoothing_kernel = SchoenbergQuinticSplineKernel{2}() +# To be set via `trixi_include` +perturb_coordinates = true fluid = RectangularShape(particle_spacing, (n_particles_xy, n_particles_xy), (0.0, 0.0), - coordinates_perturbation=0.2, # To avoid stagnant streamlines when not using TVF. + # Perturb particle coordinates to avoid stagnant streamlines without TVF + coordinates_perturbation=perturb_coordinates ? 0.2 : nothing, # To avoid stagnant streamlines when not using TVF. density=fluid_density, pressure=initial_pressure_function, velocity=initial_velocity_function) +density_calculator = SummationDensity() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, + density_calculator=density_calculator, transport_velocity=TransportVelocityAdami(background_pressure), viscosity=ViscosityAdami(; nu)) diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index cd9d18159..99b83ca03 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -9,6 +9,9 @@ particle_spacings = [0.02, 0.01, 0.005] tspan = (0.0, 5.0) reynolds_number = 100.0 +density_calculators = [SummationDensity(), ContinuityDensity()] +perturb_coordinates = [false, true] + function compute_l1v_error(v, u, t, system) v_analytical_avg = 0.0 v_avg = 0.0 @@ -60,10 +63,16 @@ function diff_p_loc_p_avg(v, u, t, system) return v[end, :] .- p_avg_tot end -for particle_spacing in particle_spacings +for density_calculator in density_calculators, perturbation in perturb_coordinates, + particle_spacing in particle_spacings + n_particles_xy = round(Int, 1.0 / particle_spacing) - output_directory = joinpath("out_tgv", + name_density_calculator = density_calculator isa SummationDensity ? + "summation_density" : "continuity_density" + name_perturbation = perturbation ? "perturbed" : "" + + output_directory = joinpath("out_tgv", "$(name_density_calculator)_$name_perturbation", "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)") saving_callback = SolutionSavingCallback(dt=0.02, output_directory=output_directory, @@ -79,6 +88,8 @@ for particle_spacing in particle_spacings # Import variables into scope trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "taylor_green_vortex_2d.jl"), + density_calculator=density_calculator, + perturb_coordinates=perturbation, particle_spacing=particle_spacing, reynolds_number=reynolds_number, tspan=tspan, saving_callback=saving_callback, pp_callback=pp_callback) end From 83e5c3bdeb34c07fa9cdd5b64dfaab7d6491aeb2 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 1 Dec 2024 07:36:23 +0100 Subject: [PATCH 2/5] Bump crate-ci/typos from 1.26.8 to 1.28.1 (#678) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.26.8 to 1.28.1. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.26.8...v1.28.1) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 27a81105f..3d5579157 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.26.8 + uses: crate-ci/typos@v1.28.1 From a1cc9d611b26e353ee100323a3ac95ebb7966627 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 2 Dec 2024 10:19:31 +0100 Subject: [PATCH 3/5] Bump codecov/codecov-action from 4 to 5 (#679) Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 4 to 5. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v4...v5) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 06c868258..cae43ee26 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -95,7 +95,7 @@ jobs: - name: Upload coverage report to Codecov # Only run coverage in one Job (Ubuntu and latest Julia version) if: matrix.os == 'ubuntu-latest' && matrix.version == '1' - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: files: lcov.info fail_ci_if_error: true From c93ab39970eacd67ca8e98a484db6dc814380473 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Mon, 2 Dec 2024 11:53:35 +0100 Subject: [PATCH 4/5] Plot recipes: support automatic marker size with `xlims` and `ylims` (#677) * add `xlims` and `ylmis` * implement suggestions --- src/visualization/recipes_plots.jl | 34 +++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index 53b588e34..a14de2948 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -10,7 +10,8 @@ RecipesBase.@recipe function f(sol::TrixiParticlesODESolution) end RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; - size=(600, 400)) # Default size + size=(600, 400), # Default size + xlims=(-Inf, Inf), ylims=(-Inf, Inf)) systems_data = map(semi.systems) do system u = wrap_u(u_ode, system, semi) coordinates = active_coordinates(u, system) @@ -25,6 +26,14 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; x_min, y_min = minimum(coordinates, dims=2) .- 0.5particle_spacing x_max, y_max = maximum(coordinates, dims=2) .+ 0.5particle_spacing + # `x_min`, `x_max`, etc. are used to automatically set the marker size. + # When `xlims` or `ylims` are passed explicitly, we have to update these to get the correct marker size. + isfinite(first(xlims)) && (x_min = xlims[1]) + isfinite(last(xlims)) && (x_max = xlims[2]) + + isfinite(first(ylims)) && (y_min = ylims[1]) + isfinite(last(ylims)) && (y_max = ylims[2]) + return (; x, y, x_min, x_max, y_min, y_max, particle_spacing, label=timer_name(system)) end @@ -32,7 +41,8 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; return (semi, systems_data...) end -RecipesBase.@recipe function f((initial_conditions::InitialCondition)...) +RecipesBase.@recipe function f((initial_conditions::InitialCondition)...; + xlims=(Inf, Inf), ylims=(Inf, Inf)) idx = 0 ics = map(initial_conditions) do ic x = collect(ic.coordinates[1, :]) @@ -46,6 +56,14 @@ RecipesBase.@recipe function f((initial_conditions::InitialCondition)...) x_min, y_min = minimum(ic.coordinates, dims=2) .- 0.5particle_spacing x_max, y_max = maximum(ic.coordinates, dims=2) .+ 0.5particle_spacing + # `x_min`, `x_max`, etc. are used to automatically set the marker size. + # When `xlims` or `ylims` are passed explicitly, we have to update these to get the correct marker size. + isfinite(first(xlims)) && (x_min = xlims[1]) + isfinite(last(xlims)) && (x_max = xlims[2]) + + isfinite(first(ylims)) && (y_min = ylims[1]) + isfinite(last(ylims)) && (y_max = ylims[2]) + idx += 1 return (; x, y, x_min, x_max, y_min, y_max, particle_spacing, @@ -56,12 +74,22 @@ RecipesBase.@recipe function f((initial_conditions::InitialCondition)...) end RecipesBase.@recipe function f(::Union{InitialCondition, Semidiscretization}, - data...; zcolor=nothing, size=(600, 400), colorbar_title="") + data...; zcolor=nothing, size=(600, 400), colorbar_title="", + xlims=(Inf, Inf), ylims=(Inf, Inf)) x_min = minimum(obj.x_min for obj in data) x_max = maximum(obj.x_max for obj in data) + y_min = minimum(obj.y_min for obj in data) y_max = maximum(obj.y_max for obj in data) + # `x_min`, `x_max`, etc. are used to automatically set the marker size. + # When `xlims` or `ylims` are passed explicitly, we have to update these to get the correct marker size. + isfinite(first(xlims)) && (x_min = xlims[1]) + isfinite(last(xlims)) && (x_max = xlims[2]) + + isfinite(first(ylims)) && (y_min = ylims[1]) + isfinite(last(ylims)) && (y_max = ylims[2]) + # Note that this assumes the plot area to be ~10% smaller than `size`, # which is the case when showing a single plot with the legend inside. # With the legend outside, this is no longer the case, so the `markersize` has to be From f051b55922aceab34ebe436ecf90a8a0b9d5e6d6 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 11 Dec 2024 06:39:53 +0100 Subject: [PATCH 5/5] Improve Test Coverage (#674) * refactoring and improve coverage * add test * cleanup * format * Update test/schemes/boundary/dummy_particles/dummy_particles.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update test/schemes/boundary/dummy_particles/dummy_particles.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * Update test/schemes/boundary/dummy_particles/dummy_particles.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> * fix --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- test/callbacks/postprocess.jl | 12 + .../dummy_particles/dummy_particles.jl | 273 +++++++++++------- 2 files changed, 173 insertions(+), 112 deletions(-) diff --git a/test/callbacks/postprocess.jl b/test/callbacks/postprocess.jl index 372a7b453..551392000 100644 --- a/test/callbacks/postprocess.jl +++ b/test/callbacks/postprocess.jl @@ -1,4 +1,16 @@ @testset verbose=true "PostprocessCallback" begin + @testset verbose=true "errors" begin + error_str1 = "`funcs` cannot be empty" + @test_throws ArgumentError(error_str1) PostprocessCallback(interval=10, + write_file_interval=0) + + error_str2 = "setting both `interval` and `dt` is not supported" + @test_throws ArgumentError(error_str2) PostprocessCallback(interval=10, + write_file_interval=0, + dt=0.1, + another_function=(v, u, t, system) -> 1) + end + @testset verbose=true "show" begin function example_function(v, u, t, system) return 0 diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index a72ac409c..dc8491cbb 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -1,142 +1,191 @@ @testset verbose=true "Dummy Particles" begin @testset "show" begin - boundary_model = BoundaryModelDummyParticles([1000.0], [1.0], + boundary_model = BoundaryModelDummyParticles([1000.0], + [1.0], SummationDensity(), - SchoenbergCubicSplineKernel{2}(), 0.1) + SchoenbergCubicSplineKernel{2}(), + 0.1) - show_compact = "BoundaryModelDummyParticles(SummationDensity, Nothing)" - @test repr(boundary_model) == show_compact + expected_repr = "BoundaryModelDummyParticles(SummationDensity, Nothing)" + @test repr(boundary_model) == expected_repr end - @testset "Viscosity Adami: Wall Velocity" begin + @testset "Viscosity Adami/Bernoulli: Wall Velocity" begin particle_spacing = 0.1 # Boundary particles in fluid compact support - boundary_1 = RectangularShape(particle_spacing, (10, 1), (0.0, 0.2), density=257.0) - boundary_2 = RectangularShape(particle_spacing, (10, 1), (0.0, 0.1), density=257.0) + boundary_1 = RectangularShape(particle_spacing, + (10, 1), + (0.0, 0.2), + density=257.0) + boundary_2 = RectangularShape(particle_spacing, + (10, 1), + (0.0, 0.1), + density=257.0) # Boundary particles out of fluid compact support - boundary_3 = RectangularShape(particle_spacing, (10, 1), (0, 0), density=257.0) + boundary_3 = RectangularShape(particle_spacing, + (10, 1), + (0.0, 0.0), + density=257.0) boundary = union(boundary_1, boundary_2, boundary_3) particles_in_compact_support = length(boundary_1.mass) + length(boundary_2.mass) - fluid = RectangularShape(particle_spacing, (16, 5), (-0.3, 0.3), density=257.0, + # Define fluid particles + fluid = RectangularShape(particle_spacing, + (16, 5), + (-0.3, 0.3), + density=257.0, loop_order=:x_first) + # Simulation parameters smoothing_kernel = SchoenbergCubicSplineKernel{2}() smoothing_length = 1.2 * particle_spacing viscosity = ViscosityAdami(nu=1e-6) - state_equation = StateEquationCole(sound_speed=10, reference_density=257, + state_equation = StateEquationCole(sound_speed=10.0, + reference_density=257.0, exponent=7) - boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass, - state_equation=state_equation, - AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length, - viscosity=viscosity) - - boundary_system = BoundarySPHSystem(boundary, boundary_model) - - fluid_system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), - state_equation, - smoothing_kernel, smoothing_length) - - neighborhood_search = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius=1.0, - eachpoint=TrixiParticles.eachparticle(fluid_system)) - - velocities = [[0; -1], [1; 1], [-1; 0], [0.7; 0.2], [0.3; 0.8]] - - @testset "Wall Velocity $v_fluid" for v_fluid in velocities - viscosity = boundary_system.boundary_model.viscosity - volume = boundary_system.boundary_model.cache.volume - - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - boundary_system.boundary_model.viscosity) - TrixiParticles.boundary_pressure_extrapolation!(boundary_model, - boundary_system, - fluid_system, - boundary.coordinates, - fluid.coordinates, v_fluid, - v_fluid .* - ones(size(fluid.coordinates)), - neighborhood_search) - - for particle in TrixiParticles.eachparticle(boundary_system) - if volume[particle] > eps() - TrixiParticles.compute_wall_velocity!(viscosity, boundary_system, - boundary.coordinates, particle) + # Define pressure extrapolation methods to test + pressure_extrapolations = [ + AdamiPressureExtrapolation(), + BernoulliPressureExtrapolation() + ] + + for pressure_extrapolation in pressure_extrapolations + @testset "Pressure Extrapolation: $(typeof(pressure_extrapolation))" begin + # Create boundary and fluid systems + boundary_model = BoundaryModelDummyParticles(boundary.density, + boundary.mass, + state_equation=state_equation, + pressure_extrapolation, + smoothing_kernel, + smoothing_length, + viscosity=viscosity) + boundary_system = BoundarySPHSystem(boundary, boundary_model) + fluid_system = WeaklyCompressibleSPHSystem(fluid, + SummationDensity(), + state_equation, + smoothing_kernel, + smoothing_length) + + neighborhood_search = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius=1.0, + eachpoint=TrixiParticles.eachparticle(fluid_system)) + + velocities = [ + [0.0; -1.0], + [1.0; 1.0], + [-1.0; 0.0], + [0.7; 0.2], + [0.3; 0.8] + ] + + @testset "Wall Velocity $v_fluid" for v_fluid in velocities + viscosity = boundary_system.boundary_model.viscosity + volume = boundary_system.boundary_model.cache.volume + + # Reset cache and perform pressure extrapolation + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + boundary_system.boundary_model.viscosity) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, + v_fluid, + v_fluid .* + ones(size(fluid.coordinates)), + neighborhood_search) + + # Compute wall velocities + for particle in TrixiParticles.eachparticle(boundary_system) + if volume[particle] > eps() + TrixiParticles.compute_wall_velocity!(viscosity, + boundary_system, + boundary.coordinates, + particle) + end + end + + # Expected wall velocities + v_wall = zeros(size(boundary.coordinates)) + v_wall[:, 1:particles_in_compact_support] .= -v_fluid + + @test isapprox(boundary_system.boundary_model.cache.wall_velocity, + v_wall) end - end - - v_wall = zeros(size(boundary.coordinates)) - v_wall[:, 1:particles_in_compact_support] .= -v_fluid - - @test isapprox(boundary_system.boundary_model.cache.wall_velocity, v_wall) - end - - scale_v = [1, 0.5, 0.7, 1.8, 67.5] - @testset "Wall Velocity Staggerd: Factor $scale" for scale in scale_v - viscosity = boundary_system.boundary_model.viscosity - volume = boundary_system.boundary_model.cache.volume - - # For a constant velocity profile (each fluid particle has the same velocity) - # the wall velocity is `v_wall = -v_fluid` (see eq. 22 in Adami_2012). - # Thus, generate a staggered velocity profile to test the smoothed velocity field - # for a variable velocity profile. - v_fluid = zeros(size(fluid.coordinates)) - for i in TrixiParticles.eachparticle(fluid_system) - if mod(i, 2) == 1 - v_fluid[:, i] .= scale - end - end - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - boundary_system.boundary_model.viscosity) - TrixiParticles.boundary_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, - boundary.coordinates, - fluid.coordinates, v_fluid, - v_fluid, - neighborhood_search) - - for particle in TrixiParticles.eachparticle(boundary_system) - if volume[particle] > eps() - TrixiParticles.compute_wall_velocity!(viscosity, boundary_system, - boundary.coordinates, particle) + scale_factors = [1.0, 0.5, 0.7, 1.8, 67.5] + + # For a constant velocity profile (each fluid particle has the same velocity), + # the wall velocity is `v_wall = -v_fluid` (see eq. 22 in Adami_2012). + # With a staggered velocity profile, we can test the smoothed velocity field + # for a variable velocity profile. + @testset "Wall Velocity Staggered: Factor $scale" for scale in scale_factors + viscosity = boundary_system.boundary_model.viscosity + volume = boundary_system.boundary_model.cache.volume + + # Create a staggered velocity profile + v_fluid = zeros(size(fluid.coordinates)) + for i in TrixiParticles.eachparticle(fluid_system) + if mod(i, 2) == 1 + v_fluid[:, i] .= scale + end + end + + # Reset cache and perform pressure extrapolation + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + boundary_system.boundary_model.viscosity) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, + v_fluid, + v_fluid, + neighborhood_search) + + # Compute wall velocities + for particle in TrixiParticles.eachparticle(boundary_system) + if volume[particle] > eps() + TrixiParticles.compute_wall_velocity!(viscosity, + boundary_system, + boundary.coordinates, + particle) + end + end + + # Expected wall velocities + v_wall = zeros(size(boundary.coordinates)) + + # First boundary row + for i in 1:length(boundary_1.mass) + if mod(i, 2) == 1 + # Particles with a diagonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.42040669416720744 * scale + else + # Particles with an orthogonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.5795933058327924 * scale + end + end + + # Second boundary row + for i in (length(boundary_1.mass) + 1):particles_in_compact_support + if mod(i, 2) == 1 + # Particles with a diagonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.12101100073462243 * scale + else + # Particles with an orthogonal distance to a fluid particle with v_fluid > 0.0 + v_wall[:, i] .= -0.8789889992653775 * scale + end + end + + @test isapprox(boundary_system.boundary_model.cache.wall_velocity, + v_wall) end end - - v_wall = zeros(size(boundary.coordinates)) - - # First boundary row - for i in 1:length(boundary_1.mass) - if mod(i, 2) == 1 - - # Particles with a diagonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.42040669416720744 * scale - else - - # Particles with a orthogonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.5795933058327924 * scale - end - end - - # Second boundary row - for i in (length(boundary_1.mass) + 1):particles_in_compact_support - if true == mod(i, 2) - - # Particles with a diagonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.12101100073462243 * scale - else - - # Particles with a orthogonal distance to a fluid particle with `v_fluid > 0.0` - v_wall[:, i] .= -0.8789889992653775 * scale - end - end - - @test isapprox(boundary_system.boundary_model.cache.wall_velocity, v_wall) end end end