From cca739d6287a4f40abe340242bec438655a339bc Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 12 Aug 2024 11:43:41 +0200 Subject: [PATCH] Cleanup examples (#592) * extend hydrostatic_water_column tests * revert * remove surface tension example with dam break * simplify accelerated_tank example * fix * fix * fix * fix accelerated_tank * reduce run time * run more tests for EDAC * update * fix edac tests * format * combine them * cleanup --- examples/fluid/accelerated_tank_2d.jl | 50 +---- .../fluid/dam_break_2d_surface_tension.jl | 31 --- examples/fluid/hydrostatic_water_column_2d.jl | 12 +- .../fluid/hydrostatic_water_column_edac_2d.jl | 35 --- src/general/semidiscretization.jl | 32 ++- test/examples/examples.jl | 204 ++++++++++++------ 6 files changed, 171 insertions(+), 193 deletions(-) delete mode 100644 examples/fluid/dam_break_2d_surface_tension.jl delete mode 100644 examples/fluid/hydrostatic_water_column_edac_2d.jl diff --git a/examples/fluid/accelerated_tank_2d.jl b/examples/fluid/accelerated_tank_2d.jl index a179a035c..efaf91cd4 100644 --- a/examples/fluid/accelerated_tank_2d.jl +++ b/examples/fluid/accelerated_tank_2d.jl @@ -4,30 +4,9 @@ using TrixiParticles using OrdinaryDiffEq -# ========================================================================================== -# ==== Resolution +# Resolution fluid_particle_spacing = 0.05 -# Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 3 - -# ========================================================================================== -# ==== Experiment Setup -tspan = (0.0, 1.0) - -# Boundary geometry and initial fluid particle positions -initial_fluid_size = (1.0, 0.9) -tank_size = (1.0, 1.0) - -fluid_density = 1000.0 -sound_speed = 10.0 -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7, clip_negative_pressure=false) - -tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, - acceleration=(0.0, -9.81), state_equation=state_equation, - n_layers=boundary_layers, spacing_ratio=1.0) - # Function for moving boundaries movement_function(t) = SVector(0.0, 0.5 * 9.81 * t^2) @@ -35,32 +14,7 @@ is_moving(t) = true boundary_movement = BoundaryMovement(movement_function, is_moving) -# Import the setup from `hydrostatic_water_column_2d.jl` trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), fluid_particle_spacing=fluid_particle_spacing, movement=boundary_movement, - acceleration=(0.0, 0.0), tank=tank, semi=nothing, ode=nothing, - sol=nothing) # Overwrite `sol` assignment to skip time integration - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system) -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=50) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") - -callbacks = CallbackSet(info_callback, saving_callback) - -# Use a Runge-Kutta method with automatic (error based) time step size control. -# Limiting of the maximum stepsize is necessary to prevent crashing. -# When particles are approaching a wall in a uniform way, they can be advanced -# with large time steps. Close to the wall, the stepsize has to be reduced drastically. -# Sometimes, the method fails to do so because forces become extremely large when -# fluid particles are very close to boundary particles, and the time integration method -# interprets this as an instability. -sol = solve(ode, RDPK3SpFSAL35(), - abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing - save_everystep=false, callback=callbacks); + tspan=(0.0, 1.0), system_acceleration=(0.0, 0.0)) diff --git a/examples/fluid/dam_break_2d_surface_tension.jl b/examples/fluid/dam_break_2d_surface_tension.jl deleted file mode 100644 index d83eb9865..000000000 --- a/examples/fluid/dam_break_2d_surface_tension.jl +++ /dev/null @@ -1,31 +0,0 @@ -# This example shows how surface tension can be applied to existing cases -# and which parameters have to be changed! -using TrixiParticles - -fluid_density = 1000.0 - -H = 0.6 -fluid_particle_spacing = H / 60 - -# Set the surface tension to a value that is accurate in your case. -# Note: This usually requires calibration to be physically accurate! -surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.025) - -# `density_diffusion` is deactivated since the interaction with the surface tension model can -# cause stability problems. -# `adhesion_coefficient` needs to be set to a value so that the fluid doesn't separate -# from the boundary. -# Note: The surface tension model leads to an increase in compressibility of the fluid, -# which needs to be rectified by an increase of the `sound_speed`. -# Note: The Wendland Kernels don't work very well here since the `SurfaceTensionAkinci` -# model is optimized for a compact support of `2 * particle_spacing`, which would result -# in a too small `smoothing_length` for the Wendland Kernel functions. -# Note: Adhesion will result in additional friction at the boundary. -trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), - surface_tension=surface_tension, - fluid_particle_spacing=fluid_particle_spacing, - smoothing_kernel=SchoenbergCubicSplineKernel{2}(), - smoothing_length=1.0 * fluid_particle_spacing, - correction=AkinciFreeSurfaceCorrection(fluid_density), - density_diffusion=nothing, adhesion_coefficient=0.05, - sound_speed=100.0, tspan=(0.0, 2.0)) diff --git a/examples/fluid/hydrostatic_water_column_2d.jl b/examples/fluid/hydrostatic_water_column_2d.jl index b35b6933f..86d392c5b 100644 --- a/examples/fluid/hydrostatic_water_column_2d.jl +++ b/examples/fluid/hydrostatic_water_column_2d.jl @@ -23,21 +23,25 @@ state_equation = StateEquationCole(; sound_speed, reference_density=fluid_densit exponent=7, clip_negative_pressure=false) tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, - n_layers=boundary_layers, - acceleration=(0.0, -gravity), state_equation=state_equation) + n_layers=boundary_layers, acceleration=(0.0, -gravity), + state_equation=state_equation) # ========================================================================================== # ==== Fluid smoothing_length = 1.2 * fluid_particle_spacing smoothing_kernel = SchoenbergCubicSplineKernel{2}() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +alpha = 0.02 +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) fluid_density_calculator = ContinuityDensity() + +# This is to set acceleration with `trixi_include` +system_acceleration = (0.0, -gravity) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, - acceleration=(0.0, -gravity), + acceleration=system_acceleration, source_terms=nothing) # ========================================================================================== diff --git a/examples/fluid/hydrostatic_water_column_edac_2d.jl b/examples/fluid/hydrostatic_water_column_edac_2d.jl deleted file mode 100644 index 701a1ccdb..000000000 --- a/examples/fluid/hydrostatic_water_column_edac_2d.jl +++ /dev/null @@ -1,35 +0,0 @@ -using TrixiParticles - -# Import the setup from `hydrostatic_water_column_2d.jl` -trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), - fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 0.9), - tank_size=(1.0, 1.0), state_equation=nothing, fluid_system=nothing, - semi=nothing, ode=nothing, sol=nothing) # Overwrite `sol` assignment to skip time integration - -# ========================================================================================== -# ==== Experiment Setup -tspan = (0.0, 1.0) - -# ========================================================================================== -# ==== Fluid -alpha = 0.02 -viscosity = ViscosityAdami(nu=alpha * smoothing_length * sound_speed / 8) - -fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, smoothing_length, - sound_speed, viscosity=viscosity, - density_calculator=ContinuityDensity(), - acceleration=(0.0, -gravity)) - -# ========================================================================================== -# ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system) -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=50) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") - -callbacks = CallbackSet(info_callback, saving_callback) - -# Use a Runge-Kutta method with automatic (error based) time step size control -sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks); diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 112131dc9..600d573d5 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -715,19 +715,27 @@ function update_nhs!(neighborhood_search, points_moving=(true, neighbor.ismoving[])) end +# This function is the same as the one below to avoid ambiguous dispatch when using `Union` function update_nhs!(neighborhood_search, - system::BoundarySPHSystem, - neighbor::Union{FluidSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem}, - u_system, u_neighbor) - # Don't update. This NHS is never used. - return neighborhood_search + system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, + neighbor::FluidSystem, u_system, u_neighbor) + # Depending on the density calculator of the boundary model, this NHS is used for + # - kernel summation (`SummationDensity`) + # - continuity equation (`ContinuityDensity`) + # - pressure extrapolation (`AdamiPressureExtrapolation`) + # + # Boundary coordinates only change over time when `neighbor.ismoving[]`. + # The current coordinates of fluids and solids change over time. + update!(neighborhood_search, system, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + points_moving=(system.ismoving[], true)) end +# This function is the same as the one above to avoid ambiguous dispatch when using `Union` function update_nhs!(neighborhood_search, system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, - neighbor::Union{FluidSystem, TotalLagrangianSPHSystem}, - u_system, u_neighbor) + neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor) # Depending on the density calculator of the boundary model, this NHS is used for # - kernel summation (`SummationDensity`) # - continuity equation (`ContinuityDensity`) @@ -773,6 +781,14 @@ function update_nhs!(neighborhood_search, points_moving=(true, false)) end +function update_nhs!(neighborhood_search, + system::BoundarySPHSystem, + neighbor::FluidSystem, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search +end + function update_nhs!(neighborhood_search, system::BoundaryDEMSystem, neighbor::Union{DEMSystem, BoundaryDEMSystem}, diff --git a/test/examples/examples.jl b/test/examples/examples.jl index dac2a2979..c962b5912 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -2,6 +2,96 @@ # but without checking the correctness of the solution. @testset verbose=true "Examples" begin @testset verbose=true "Fluid" begin + @trixi_testset "fluid/hydrostatic_water_column_2d.jl" begin + # Import variables into scope + trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "hydrostatic_water_column_2d.jl"), + fluid_system=nothing, sol=nothing, semi=nothing, ode=nothing) + + hydrostatic_water_column_tests = Dict( + "WCSPH default" => (), + "WCSPH with source term damping" => (source_terms=SourceTermDamping(damping_coefficient=1e-4),), + "WCSPH with SummationDensity" => (fluid_density_calculator=SummationDensity(), + clip_negative_pressure=true), + "WCSPH with ViscosityAdami" => ( + # from 0.02*10.0*1.2*0.05/8 + viscosity=ViscosityAdami(nu=0.0015),), + "WCSPH with ViscosityMorris" => ( + # from 0.02*10.0*1.2*0.05/8 + viscosity=ViscosityMorris(nu=0.0015),), + "WCSPH with ViscosityAdami and SummationDensity" => ( + # from 0.02*10.0*1.2*0.05/8 + viscosity=ViscosityAdami(nu=0.0015), + fluid_density_calculator=SummationDensity(), + clip_negative_pressure=true), + "WCSPH with ViscosityMorris and SummationDensity" => ( + # from 0.02*10.0*1.2*0.05/8 + viscosity=ViscosityMorris(nu=0.0015), + fluid_density_calculator=SummationDensity(), + clip_negative_pressure=true), + "WCSPH with smoothing_length=1.3" => (smoothing_length=1.3,), + "WCSPH with SchoenbergQuarticSplineKernel" => (smoothing_length=1.1, + smoothing_kernel=SchoenbergQuarticSplineKernel{2}()), + "WCSPH with SchoenbergQuinticSplineKernel" => (smoothing_length=1.1, + smoothing_kernel=SchoenbergQuinticSplineKernel{2}()), + "WCSPH with WendlandC2Kernel" => (smoothing_length=3.0, + smoothing_kernel=WendlandC2Kernel{2}()), + "WCSPH with WendlandC4Kernel" => (smoothing_length=3.5, + smoothing_kernel=WendlandC4Kernel{2}()), + "WCSPH with WendlandC6Kernel" => (smoothing_length=4.0, + smoothing_kernel=WendlandC6Kernel{2}()), + "EDAC with source term damping" => (source_terms=SourceTermDamping(damping_coefficient=1e-4), + fluid_system=EntropicallyDampedSPHSystem(tank.fluid, + smoothing_kernel, + smoothing_length, + sound_speed, + viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, + -gravity))), + "EDAC with SummationDensity" => (fluid_system=EntropicallyDampedSPHSystem(tank.fluid, + smoothing_kernel, + smoothing_length, + sound_speed, + viscosity=viscosity, + density_calculator=SummationDensity(), + acceleration=(0.0, + -gravity)),), + "EDAC with ViscosityAdami" => (fluid_system=EntropicallyDampedSPHSystem(tank.fluid, + smoothing_kernel, + smoothing_length, + sound_speed, + viscosity=ViscosityAdami(nu=0.0015), + density_calculator=ContinuityDensity(), + acceleration=(0.0, + -gravity)),), + "EDAC with ViscosityMorris" => (fluid_system=EntropicallyDampedSPHSystem(tank.fluid, + smoothing_kernel, + smoothing_length, + sound_speed, + viscosity=ViscosityMorris(nu=0.0015), + density_calculator=ContinuityDensity(), + acceleration=(0.0, + -gravity)),), + ) + + for (test_description, kwargs) in hydrostatic_water_column_tests + @testset "$test_description" begin + println("═"^100) + println("$test_description") + + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "hydrostatic_water_column_2d.jl"); + kwargs...) + + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + end + end + @trixi_testset "fluid/oscillating_drop_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", @@ -12,34 +102,6 @@ @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/hydrostatic_water_column_2d.jl" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "hydrostatic_water_column_2d.jl")) - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - - @trixi_testset "fluid/hydrostatic_water_column_2d.jl with source term damping" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "hydrostatic_water_column_2d.jl"), - source_terms=SourceTermDamping(; - damping_coefficient=1e-4)) - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - - @trixi_testset "fluid/hydrostatic_water_column_2d.jl with SummationDensity" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "hydrostatic_water_column_2d.jl"), - fluid_density_calculator=SummationDensity(), - clip_negative_pressure=true) - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - @trixi_testset "fluid/hydrostatic_water_column_3d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", @@ -60,14 +122,6 @@ @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/hydrostatic_water_column_edac_2d.jl" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "hydrostatic_water_column_edac_2d.jl")) - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - @trixi_testset "fluid/accelerated_tank_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), joinpath(examples_dir(), "fluid", @@ -77,34 +131,59 @@ end @trixi_testset "fluid/dam_break_2d.jl" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "dam_break_2d.jl"), tspan=(0.0, 0.1)) [ - r"┌ Info: The desired tank length in y-direction .*\n", - r"└ New tank length in y-direction.*\n", - ] - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 + # Import variables into scope + trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + boundary_layers=1, spacing_ratio=3, sol=nothing, semi=nothing, + ode=nothing) + + dam_break_tests = Dict( + "default" => (), + "with SummationDensity" => (fluid_density_calculator=SummationDensity(), + clip_negative_pressure=true), + "with DensityDiffusionMolteniColagrossi" => (density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),), + "no density diffusion" => (density_diffusion=nothing,), + "with KernelAbstractions" => (data_type=Array,), + "with BoundaryModelMonaghanKajtar" => (boundary_model=BoundaryModelMonaghanKajtar(gravity, + spacing_ratio, + boundary_particle_spacing, + tank.boundary.mass), + boundary_layers=1, spacing_ratio=3), + "with SurfaceTensionAkinci" => (surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.025), + fluid_particle_spacing=0.5 * + fluid_particle_spacing, + smoothing_kernel=SchoenbergCubicSplineKernel{2}(), + smoothing_length=0.5 * + fluid_particle_spacing, + correction=AkinciFreeSurfaceCorrection(fluid_density), + density_diffusion=nothing, + adhesion_coefficient=0.05, + sound_speed=100.0), + ) + + for (test_description, kwargs) in dam_break_tests + @testset "$test_description" begin + println("═"^100) + println("$test_description") + + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d.jl"); + tspan=(0, 0.1), kwargs...) [ + r"┌ Info: The desired tank length in y-direction .*\n", + r"└ New tank length in y-direction.*\n", + ] + + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + end end @trixi_testset "fluid/dam_break_oil_film_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_oil_film_2d.jl"), - tspan=(0.0, 0.1)) [ - r"┌ Info: The desired tank length in y-direction .*\n", - r"└ New tank length in y-direction.*\n", - ] - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - - @trixi_testset "fluid/dam_break_2d.jl with KernelAbstractions.jl" begin - # Emulate the GPU code on the CPU by passing `data_type = Array` - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "dam_break_2d.jl"), tspan=(0.0, 0.1), - data_type=Array) [ + tspan=(0.0, 0.05)) [ r"┌ Info: The desired tank length in y-direction .*\n", r"└ New tank length in y-direction.*\n", ] @@ -147,15 +226,6 @@ @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/dam_break_2d_surface_tension.jl" begin - @test_nowarn_mod trixi_include(@__MODULE__, - joinpath(examples_dir(), "fluid", - "dam_break_2d_surface_tension.jl"), - tspan=(0.0, 0.1)) - @test sol.retcode == ReturnCode.Success - @test count_rhs_allocations(sol, semi) == 0 - end - @trixi_testset "fluid/sphere_surface_tension_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid",