diff --git a/examples/fluid/calibration_water_drop_shape.jl b/examples/fluid/calibration_water_drop_shape.jl new file mode 100644 index 000000000..916033223 --- /dev/null +++ b/examples/fluid/calibration_water_drop_shape.jl @@ -0,0 +1,248 @@ +# In this example we try to approach the static shape of a water droplet on a horizontal plane. +# The shape of a static droplet can be calculated from the Young-Laplace equation. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0001 + +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 5.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (0.05, 0.01) + +fluid_density = 1000.0 +sound_speed = 50 + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.0015 + +sphere1_center = (0.5, sphere_radius) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere()) + +# ========================================================================================== +# ==== Fluid +# fluid_smoothing_length = 1.2 * fluid_particle_spacing +# fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_smoothing_length = 2.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_smoothing_length_2 = 2.5 * fluid_particle_spacing +fluid_smoothing_kernel_2 = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +# water at 20C +#nu=0.00089 +nu = 0.00089 +# too much 0.00089 -> 0.00045 -> 0.0002 -> 0.0001 +# no impact 0.00005 + +# the following term only holds for 2d sims +# alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ViscosityAdami(nu=nu) +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +# with increased smoothing length surface_tension is too small +sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.008), + correction=AkinciFreeSurfaceCorrection(fluid_density)) + +# 0.001 +# 0.2 > 90 +# 0.15 > 90 +# 0.125 > 90 +# 0.11 60-90 +# 0.1 ~ 45 + +# 0.0008 +# 0.11 ~45-50 +# 0.1125 ~90 +# 0.115 ~90 +# 0.12 > 90 + +# 0.0005 +#0.11 ~100 + +# 0.00025 +# 0.115 ~60 grad aber zu tief in der mitte 0.006 +# 0.115 and increase nu to 0.0002 + +# 0.0001 deg 90 (x-axis: 2mm, y-axis: 1.8mm) +# start x-axis ~1.7mm y-axis: 2.2mm ~110deg distance to boundary to large at adhesion coefficient 1.0 +# increase adhesion coefficient to 1.1 x-axis ~1.8mm y-axis: 2mm ~100deg distance still to high +# increase adhesion coefficient to 1.2 x-axis ~2.8 y-axis: 1.7 distance still to high +# decrease adhesion coefficient to 1.15 x-axis 2.4mm y-axis: 1.9 distance high +# decrease adhesion coefficient to 1.125 x-axis 2-2.2mm y-axis: 1.9 distance high ~90+ +# decrease surface tension coefficient from 0.115 to 0.11 x-axis: 2.6mm y-axis: 1.7 -> invalid too much adhesion +# increase surface tension coefficient from 0.11 to 0.1125 x-axis: 2-2.2mm y-axis: 1.9 +# increase surface tension coefficient from 0.1125 to 0.12 x-axis: 1.9-2mm y-axis: 2.0 +# increase viscosity from 0.0001 to 0.0002 x-axis: 1.8mm y-axis:2.2 +# increase adhesion coefficient from 1.125 to 1.15 x-axis: 1.8mm y-axis:2.1 +# increase adhesion coefficient from 1.15 to 1.2 x-axis: 1.8mm y-axis:2.1 +# increase adhesion coefficient from 1.2 to 1.25 x-axis: 2.6mm y-axis:1.7 +# increase viscosity from 0.0002 to 0.0003 x-axis: 2.1mm y-axis:1.9 <=== (adh=1.25, surft=0.12, nu=0.0003) +# increase adhesion coefficient from 1.25 to 1.275 x-axis: 2.6mm y-axis:1.7 +# decrease adhesion coefficient from 1.275 to 1.26 x-axis: 2.2mm y-axis:1.8 +# decrease adhesion coefficient from 1.26 to 1.255 x-axis: 1.8-2.6mm y-axis:1.8 +# increase viscosity from 0.0003 to 0.00031 x-axis: 2.2-2.8mm y-axis:1.7 +# increase viscosity from 0.00031 to 0.00032 x-axis: 2.2-2.8mm y-axis:1.7 + +# 0.0001 deg 75 (x-axis: 2.2mm, y-axis: 1.6mm) +#(adh=1.25, surft=0.12, nu=0.0003) +# decrease surft=0.11 x-axis: 2.4-2.6mm y-axis:1.6 +# decrease adh=1.2 x-axis: 2.2-2.4mm y-axis:1.8 +# increase viscosity = 0.00035 x-axis: 2.4-3.2mm y-axis:1.6 +# increase viscosity = 0.0004 x-axis: 2.6-2.8mm y-axis:1.4 +# increase viscosity = 0.00045 x-axis: 2.4-3.2mm y-axis:1.4 +# decrease adh=1.15 x-axis: 1.8mm y-axis:2.1 +# increase adh=1.175 x-axis: 1.8mm y-axis:2.2 +# increase adh=1.19 x-axis: 1.8mm y-axis:2.1 +#(adh=1.25, surft=0.12, nu=0.0003) +# increase adh=1.3 and x-axis: 2.4mm y-axis:1.8 +# decrease adh=1.275 and x-axis: 2.4mm y-axis:1.8 +# decrease adh=1.26 and x-axis: 2.4mm y-axis:1.8 +# decrease adh=1.255 and x-axis: 2.4mm y-axis:1.8 +# decrease smoothing length from 4 -> 3.75 x-axis: 2mm y-axis:1.9 + +# 60deg x=2.4, y=1.3 +# adh = 1.15, surft=0.115, nu=0.0003, h=3.75, x=1.7-1.8, y=2.2 +# adh = 1.2, surft=0.115, nu=0.0003, h=3.75, x=2-2.2, y=1.8 +# adh = 1.25, surft=0.115, nu=0.0003, h=3.75, x=2.2, y=1.8 +# adh = 1.3, surft=0.115, nu=0.0003, h=3.75, x=2.6, y=1.6 +# adh = 1.3, surft=0.12, nu=0.0003, h=3.75, x=2.2-2.6, y=1.6 +# adh = 1.3, surft=0.125, nu=0.0003, h=3.75, x=2, y=1.8 +# adh = 1.35, surft=0.125, nu=0.0003, h=3.75, x=2.5, y=1.6 +# adh = 1.35, surft=0.13, nu=0.0003, h=3.75, x=2-2.2, y=1.8 +# adh = 1.4, surft=0.13, nu=0.0003, h=3.75, x=2-2.2, y=1.8 +# adh = 1.45, surft=0.13, nu=0.0003, h=3.75, x=2.6, y=1.6 +# not possible with smoothing_length -> lower to 3 and switch from C4 to C2 +# adh = 1.45, surft=0.13, nu=0.0003, h=3.0, x=1.4, y=2.6 +# adh = 1.45, surft=0.12, nu=0.0003, h=3.0, x=1.4, y=2.4 +# adh = 1.5, surft=0.12, nu=0.0003, h=3.0, x=1.6, y=2.6 +# adh = 1.5, surft=0.11, nu=0.0003, h=3.0, x=1.4, y=2.4 +# adh = 1.5, surft=0.10, nu=0.0003, h=3.0, x=1.4, y=2.4 +# adh = 1.5, surft=0.05, nu=0.0003, h=3.0, x=3.5, y=2.0 +# adh = 1.4, surft=0.075, nu=0.0003, h=3.0, x=3, y=2.0 +# adh = 1.4, surft=0.075, nu=0.0003, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.05, nu=0.0003, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.01, nu=0.0003, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.01, nu=0.0001, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.001, nu=0.00089, h=2.0, x=3, y=2.0 crap +# adh = 1.2, surft=0.002, nu=0.00089, h=2.0, x=3, y=2.0 crap +# adh = 1.2, surft=0.002, nu=0.00089, h=2.5, x=2.5, y=1.4 close but too much adh +# adh = 0.6, surft=0.002, nu=0.00089, h=2.5, x=2.7, y=1.3 close but too much adh +# adh = 0.3, surft=0.002, nu=0.00089, h=2.5, x=2.8, y=1.2 +# adh = 0.15, surft=0.002, nu=0.00089, h=2.5, x=2.8, y=1.2 +# adh = 0.075, surft=0.002, nu=0.00089, h=2.5, x=3.4, y=1.2 +# adh = 0.075, surft=0.004, nu=0.00089, h=2.5, x=2.8, y=1.4 +# adh = 0.075, surft=0.01, nu=0.00089, h=2.5, x=2.2, y=1.7 +# adh = 0.075, surft=0.007, nu=0.00089, h=2.5, x=2.6, y=1.4 <- okish but form needs more adh +# adh = 0.075, surft=0.008, nu=0.00089, h=2.5, x=2.6, y=1.4 +# adh = 0.1, surft=0.01, nu=0.00089, h=2.5, x=2.4, y=1.6 <- okish but form needs more adh +# adh = 0.15, surft=0.01, nu=0.00089, wall_nu=0.00089, h=2.5, x=2.4, y=1.8 +# adh = 0.15, surft=0.01, nu=0.00089, h=2.5, wall_nu=0.5*0.00089, x=2.2, y=1.7 <-lower wall viscosity +# adh = 0.15, surft=0.01, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.6, y=1.6 <-lower wall viscosity +# adh = 0.15, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.2, y=1.8 +# adh = 0.2, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.2, y=1.8 +# adh = 0.2, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.6, y=1.7 +# adh = 0.25, surft=0.014, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.4, y=1.6 +# adh = 0.25, surft=0.014, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.8, y=1.5 diverged for up to 3s +# adh = 0.25, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 diverged at 2.722 <---- +# adh = 0.25, surft=0.016, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 +# adh = 0.1, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2.2, y=1.7 <- 75d +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2.2, y=1.7 +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.3, y=1.7 +# adh = 0.2, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.6, y=1.6 +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.00089, x=2.6, y=1.5 +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=5 * 0.00089, x=2.4, y=1.6 +# adh = 0.13, surft=0.011, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.6, y=1.5 +# adh = 0.11, surft=0.011, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.4, y=1.6 + +# 75deg (x-axis: 2.2mm, y-axis: 1.6mm) +# adh = 0.20, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2, y=1.7 <- 90d +# adh = 0.1, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2.2, y=1.7 +# adh = 0.1, surft=0.0125, nu=0.00089, h=2.5, wall_nu=7*0.00089, x=2.2, y=1.6 +# adh = 0.08, surft=0.011, nu=0.00089, h=2.5, wall_nu=7*0.00089, x=2.3, y=1.6 +# adh = 0.07, surft=0.011, nu=0.00089, h=2.5, wall_nu=7*0.00089, x=2.2, y=1.6 +# adh = 0.07, surft=0.011, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.3, y=1.6 +# adh = 0.05, surft=0.008, nu=0.00089, h=2.5, wall_nu=2*0.00089, x=2.2, y=1.6 + +# 90deg (x-axis: 2mm, y-axis: 1.8mm) +# adh = 0.15, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.2, y=1.7 after capallariy test +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.4*0.00089, x=2.2, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.00089, x=2.2, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=2*0.00089, x=2.2, y=1.7 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=4*0.00089, x=2.2, y=1.7 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2, y=1.7 <- does not work in capallariy test +# adh = 0.04, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.1, y=1.7 <- further testing showed surface tension is about 10-20% too high +# adh = 0.04, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 +# adh = 0.03, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 **** <- seems to be still too high +# adh = 0.02, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.8 +# adh = 0.02, surft=0.011, nu=0.00089, h=2.5, wall_nu=15*0.00089, x=2.2, y=1.7 +# adh = 0.01, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 +# adh = 0.01, surft=0.008, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.2, y=1.6 +# adh = 0.01, surft=0.008, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.6 + +# sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel_2, +# fluid_smoothing_length_2, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel_2, + fluid_smoothing_length_2, + viscosity=ViscosityAdami(nu=10 * nu)) + +# adhesion_coefficient = 1.0 and surface_tension_coefficient=0.01 for perfect wetting +# adhesion_coefficient = 0.001 and surface_tension_coefficient=2.0 for no wetting +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + adhesion_coefficient=0.01) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(boundary_system, sphere_surface_tension) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", + write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-4, maxiters=1E12, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/dam_break_2_phase_2d.jl b/examples/fluid/dam_break_2_phase_2d.jl new file mode 100644 index 000000000..8c1196023 --- /dev/null +++ b/examples/fluid/dam_break_2_phase_2d.jl @@ -0,0 +1,104 @@ +# 2D dam break simulation with an oil film on top + +using TrixiParticles +using OrdinaryDiffEq + +# Size parameters +H = 0.6 +W = 2 * H + +gravity = 9.81 +tspan = (0.0, 2.0) + +# Resolution +fluid_particle_spacing = H / 60 + +# Numerical settings +smoothing_length = 3.5 * fluid_particle_spacing +sound_speed = 100 + +# physical values +nu_water = 8.9E-7 +nu_air = 1.544E-5 +nu_ratio = nu_water / nu_air + +nu_sim_air = 0.02 * smoothing_length * sound_speed +nu_sim_water = nu_ratio * nu_sim_air + +air_viscosity = ViscosityMorris(nu=nu_sim_air) +water_viscosity = ViscosityMorris(nu=nu_sim_water) + +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + sol=nothing, fluid_particle_spacing=fluid_particle_spacing, + viscosity=water_viscosity, smoothing_length=smoothing_length, + gravity=gravity, tspan=tspan, + density_diffusion=nothing, + sound_speed=sound_speed, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.08), + cfl=0.8, + tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing, + 2.6 * H)) + +# TODO: broken (fixed?) +# trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), +# sol=nothing, fluid_particle_spacing=fluid_particle_spacing, +# viscosity=ViscosityMorris(nu=nu), smoothing_length=smoothing_length, +# gravity=gravity, tspan=tspan, +# density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1), +# sound_speed=sound_speed, +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), +# correction=AkinciFreeSurfaceCorrection(1000.0)) + +# ========================================================================================== +# ==== Setup oil layer + +oil_size = (tank_size[1], 1.5 * H) +oil_size2 = (tank_size[1] - W, H) +oil_density = 1.0 + +oil = RectangularShape(fluid_particle_spacing, + round.(Int, oil_size ./ fluid_particle_spacing), + zeros(length(oil_size)), density=oil_density) + +oil2 = RectangularShape(fluid_particle_spacing, + round.(Int, oil_size2 ./ fluid_particle_spacing), + (W, 0.0), density=oil_density) + +# move on top of the water +for i in axes(oil.coordinates, 2) + oil.coordinates[:, i] .+= [0.0, H] +end + +oil = union(oil, oil2) + +oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, + StateEquationCole(; sound_speed, + reference_density=oil_density, + exponent=1, + clip_negative_pressure=false), + smoothing_kernel, + smoothing_length, viscosity=air_viscosity, + #density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) + +# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, +# StateEquationIdealGas(; gas_constant=287.0, temperature=293.0, gamma=1.4), +# smoothing_kernel, +# smoothing_length, viscosity=air_viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, oil_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing)) +ode = semidiscretize(semi, tspan, data_type=nothing) + +# sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +# dt=1.0, # This is overwritten by the stepsize callback +# save_everystep=false, callback=callbacks); +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # 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); diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 117dde772..73c590b0d 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -35,7 +35,8 @@ tank_size = (floor(5.366 * H / boundary_particle_spacing) * boundary_particle_sp fluid_density = 1000.0 sound_speed = 20 * sqrt(gravity * H) state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=1, clip_negative_pressure=false) + exponent=1, clip_negative_pressure=false, + background_pressure=0.0) tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, n_layers=boundary_layers, spacing_ratio=spacing_ratio, @@ -59,9 +60,11 @@ density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, + reference_particle_spacing=fluid_particle_spacing, density_diffusion=density_diffusion, acceleration=(0.0, -gravity), correction=nothing, - surface_tension=nothing) + surface_tension=nothing, + reference_particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Boundary diff --git a/examples/fluid/dam_break_2phase_2d.jl b/examples/fluid/dam_break_2phase_2d.jl new file mode 100644 index 000000000..f6ca1c3c7 --- /dev/null +++ b/examples/fluid/dam_break_2phase_2d.jl @@ -0,0 +1,95 @@ +# 2D dam break simulation with an air layer on top + +using TrixiParticles +using OrdinaryDiffEq + +# Size parameters +H = 0.6 +W = 2 * H + +gravity = 9.81 +tspan = (0.0, 2.0) + +# Resolution +# Note: The resolution is very coarse. A better result is obtained with H/60 (which takes over 1 hour) +fluid_particle_spacing = H / 20 + +# Numerical settings +smoothing_length = 3.5 * fluid_particle_spacing +#sound_speed = 100 +# when using the Ideal gas equation +sound_speed = 400 + +# physical values +nu_water = 8.9E-7 +nu_air = 1.544E-5 +nu_ratio = nu_water / nu_air + +nu_sim_air = 0.02 * smoothing_length * sound_speed +nu_sim_water = nu_ratio * nu_sim_air + +air_viscosity = ViscosityMorris(nu=nu_sim_air) +water_viscosity = ViscosityMorris(nu=nu_sim_water) + +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + sol=nothing, fluid_particle_spacing=fluid_particle_spacing, + viscosity=water_viscosity, smoothing_length=smoothing_length, + gravity=gravity, tspan=tspan, density_diffusion=nothing, + sound_speed=sound_speed, cfl=0.8, exponent=7, background_pressure=0.0, + tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing, + 2.6 * H)) + +# ========================================================================================== +# ==== Setup air_system layer + +air_size = (tank_size[1], 1.5 * H) +air_size2 = (tank_size[1] - W, H) +air_density = 1.0 + +# Air above the initial water +air_system = RectangularShape(fluid_particle_spacing, + round.(Int, air_size ./ fluid_particle_spacing), + zeros(length(air_size)), density=air_density) + +# Air for the rest of the empty volume +air_system2 = RectangularShape(fluid_particle_spacing, + round.(Int, air_size2 ./ fluid_particle_spacing), + (W, 0.0), density=air_density) + +# move on top of the water +for i in axes(air_system.coordinates, 2) + air_system.coordinates[:, i] .+= [0.0, H] +end + +air_system = union(air_system, air_system2) + +# We use background_pressure here to prevent negative pressure regions in the gas phase. +air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator, + StateEquationCole(; sound_speed, + reference_density=air_density, + exponent=1, + clip_negative_pressure=false, + background_pressure=1000), + smoothing_kernel, smoothing_length, + viscosity=air_viscosity, + acceleration=(0.0, -gravity)) + +# air_system_system = WeaklyCompressibleSPHSystem(air_system, fluid_density_calculator, +# StateEquationIdealGas(; gamma=1.4, +# gas_constant=287.0, +# temperature=293.15), +# smoothing_kernel, smoothing_length, +# viscosity=air_viscosity, +# acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, air_system_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing)) +ode = semidiscretize(semi, tspan, data_type=nothing) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # 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); diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index ecfe66481..11ce6a54d 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -17,21 +17,39 @@ fluid_particle_spacing = H / 60 smoothing_length = 3.5 * fluid_particle_spacing sound_speed = 20 * sqrt(gravity * H) -nu = 0.02 * smoothing_length * sound_speed / 8 -oil_viscosity = ViscosityMorris(nu=10 * nu) +# Physical values +nu_water = 8.9E-7 +nu_oil = 6E-5 +nu_ratio = nu_water / nu_oil +nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) +nu_sim_water = nu_ratio * nu_sim_oil + +oil_viscosity = ViscosityMorris(nu=nu_sim_oil) +# Physical values +nu_water = 8.9E-7 +nu_oil = 6E-5 +nu_ratio = nu_water / nu_oil + +nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) +nu_sim_water = nu_ratio * nu_sim_oil + +oil_viscosity = ViscosityMorris(nu=nu_sim_oil) + +# TODO: broken if both are set to surface tension trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, - viscosity=ViscosityMorris(nu=nu), smoothing_length=smoothing_length, - gravity=gravity, - density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1), - sound_speed=sound_speed) + viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, + gravity=gravity, tspan=tspan, density_diffusion=nothing, + sound_speed=sound_speed, prefix="") # ========================================================================================== # ==== Setup oil layer oil_size = (W, 0.1 * H) oil_density = 700.0 +oil_eos = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1, + clip_negative_pressure=false) oil = RectangularShape(fluid_particle_spacing, round.(Int, oil_size ./ fluid_particle_spacing), @@ -43,12 +61,19 @@ for i in axes(oil.coordinates, 2) end oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, - state_equation, smoothing_kernel, + oil_eos, smoothing_kernel, smoothing_length, viscosity=oil_viscosity, - density_diffusion=density_diffusion, acceleration=(0.0, -gravity), - surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), - correction=AkinciFreeSurfaceCorrection(oil_density)) + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01), + correction=AkinciFreeSurfaceCorrection(oil_density), + reference_particle_spacing=fluid_particle_spacing) + +# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, +# oil_eos, smoothing_kernel, +# smoothing_length, viscosity=oil_viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03), +# reference_particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Simulation diff --git a/examples/fluid/drops_on_box.jl b/examples/fluid/drops_on_box.jl new file mode 100644 index 000000000..54974e4dd --- /dev/null +++ b/examples/fluid/drops_on_box.jl @@ -0,0 +1,143 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (1.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +box = RectangularShape(fluid_particle_spacing, (300, 125), (0.485, 0.0), + density=fluid_density) + +# box = SphereShape(fluid_particle_spacing, 1.0, (0.5, 0.0), +# fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0)) + +sphere_radius = 0.0025 + +sphere1_center = (0.5, 0.05) +sphere2_center = (0.5, 0.1) +sphere3_center = (0.5, 0.15) +sphere4_center = (0.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere3 = SphereShape(fluid_particle_spacing, sphere_radius, sphere3_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere4 = SphereShape(fluid_particle_spacing, sphere_radius, sphere4_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) + +water = union(sphere1, sphere2, sphere3, sphere4) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 1.0 * fluid_particle_spacing +fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.00089 +alpha = 0.75 * 8 * nu / (fluid_smoothing_length * sound_speed) +# viscosity = ViscosityAdami(nu=nu) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) + +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +sphere_surface_tension = WeaklyCompressibleSPHSystem(water, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.1 * + 0.011), + correction=AkinciFreeSurfaceCorrection(fluid_density), + reference_particle_spacing=fluid_particle_spacing) + +# sphere_surface_tension = WeaklyCompressibleSPHSystem(water, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity)) + +# sphere_surface_tension2 = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.011), +# correction=AkinciFreeSurfaceCorrection(fluid_density)) + +# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01), +# correction=AkinciFreeSurfaceCorrection(fluid_density)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=0.5 * nu)) + +boundary_model_box = BoundaryModelDummyParticles(box.density, box.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=ViscosityAdami(nu=0.5 * nu)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + adhesion_coefficient=0.01) + +boundary_system2 = BoundarySPHSystem(box, boundary_model_box, + adhesion_coefficient=0.01) + +# boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# boundary_system2 = BoundarySPHSystem(box, boundary_model_box) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(boundary_system, boundary_system2, sphere_surface_tension) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.0025, output_directory="out", + prefix="four_drop_terminal_m20_90d_01surft_05wallnu_075viscmon_001adh", + write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 + reltol=1e-3, # Default reltol is 1e-3 + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index cd41d6da8..18b19cfb2 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -40,7 +40,7 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, # ========================================================================================== # ==== Fluid -fluid_smoothing_length = 1.0 * fluid_particle_spacing +fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps() fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() fluid_density_calculator = ContinuityDensity() @@ -50,13 +50,13 @@ alpha = 8 * nu / (fluid_smoothing_length * sound_speed) viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) -sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, - state_equation, fluid_smoothing_kernel, +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, fluid_smoothing_length, - viscosity=viscosity, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity), surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), - correction=AkinciFreeSurfaceCorrection(fluid_density)) + reference_particle_spacing=fluid_particle_spacing) sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, state_equation, fluid_smoothing_kernel, @@ -75,16 +75,16 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar viscosity=ViscosityAdami(nu=wall_viscosity)) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, - adhesion_coefficient=0.001) + adhesion_coefficient=1.0) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(boundary_system, sphere_surface_tension, sphere) +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) -saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", - write_meta_data=true) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="", write_meta_data=true) callbacks = CallbackSet(info_callback, saving_callback) diff --git a/examples/fluid/falling_water_spheres_3d.jl b/examples/fluid/falling_water_spheres_3d.jl index a4ebb799e..078defee7 100644 --- a/examples/fluid/falling_water_spheres_3d.jl +++ b/examples/fluid/falling_water_spheres_3d.jl @@ -3,30 +3,30 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -fluid_particle_spacing = 0.008 +fluid_particle_spacing = 0.005 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -nu = 0.01 +nu = 0.001 fluid_density = 1000.0 sound_speed = 50 sphere1_radius = 0.05 -sphere1_center = (0.5, 0.5, 0.2) -sphere2_center = (1.5, 0.5, 0.2) +sphere1_center = (0.5, 0.5, 0.075) +sphere2_center = (1.5, 0.5, 0.075) sphere1 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere1_center, - fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0)) + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0)) sphere2 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere2_center, - fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0)) + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0)) # `compact_support` needs to be `2.0 * particle_spacing` to be correct fluid_smoothing_length = 1.0 * fluid_particle_spacing trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"), - fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.2), + fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.1), initial_fluid_size=(0.0, 0.0, 0.0), tank_size=(2.0, 1.0, 0.1), sound_speed=sound_speed, faces=(true, true, true, true, true, false), @@ -34,4 +34,4 @@ trixi_include(@__MODULE__, fluid_smoothing_length=fluid_smoothing_length, fluid_smoothing_kernel=SchoenbergCubicSplineKernel{3}(), nu=nu, alpha=10 * nu / (fluid_smoothing_length * sound_speed), - surface_tension_coefficient=1.5, adhesion_coefficient=0.5) + surface_tension_coefficient=10, adhesion_coefficient=0.1) diff --git a/examples/fluid/falling_water_spheres_morris_2d.jl b/examples/fluid/falling_water_spheres_morris_2d.jl new file mode 100644 index 000000000..92df08eb2 --- /dev/null +++ b/examples/fluid/falling_water_spheres_morris_2d.jl @@ -0,0 +1,110 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.05 + +sphere1_center = (0.5, 0.2) +sphere2_center = (1.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -1.0)) +# sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, +# fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -1.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.005 +alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +# sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, +# fluid_smoothing_length, +# sound_speed, viscosity=viscosity, +# density_calculator=ContinuityDensity(), +# acceleration=(0.0, -gravity), +# reference_particle_spacing=fluid_particle_spacing, +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728)) + + +sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, SummationDensity(), + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=ArtificialViscosityMonaghan(alpha=alpha, + beta=0.0), + surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, + fluid_smoothing_length), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient= 0.0728), + correction=GradientCorrection(), reference_particle_spacing=fluid_particle_spacing, + ) + +# sphere = EntropicallyDampedSPHSystem(sphere2, fluid_smoothing_kernel, +# fluid_smoothing_length, +# sound_speed, viscosity=viscosity, +# density_calculator=ContinuityDensity(), +# acceleration=(0.0, -gravity), +# reference_particle_spacing=fluid_particle_spacing, +# surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, +# fluid_smoothing_length)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(sphere_surface_tension, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.001, output_directory="out", + prefix="with_correction_no_heur", write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-7, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/periodic_field_2d.jl b/examples/fluid/periodic_field_2d.jl new file mode 100644 index 000000000..692696f5a --- /dev/null +++ b/examples/fluid/periodic_field_2d.jl @@ -0,0 +1,213 @@ +using TrixiParticles +using OrdinaryDiffEq + +# Domain size +domain_width = 0.5 +domain_height = 1.0 +no_particles = 50 + +# Particle spacing +particle_spacing = domain_width / no_particles # Adjust for resolution + +# Numerical settings +smoothing_length = 1.2 * particle_spacing +sound_speed = 20.0 + +# Fluid properties +fluid_density = 1.0 +# No gravity +gravity = (0.0, 0.0) + +# Time span +tspan = (0.0, 1.0) + +rect_size = (domain_width, domain_width) + +color0 = RectangularShape(particle_spacing, + round.(Int, rect_size ./ particle_spacing), + zeros(length(rect_size)), density=fluid_density) + +color1 = RectangularShape(particle_spacing, + round.(Int, rect_size ./ particle_spacing), + (0.0, domain_width), density=fluid_density) + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, + clip_negative_pressure=true) + +wcsph_color0 = WeaklyCompressibleSPHSystem(color0, SummationDensity(), + state_equation, SchoenbergCubicSplineKernel{2}(), + smoothing_length, + reference_particle_spacing=particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=1.0), + color_value=0) + +wcsph_color1 = WeaklyCompressibleSPHSystem(color1, SummationDensity(), + state_equation, SchoenbergCubicSplineKernel{2}(), + smoothing_length, + reference_particle_spacing=particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=1.0), + color_value=1) + +periodic_box = PeriodicBox(min_corner=[0.0, 0.0], max_corner=[domain_width, domain_height]) +semi = Semidiscretization(wcsph_color0, wcsph_color1, + neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing, + periodic_box=periodic_box)) +ode = semidiscretize(semi, tspan, data_type=nothing) + +info_callback = InfoCallback(interval=100) + +solution_prefix = "" +saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix) + +# This can be overwritten with `trixi_include` +extra_callback = nothing + +use_reinit = false +stepsize_callback = StepsizeCallback(cfl=0.9) + +callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); + +# # Create particle coordinates +# nx = Int(domain_width / particle_spacing) + 1 +# ny = Int(domain_height / particle_spacing) + 1 +# x_coords = range(0.0, stop=domain_width, length=nx) +# y_coords = range(0.0, stop=domain_height, length=ny) + +# coordinates = [ [x, y] for y in y_coords, x in x_coords ] +# coordinates = reduce(vcat, coordinates)' + +# # Assign colors based on y-coordinate +# colors = [ coord[2] > 0.5 ? 1 : 0 for coord in eachcol(coordinates) ] + +# # Create particle properties +# n_particles = size(coordinates, 2) +# particle_mass = fluid_density * particle_spacing^2 + +# particles = Particles( +# coordinates = coordinates, +# velocities = zeros(2, n_particles), +# masses = fill(particle_mass, n_particles), +# densities = fill(fluid_density, n_particles), +# colors = colors +# ) + +# # Initialize random velocities +# # Internal energy per particle +# internal_energy = 0.5 * sound_speed^2 + +# # Desired kinetic energy per particle +# desired_ke_per_particle = internal_energy * 1e-6 + +# # Generate random velocities +# using Random +# Random.seed!(1234) # For reproducibility + +# velocities = zeros(2, n_particles) +# for i in 1:n_particles +# # Random velocity direction +# theta = 2 * π * rand() +# # Random velocity magnitude +# v_mag = sqrt(2 * desired_ke_per_particle) +# velocities[:, i] = v_mag * [cos(theta), sin(theta)] +# end + +# # Assign velocities to particles +# particles.velocities = velocities + +# # Use appropriate density calculator +# fluid_density_calculator = ContinuityDensity() + +# # Exclude viscosity +# viscosity = nothing + +# # Create the fluid system +# fluid_system = WeaklyCompressibleSPHSystem( +# particles, +# fluid_density_calculator, +# StateEquationCole( +# sound_speed = sound_speed, +# reference_density = fluid_density, +# exponent = 7.0, +# clip_negative_pressure = false +# ), +# SmoothingKernelCubicSpline(), +# smoothing_length, +# viscosity = viscosity, +# acceleration = gravity, +# surface_tension = SurfaceTensionMomentumMorris() +# ) + +# # Define the periodic box matching the simulation domain +# periodic_box = PeriodicBox( +# min_corner = [0.0, 0.0], +# max_corner = [domain_width, domain_height] +# ) + +# # Configure the neighborhood search with the periodic box +# neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box) + +# # Set up the semidiscretization and ODE problem +# semi = Semidiscretization( +# fluid_system, +# neighborhood_search = neighborhood_search +# ) + +# ode = semidiscretize(semi, tspan) + +# # Define callbacks to record kinetic and internal energy +# kinetic_energies = Float64[] +# internal_energies = Float64[] +# times = Float64[] + +# function record_energies!(integrator) +# u = integrator.u +# v = integrator.cache + +# # Compute kinetic energy +# velocities = v.velocities +# masses = particles.masses +# ke = 0.5 * sum(masses .* sum(velocities .^ 2, dims=1)) +# push!(kinetic_energies, ke) + +# # Compute internal energy +# densities = v.densities +# pressures = v.pressures +# ie = sum(pressures .* masses / (fluid_density * sound_speed^2)) +# push!(internal_energies, ie) + +# push!(times, integrator.t) +# return +# end + +# cb = CallbackSet(SavingCallback(record_energies!)) + +# # Time integration settings +# dtmax = 0.25 * particle_spacing / (sound_speed + maximum(abs.(velocities))) + +# # Solve the ODE problem +# sol = solve( +# ode, +# RDPK3SpFSAL35(), +# callback = cb, +# save_everystep = false, +# abstol = 1e-6, +# reltol = 1e-6, +# dtmax = dtmax +# ) + +# # Compute the ratio of kinetic energy to internal energy +# energy_ratios = kinetic_energies ./ internal_energies + +# # Plot the ratio over time +# plot( +# times, energy_ratios, +# xlabel = "Time", +# ylabel = "KE / Internal Energy", +# title = "Parasitic Currents Test", +# legend = false +# ) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index af558133f..5254010eb 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -5,16 +5,16 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -particle_spacing = 0.05 +particle_spacing = 0.02 # Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 3 +boundary_layers = 5 # Make sure that the kernel support of fluid particles at an open boundary is always # fully sampled. # Note: Due to the dynamics at the inlets and outlets of open boundaries, # it is recommended to use `open_boundary_layers > boundary_layers` -open_boundary_layers = 6 +open_boundary_layers = 8 # ========================================================================================== # ==== Experiment Setup @@ -34,12 +34,11 @@ fluid_density = 1000.0 # For this particular example, it is necessary to have a background pressure. # Otherwise the suction at the outflow is to big and the simulation becomes unstable. -pressure = 1000.0 +pressure = 0.0 -sound_speed = 10 * prescribed_velocity +sound_speed = 20 * prescribed_velocity -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7, background_pressure=pressure) +state_equation = nothing pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, pressure=pressure, n_layers=boundary_layers, @@ -48,7 +47,7 @@ pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_densi # Shift pipe walls in negative x-direction for the inflow pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers -n_buffer_particles = 4 * pipe.n_particles_per_dimension[2] +n_buffer_particles = 40 * pipe.n_particles_per_dimension[2] # ========================================================================================== # ==== Fluid @@ -67,6 +66,8 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi buffer_size=n_buffer_particles) # Alternatively the WCSPH scheme can be used +# state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, +# exponent=1, background_pressure=pressure) # alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) # viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) @@ -77,33 +78,34 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi # ========================================================================================== # ==== Open Boundary -function velocity_function(pos, t) +function velocity_function2d(pos, t) # Use this for a time-dependent inflow velocity # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) return SVector(prescribed_velocity, 0.0) end -inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, - open_boundary_layers, density=fluid_density, particle_spacing) +inflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, domain_size[2]]), + plane_normal=flow_direction, open_boundary_layers, + density=fluid_density, particle_spacing, boundary_type=:inflow) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelLastiwka(), + boundary_model=BoundaryModelTafuni(), buffer_size=n_buffer_particles, reference_density=fluid_density, reference_pressure=pressure, - reference_velocity=velocity_function) + reference_velocity=velocity_function2d) -outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), - flow_direction, open_boundary_layers, density=fluid_density, - particle_spacing) +outflow = BoundaryZone(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), + plane_normal=-flow_direction, open_boundary_layers, + density=fluid_density, particle_spacing, boundary_type=:outflow) open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=BoundaryModelLastiwka(), + boundary_model=BoundaryModelTafuni(), buffer_size=n_buffer_particles, reference_density=fluid_density, reference_pressure=pressure, - reference_velocity=velocity_function) + reference_velocity=velocity_function2d) # ========================================================================================== # ==== Boundary @@ -111,7 +113,7 @@ open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, AdamiPressureExtrapolation(), state_equation=state_equation, - #viscosity=ViscosityAdami(nu=1e-4), + viscosity=ViscosityAdami(nu=1e-4), smoothing_kernel, smoothing_length) boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl new file mode 100644 index 000000000..71506c37b --- /dev/null +++ b/examples/fluid/pipe_flow_3d.jl @@ -0,0 +1,46 @@ +# 3D channel flow simulation with open boundaries. +using TrixiParticles +using OrdinaryDiffEq + +tspan = (0.0, 2.0) + +# load variables +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + sol=nothing) + +function velocity_function3d(pos, t) + # Use this for a time-dependent inflow velocity + # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) + + return SVector(prescribed_velocity, 0.0, 0.0) +end + +domain_size = (1.0, 0.4, 0.4) + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) + +pipe3d = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true, true, true)) + +flow_direction = [1.0, 0.0, 0.0] + +# setup simulation +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + domain_size=domain_size, flow_direction=flow_direction, + pipe=pipe3d, + n_buffer_particles=4 * pipe3d.n_particles_per_dimension[2]^2, + smoothing_kernel=WendlandC2Kernel{3}(), + reference_velocity=velocity_function3d, + inflow=InFlow(; + plane=([0.0, 0.0, 0.0], + [0.0, domain_size[2], 0.0], + [0.0, 0.0, domain_size[3]]), flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing), + outflow=OutFlow(; + plane=([domain_size[1], 0.0, 0.0], + [domain_size[1], domain_size[2], 0.0], + [domain_size[1], 0.0, domain_size[3]]), + flow_direction, open_boundary_layers, density=fluid_density, + particle_spacing)) diff --git a/examples/fluid/rain_flat.jl b/examples/fluid/rain_flat.jl new file mode 100644 index 000000000..b0d4296be --- /dev/null +++ b/examples/fluid/rain_flat.jl @@ -0,0 +1,128 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.00004 # x75 visc +# fluid_particle_spacing = 0.0001 # no improvement +boundary_particle_spacing = 0.00005 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 0.02) +fluid_density = 1000.0 +sound_speed = 100 + +# Boundary geometry and initial fluid particle positions +plate_length = 0.01 +plate_width = 0.01 + +sphere_radius = 0.0005 +#terminal_velocity = 75 * 5.8 # match Re +terminal_velocity = 5.8 # real +#terminal_velocity = 1.0 + +boundary_thickness = 4 * boundary_particle_spacing + +# ground floor +plate_size = (plate_length, plate_width, boundary_thickness) + +plate = RectangularShape(boundary_particle_spacing, + round.(Int, plate_size ./ boundary_particle_spacing), + (0.0, 0.0, 0.0), density=fluid_density) + +sphere1_center = (0.005, 0.005, sphere_radius + boundary_thickness) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -terminal_velocity)) + +# box = RectangularTank(fluid_particle_spacing, (0.3, 0.125, 0.0), (0.35, 0.075, 0.2), fluid_density, +# n_layers=8, spacing_ratio=spacing_ratio, acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, +# faces=(false, false, true, false, true, true), +# min_coordinates=(tank_size[1], 0.0, 0.0)) + +# # move to the end +# for i in axes(end_tank.boundary.coordinates, 2) +# end_tank.boundary.coordinates[:, i] .+= [tank_size[1], 0.0] +# end + +#tank = union(tank, tank_end, tank_side_mz, tank_side_pz, tank_side_mz_end, tank_side_pz_end) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() +# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +# kinematic viscosity of water at 20C +nu_water = 8.9E-7 + +# Morris has a higher velocity with same viscosity. +# viscosity = ViscosityMorris(nu=75*nu_water) +viscosity = ViscosityAdami(nu=20*nu_water) + +density_diffusion = DensityDiffusionAntuono(sphere1, delta=0.1) + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +fluid_system = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, 0.0, -gravity), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728, + free_surface_threshold=0.6), # 0.55 too many # 0.8 even more + reference_particle_spacing=fluid_particle_spacing, + surface_normal_method=ColorfieldSurfaceNormal(smoothing_kernel, + smoothing_length, boundary_contact_threshold=1.0)) + + +# fluid_system = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, 0.0, -gravity)) + +# fluid_system = EntropicallyDampedSPHSystem(sphere1, smoothing_kernel, +# smoothing_length, +# sound_speed, +# viscosity=viscosity, +# density_calculator=fluid_density_calculator, +# reference_particle_spacing=fluid_particle_spacing, +# acceleration=(0.0, 0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(plate.density, plate.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system = BoundarySPHSystem(plate, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.00005, prefix="rainfall_morris_alpha001_h40micro_nu20_lowerdt") +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-4, # Limit stepsize to prevent crashing + dt=1e-8, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl index 9006ebf0e..2e8b6efe4 100644 --- a/examples/fluid/sphere_surface_tension_2d.jl +++ b/examples/fluid/sphere_surface_tension_2d.jl @@ -5,7 +5,9 @@ using OrdinaryDiffEq fluid_density = 1000.0 -particle_spacing = 0.1 +particle_spacing = 0.05 +# Use a higher resolution for a better result +# particle_spacing = 0.025 # Note: Only square shapes will result in a sphere. # Furthermore, changes of the coefficients might be necessary for higher resolutions or larger squares. @@ -20,9 +22,10 @@ state_equation = StateEquationCole(; sound_speed, reference_density=fluid_densit # smoothing_kernel = WendlandC2Kernel{2}() # nu = 0.01 -smoothing_length = 1.0 * particle_spacing -smoothing_kernel = SchoenbergCubicSplineKernel{2}() -nu = 0.025 +smoothing_length = 3.5 * particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() +# nu = 0.001 # SurfaceTensionMomentumMorris +nu = 0.05 # SurfaceTensionMorris fluid = RectangularShape(particle_spacing, round.(Int, fluid_size ./ particle_spacing), zeros(length(fluid_size)), density=fluid_density) @@ -30,25 +33,51 @@ fluid = RectangularShape(particle_spacing, round.(Int, fluid_size ./ particle_sp alpha = 8 * nu / (smoothing_length * sound_speed) source_terms = SourceTermDamping(; damping_coefficient=0.5) fluid_system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), - state_equation, smoothing_kernel, + state_equation, fluid_smoothing_kernel, smoothing_length, viscosity=ArtificialViscosityMonaghan(alpha=alpha, beta=0.0), - surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), - correction=AkinciFreeSurfaceCorrection(fluid_density), - source_terms=source_terms) + surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, + smoothing_length), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50 * + 0.0728), + correction=GradientCorrection(), reference_particle_spacing=particle_spacing, + ) + +# fluid_system = EntropicallyDampedSPHSystem(fluid, fluid_smoothing_kernel, +# smoothing_length, +# sound_speed, +# viscosity=ViscosityMorris(nu=nu), +# density_calculator=ContinuityDensity(), +# reference_particle_spacing=particle_spacing, +# acceleration=zeros(length(fluid_size)), +# surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, +# smoothing_length), +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50 * +# 0.0728)) + +# fluid_system = EntropicallyDampedSPHSystem(fluid, fluid_smoothing_kernel, +# smoothing_length, +# sound_speed, +# viscosity=ViscosityMorris(nu=nu), +# density_calculator=ContinuityDensity(), +# reference_particle_spacing=particle_spacing, +# acceleration=zeros(length(fluid_size)), +# surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, +# smoothing_length), +# surface_tension=SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0)) # ========================================================================================== # ==== Simulation semi = Semidiscretization(fluid_system) -tspan = (0.0, 3.0) +tspan = (0.0, 50.0) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) # For overwriting via `trixi_include` -saving_callback = SolutionSavingCallback(dt=0.02) +saving_callback = SolutionSavingCallback(dt=1.0) stepsize_callback = StepsizeCallback(cfl=1.0) diff --git a/examples/fluid/sphere_surface_tension_3d.jl b/examples/fluid/sphere_surface_tension_3d.jl index 58a506294..6c94e2904 100644 --- a/examples/fluid/sphere_surface_tension_3d.jl +++ b/examples/fluid/sphere_surface_tension_3d.jl @@ -5,7 +5,7 @@ using OrdinaryDiffEq fluid_density = 1000.0 -particle_spacing = 0.1 +particle_spacing = 0.05 fluid_size = (0.9, 0.9, 0.9) sound_speed = 20.0 @@ -13,13 +13,12 @@ sound_speed = 20.0 # For all surface tension simulations, we need a compact support of `2 * particle_spacing` smoothing_length = 1.0 * particle_spacing -nu = 0.01 +nu = 0.04 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "sphere_surface_tension_2d.jl"), - dt=0.1, cfl=1.2, surface_tension_coefficient=0.1, - tspan=(0.0, 10.0), nu=nu, - alpha=10 * nu / (smoothing_length * sound_speed), - smoothing_kernel=SchoenbergCubicSplineKernel{3}(), + surface_tension_coefficient=0.5, dt=0.25, + tspan=(0.0, 100.0), nu=nu, smoothing_length=3.0 * particle_spacing, + fluid_smoothing_kernel=WendlandC2Kernel{3}(), particle_spacing=particle_spacing, sound_speed=sound_speed, fluid_density=fluid_density, fluid_size=fluid_size) diff --git a/examples/fluid/sphere_surface_tension_wall_2d.jl b/examples/fluid/sphere_surface_tension_wall_2d.jl index 3e1cbdb80..7158a5d73 100644 --- a/examples/fluid/sphere_surface_tension_wall_2d.jl +++ b/examples/fluid/sphere_surface_tension_wall_2d.jl @@ -9,6 +9,7 @@ fluid_particle_spacing = 0.0025 # ========================================================================================== # ==== Experiment Setup +gravity = 9.81 tspan = (0.0, 0.5) # Boundary geometry and initial fluid particle positions @@ -17,7 +18,10 @@ tank_size = (0.5, 0.1) fluid_density = 1000.0 sound_speed = 120.0 -sphere_radius = 0.05 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +sphere_radius = 0.04 sphere1_center = (0.25, sphere_radius) sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, @@ -25,7 +29,11 @@ sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, # ========================================================================================== # ==== Fluid -fluid_smoothing_length = 1.0 * fluid_particle_spacing +# fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps() +# fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_smoothing_length = 3.25 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() # For perfect wetting # nu = 0.0005 @@ -35,10 +43,46 @@ nu = 0.001 alpha = 8 * nu / (fluid_smoothing_length * sound_speed) # `adhesion_coefficient = 1.0` and `surface_tension_coefficient = 0.01` for perfect wetting # `adhesion_coefficient = 0.001` and `surface_tension_coefficient = 2.0` for no wetting + +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, ContinuityDensity(), +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=10.0, +# free_surface_threshold=0.75), +# reference_particle_spacing=fluid_particle_spacing, +# surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, +# fluid_smoothing_length, +# boundary_contact_threshold=0.05)) + +sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, ContinuityDensity(), + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0, + free_surface_threshold=0.75), + reference_particle_spacing=fluid_particle_spacing, + surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, + fluid_smoothing_length, + boundary_contact_threshold=0.05)) + +# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, ContinuityDensity(), +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=2.0), +# correction=AkinciFreeSurfaceCorrection(fluid_density), +# reference_particle_spacing=fluid_particle_spacing) + trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"), sphere=nothing, sphere1=sphere1, adhesion_coefficient=0.001, - wall_viscosity=4.0 * nu, surface_tension_coefficient=2.0, alpha=alpha, + wall_viscosity=4.0 * nu, surface_tension_coefficient=0.9, alpha=alpha, sound_speed=sound_speed, fluid_density=fluid_density, nu=nu, fluid_particle_spacing=fluid_particle_spacing, tspan=tspan, - tank_size=tank_size) + tank_size=tank_size, fluid_smoothing_length=fluid_smoothing_length, + sphere_surface_tension=sphere_surface_tension) diff --git a/examples/fluid/static_spheres.jl b/examples/fluid/static_spheres.jl new file mode 100644 index 000000000..e5fb8683b --- /dev/null +++ b/examples/fluid/static_spheres.jl @@ -0,0 +1,98 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 0.0 +tspan = (0.0, 0.3) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.5, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.05 + +sphere1_center = (0.5, 0.5) +sphere2_center = (1.5, 0.5) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=RoundSphere(), velocity=(0.0, 0.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=RoundSphere(), velocity=(0.0, 0.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.005 +alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + reference_particle_spacing=fluid_particle_spacing, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728)) + +sphere = EntropicallyDampedSPHSystem(sphere2, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + reference_particle_spacing=fluid_particle_spacing, + acceleration=(0.0, -gravity), + surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, + fluid_smoothing_length)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="static", write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-7, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/waterfall.jl b/examples/fluid/waterfall.jl new file mode 100644 index 000000000..516ffc6c7 --- /dev/null +++ b/examples/fluid/waterfall.jl @@ -0,0 +1,127 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.008 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 7 +spacing_ratio = 1 + +boundary_particle_spacing = fluid_particle_spacing / spacing_ratio + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 5.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.2, 0.1, 0.2) +tank_size = (0.3, 0.15, 0.20) +end_size = (0.3, 0.05, 0.20) + +fluid_density = 1000.0 +sound_speed = 40 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=1000) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=4, spacing_ratio=spacing_ratio, + acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, + faces=(true, false, true, false, true, true)) + +end_tank = RectangularTank(fluid_particle_spacing, (0.0, 0.0, 0.0), end_size, fluid_density, + n_layers=8, spacing_ratio=spacing_ratio, + acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, + faces=(false, false, true, false, true, true), + min_coordinates=(tank_size[1], 0.0, 0.0)) + +# # move to the end +# for i in axes(end_tank.boundary.coordinates, 2) +# end_tank.boundary.coordinates[:, i] .+= [tank_size[1], 0.0] +# end + +# tank = union(tank, end_tank) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() +# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +nu_water = 8.9E-7 + +# Morris has a higher velocity with same viscosity. +# viscosity = ViscosityMorris(nu=100*nu_water) +viscosity = ViscosityAdami(nu=100 * nu_water) + +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +# fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity, 0.0)) + +fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, + smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity, 0.0), + # surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), + reference_particle_spacing=fluid_particle_spacing, + buffer_size=1) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model2 = BoundaryModelDummyParticles(end_tank.boundary.density, + end_tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system2 = BoundarySPHSystem(end_tank.boundary, boundary_model2) + +outflow = BoundaryZone(; plane=([0.4, -0.2, -0.1], [0.4, -0.2, 0.3], [0.8, -0.2, 0.3]), + plane_normal=[0.0, -1.0, 0.0], open_boundary_layers=1, + density=2 * eps(), particle_spacing=fluid_particle_spacing, + boundary_type=:outflow) + +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=BasicOutlet()) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system, boundary_system2, + open_boundary_out) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="lower_visc") +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +# 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-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/waterfallv2.jl b/examples/fluid/waterfallv2.jl new file mode 100644 index 000000000..4adf63fc5 --- /dev/null +++ b/examples/fluid/waterfallv2.jl @@ -0,0 +1,170 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0025 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 5.0) +fluid_density = 1000.0 +sound_speed = 80 + +# Boundary geometry and initial fluid particle positions +canal_width = 0.1 +canal_height = 0.1 +canal_length = 0.1 +outflow_length = 0.3 + +initial_fluid_size = (canal_length, canal_height, canal_width) +boundary_thickness = 4 * fluid_particle_spacing + +# ground floor +plate = (canal_length + outflow_length, boundary_thickness, canal_width) +# end to close the canal +plate_end = (boundary_thickness, canal_height + boundary_thickness, canal_width) + +plate_side = (1.5 * canal_length, canal_height + boundary_thickness, boundary_thickness) +# goes from the end of plate_side to the end of the plate +plate_side_end = (plate[1] - 1.5 * canal_length + boundary_thickness, + 0.5 * canal_height + boundary_thickness, boundary_thickness) + +fluid = RectangularShape(fluid_particle_spacing, + round.(Int, initial_fluid_size ./ fluid_particle_spacing), + (0.0, plate[2], 0.0), density=fluid_density) + +tank = RectangularShape(fluid_particle_spacing, + round.(Int, plate ./ fluid_particle_spacing), + (0.0, 0.0, 0.0), density=fluid_density) + +tank_end = RectangularShape(fluid_particle_spacing, + round.(Int, plate_end ./ fluid_particle_spacing), + (-plate_end[1], 0.0, 0.0), density=fluid_density) + +tank_side_pz = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side ./ fluid_particle_spacing), + (-plate_end[1], 0.0, initial_fluid_size[1]), + density=fluid_density) + +tank_side_mz = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side ./ fluid_particle_spacing), + (-plate_end[1], 0.0, -plate_end[1]), density=fluid_density) + +tank_side_pz_end = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side_end ./ fluid_particle_spacing), + (plate_side[1] - boundary_thickness, 0.0, + initial_fluid_size[1]), density=fluid_density) + +tank_side_mz_end = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side_end ./ fluid_particle_spacing), + (plate_side[1] - boundary_thickness, 0.0, + -plate_end[1]), density=fluid_density) + +# box = RectangularTank(fluid_particle_spacing, (0.3, 0.125, 0.0), (0.35, 0.075, 0.2), fluid_density, +# n_layers=8, spacing_ratio=spacing_ratio, acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, +# faces=(false, false, true, false, true, true), +# min_coordinates=(tank_size[1], 0.0, 0.0)) + +# # move to the end +# for i in axes(end_tank.boundary.coordinates, 2) +# end_tank.boundary.coordinates[:, i] .+= [tank_size[1], 0.0] +# end + +tank = union(tank, tank_end, tank_side_mz, tank_side_pz, tank_side_mz_end, tank_side_pz_end) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.25 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() +# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +nu_water = 8.9E-7 + +# Morris has a higher velocity with same viscosity. +# viscosity = ViscosityMorris(nu=100*nu_water) +viscosity = ViscosityAdami(nu=50 * nu_water) + +# density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +# fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity, 0.0)) + +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity, 0.0), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50*0.0728, + free_surface_threshold=0.7), + reference_particle_spacing=fluid_particle_spacing, + surface_normal_method=ColorfieldSurfaceNormal(smoothing_kernel, + smoothing_length, + boundary_contact_threshold=0.1)) + +# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, +# smoothing_length, +# sound_speed, viscosity=viscosity, +# density_calculator=ContinuityDensity(), +# acceleration=(0.0, -gravity, 0.0), +# # surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), +# reference_particle_spacing=fluid_particle_spacing) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.density, tank.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system = BoundarySPHSystem(tank, boundary_model) + +# boundary_model2 = BoundaryModelDummyParticles(end_tank.boundary.density, +# end_tank.boundary.mass, +# state_equation=state_equation, +# boundary_density_calculator, +# smoothing_kernel, smoothing_length, +# viscosity=viscosity) + +# boundary_system2 = BoundarySPHSystem(end_tank.boundary, boundary_model2) + +# outflow = BoundaryZone(; plane=([0.4, -0.2, -0.1], [0.4, -0.2, 0.3], [0.8, -0.2, 0.3]), +# plane_normal=[0.0, -1.0, 0.0], open_boundary_layers=1, +# density=2 * eps(), particle_spacing=fluid_particle_spacing, +# boundary_type=:outflow) + +# open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, +# boundary_model=BasicOutlet()) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="waterfall_v2_less_height") +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +# 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-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d6e0c2996..face09156 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -54,8 +54,7 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow, - OutFlow + BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, BoundaryZone export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback export ContinuityDensity, SummationDensity @@ -63,12 +62,12 @@ export PenaltyForceGanzenmueller export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel -export StateEquationCole +export StateEquationCole, StateEquationIdealGas export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, BoundaryModelLastiwka + PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk @@ -83,6 +82,8 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk -export SurfaceTensionAkinci, CohesionForceAkinci +export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris, + SurfaceTensionMomentumMorris +export ColorfieldSurfaceNormal end # module diff --git a/src/general/buffer.jl b/src/general/buffer.jl index aae5ccb5c..161d774d6 100644 --- a/src/general/buffer.jl +++ b/src/general/buffer.jl @@ -20,13 +20,13 @@ function allocate_buffer(initial_condition, buffer::SystemBuffer) # Initialize particles far away from simulation domain coordinates = fill(1e16, ndims(initial_condition), buffer_size) - if all(rho -> isapprox(rho, first(initial_condition.density), atol=eps(), rtol=eps()), - initial_condition.density) - density = first(initial_condition.density) - else - throw(ArgumentError("`initial_condition.density` needs to be constant when using `SystemBuffer`")) + if !all(rho -> isapprox(rho, first(initial_condition.density), atol=eps(), rtol=eps()), + initial_condition.density) + @warn "allocated buffer is using constant values" end + density = first(initial_condition.density) + particle_spacing = initial_condition.particle_spacing buffer_ic = InitialCondition(; coordinates, density, particle_spacing) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..1057e307e 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -407,7 +407,7 @@ end function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) (; systems) = semi - return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) + return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system, semi), systems) end function drift!(du_ode, v_ode, u_ode, semi, t) @@ -805,6 +805,22 @@ function update_nhs!(neighborhood_search, return neighborhood_search end +function update_nhs!(neighborhood_search, + system::TotalLagrangianSPHSystem, + neighbor::OpenBoundarySPHSystem, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search +end + +function update_nhs!(neighborhood_search, + system::OpenBoundarySPHSystem, + neighbor::TotalLagrangianSPHSystem, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search +end + # Forward to PointNeighbors.jl function update!(neighborhood_search, system, x, y; points_moving=(true, false)) PointNeighbors.update!(neighborhood_search, x, y; points_moving) @@ -818,13 +834,49 @@ function update!(neighborhood_search, system::GPUSystem, x, y; points_moving=(tr end function check_configuration(systems) + min_color = 0 + max_color = 0 + no_fluid_and_bnd_systems = 0 + uses_surface_tension_model = false + foreach_system(systems) do system check_configuration(system, systems) + + if system isa FluidSystem && !isnothing(system.surface_tension) + uses_surface_tension_model = true + end + + if system isa FluidSystem || system isa BoundarySPHSystem + no_fluid_and_bnd_systems += 1 + if system.color < min_color + min_color = system.color + end + + if system.color > max_color + max_color = system.color + end + end + end + + if max_color == 0 && min_color == 0 && no_fluid_and_bnd_systems > 1 && + uses_surface_tension_model + throw(ArgumentError("If a surface tension model is used the values of at least one system needs to have a color different than 0.")) end end check_configuration(system, systems) = nothing +function check_configuration(fluid_system::FluidSystem, systems) + if !isnothing(fluid_system.surface_tension) + foreach_system(systems) do neighbor + if neighbor isa FluidSystem && isnothing(fluid_system.surface_tension) && + isnothing(fluid_system.surface_normal_method) + throw(ArgumentError("All `FluidSystem` need to use a surface tension model or a surface normal method.")) + end + end + end +end + function check_configuration(boundary_system::BoundarySPHSystem, systems) (; boundary_model) = boundary_system @@ -854,3 +906,13 @@ function check_configuration(system::TotalLagrangianSPHSystem, systems) "`ContinuityDensity` is not yet supported for a `TotalLagrangianSPHSystem`")) end end + +function check_configuration(system::OpenBoundarySPHSystem, systems) + (; boundary_model, boundary_zone) = system + + if boundary_model isa BoundaryModelLastiwka && + first(typeof(boundary_zone).parameters) === BidirectionalFlow + throw(ArgumentError("`BoundaryModelLastiwka` needs a specific flow direction. " * + "Please specify inflow and outflow.")) + end +end diff --git a/src/general/system.jl b/src/general/system.jl index 674555aea..347aebefb 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -147,3 +147,17 @@ end # Only for systems requiring a mandatory callback reset_callback_flag!(system) = system + +# Assuming a constant particle spacing one can calculate the number of neighbors within the +# compact support for an undisturbed particle distribution. +function number_density(::Val{D}, particle_spacing, compact_support) where {D} + throw(ArgumentError("Unsupported dimension: $D")) +end + +@inline function number_density(::Val{2}, particle_spacing, compact_support) + return floor(Int64, pi * compact_support^2 / particle_spacing^2) +end + +@inline @fastpow function number_density(::Val{3}, particle_spacing, compact_support) + return floor(Int64, 4.0 / 3.0 * pi * compact_support^3 / particle_spacing^3) +end diff --git a/src/schemes/boundary/boundary.jl b/src/schemes/boundary/boundary.jl index 6d50eec9b..a09545a9b 100644 --- a/src/schemes/boundary/boundary.jl +++ b/src/schemes/boundary/boundary.jl @@ -1,7 +1,10 @@ include("dummy_particles/dummy_particles.jl") include("system.jl") include("open_boundary/boundary_zones.jl") +include("open_boundary/mirroring.jl") include("open_boundary/method_of_characteristics.jl") include("open_boundary/system.jl") # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem` # and the `TotalLagrangianSPHSystem` and are therefore included later. + +@inline Base.ndims(boundary_model::BoundaryModelDummyParticles) = ndims(boundary_model.smoothing_kernel) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 1cd8e0dfe..af29289bb 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -16,12 +16,12 @@ Boundary model for `BoundarySPHSystem`. - `smoothing_length`: Smoothing length should be the same as for the adjacent fluid system. # Keywords -- `state_equation`: This should be the same as for the adjacent fluid system - (see e.g. [`StateEquationCole`](@ref)). -- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). -- `viscosity`: Slip (default) or no-slip condition. See description below for further - information. - +- `state_equation`: This should be the same as for the adjacent fluid system + (see e.g. [`StateEquationCole`](@ref)). +- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). +- `viscosity`: Slip (default) or no-slip condition. See description below for further + information. +- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary. # Examples ```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1) # Free-slip condition @@ -44,6 +44,7 @@ struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C} density_calculator :: DC smoothing_kernel :: K smoothing_length :: ELTYPE + number_density :: Int64 viscosity :: V correction :: COR cache :: C @@ -54,19 +55,33 @@ end function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, - state_equation=nothing, correction=nothing) + state_equation=nothing, correction=nothing, + reference_particle_spacing=0.0) pressure = initial_boundary_pressure(initial_density, density_calculator, state_equation) NDIMS = ndims(smoothing_kernel) + ELTYPE = eltype(smoothing_length) n_particles = length(initial_density) cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., create_cache_model(initial_density, density_calculator)..., - create_cache_model(correction, initial_density, NDIMS, n_particles)...) + create_cache_model(correction, initial_density, NDIMS, n_particles)..., + (; colorfield_bnd=zeros(ELTYPE, n_particles), + colorfield=zeros(ELTYPE, n_particles), + neighbor_count=zeros(ELTYPE, n_particles))...) + + number_density_ = 0 + if reference_particle_spacing > 0.0 + number_density_ = number_density(Val(ndims(boundary_model)), + reference_particle_spacing, + compact_support(smoothing_kernel, + smoothing_length)) + end return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation, density_calculator, smoothing_kernel, - smoothing_length, viscosity, correction, cache) + smoothing_length, number_density_, viscosity, + correction, cache) end @doc raw""" diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index cc9a77157..d2d40cdc6 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -1,17 +1,23 @@ +struct BidirectionalFlow end + +struct InFlow end + +struct OutFlow end + @doc raw""" - InFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) + BoundaryZone(; plane, plane_normal, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer, boundary_type=:bidirectional_flow) -Inflow boundary zone for [`OpenBoundarySPHSystem`](@ref). +Boundary zone for [`OpenBoundarySPHSystem`](@ref). -The specified plane (line in 2D or rectangle in 3D) will be extruded in upstream -direction (the direction opposite to `flow_direction`) to create a box for the boundary zone. -There are three ways to specify the actual shape of the inflow: +The specified plane (line in 2D or rectangle in 3D) will be extruded in the direction +opposite to `plane_normal` to create a box for the boundary zone. +There are three ways to specify the actual shape of the boundary zone: 1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then - be filled with inflow particles (default). + be filled with boundary particles (default). 2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D, - which is then extruded in upstream direction to create the inflow particles. + which is then extruded in the direction opposite to `plane_normal` to create the boundary particles. - In 2D, the shape must be either an initial condition with 2D coordinates, which lies on the line specified by `plane`, or an initial condition with 1D coordinates, which lies on the line specified by `plane` when a y-coordinate of `0` is added. @@ -19,7 +25,7 @@ There are three ways to specify the actual shape of the inflow: in the rectangle specified by `plane`, or an initial condition with 2D coordinates, which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added. 3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D, - which will be used for the inflow particles. + which will be used for the boundary particles. !!! note "Note" Particles outside the boundary zone box will be removed. @@ -34,7 +40,11 @@ There are three ways to specify the actual shape of the inflow: In 3D, pass three points ``(A, B, C)``, so that the rectangular inflow surface is spanned by the vectors ``\widehat{AB}`` and ``\widehat{AC}``. These two vectors must be orthogonal. -- `flow_direction`: Vector defining the flow direction. +- `plane_normal`: Vector defining the plane normal. It always points inside the fluid domain. +- `boundary_type=:bidirectional_flow`: Specify the type of the boundary. Available types are + - `:inflow` for an inflow boundary + - `:outflow` for an outflow boundary + - `:bidrectional_flow` (default) for an bidirectional flow boundary - `open_boundary_layers`: Number of particle layers in upstream direction. - `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). - `density`: Particle density (see [`InitialCondition`](@ref)). @@ -49,228 +59,149 @@ There are three ways to specify the actual shape of the inflow: ```julia # 2D plane_points = ([0.0, 0.0], [0.0, 1.0]) -flow_direction=[1.0, 0.0] +plane_normal=[1.0, 0.0] -inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, - open_boundary_layers=4) +inflow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0, + open_boundary_layers=4, boundary_type=:inflow) # 3D plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) -flow_direction=[0.0, 0.0, 1.0] +plane_normal=[0.0, 0.0, 1.0] -inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, - open_boundary_layers=4) +outflow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, density=1.0, + open_boundary_layers=4, boundary_type=:outflow) # 3D particles sampled as cylinder circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) -inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, - extrude_geometry=circle, open_boundary_layers=4) +bidirectional_flow = BoundaryZone(; plane=plane_points, plane_normal, particle_spacing=0.1, + density=1.0, extrude_geometry=circle, open_boundary_layers=4) ``` !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct InFlow{NDIMS, IC, S, ZO, ZW, FD} +struct BoundaryZone{F, NDIMS, IC, S, ZO, ZW, FD, PN} initial_condition :: IC spanning_set :: S zone_origin :: ZO zone_width :: ZW flow_direction :: FD + plane_normal :: PN - function InFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) + function BoundaryZone(; plane, plane_normal, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer, boundary_type=:bidirectional_flow) if open_boundary_layers <= 0 throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) end - # Unit vector pointing in downstream direction - flow_direction_ = normalize(SVector(flow_direction...)) - - # Sample particles in boundary zone - if isnothing(initial_condition) && isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, - density, - direction=-flow_direction_, - n_extrude=open_boundary_layers) - elseif !isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; - particle_spacing, - density, - direction=-flow_direction_, - n_extrude=open_boundary_layers) - end + # `plane_normal` always points in fluid domain + plane_normal_ = normalize(SVector(plane_normal...)) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) + if boundary_type === :bidirectional_flow + boundary_type_ = BidirectionalFlow() + flow_direction = nothing - zone_width = open_boundary_layers * initial_condition.particle_spacing + elseif boundary_type === :inflow + # Unit vector pointing in downstream direction + flow_direction = plane_normal_ - # Vectors spanning the boundary zone/box - spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + boundary_type_ = InFlow() - # First vector of `spanning_vectors` is normal to the inflow plane. - # The normal vector must point in upstream direction for an inflow boundary. - dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + elseif boundary_type === :outflow + # Unit vector pointing in downstream direction + flow_direction = -plane_normal_ - if !isapprox(abs(dot_), 1.0, atol=1e-7) - throw(ArgumentError("`flow_direction` is not normal to inflow plane")) - else - # Flip the normal vector to point in the opposite direction of `flow_direction` - spanning_set[:, 1] .*= -sign(dot_) + boundary_type_ = OutFlow() end - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + ic, spanning_set_, zone_origin, + zone_width = set_up_boundary_zone(plane, plane_normal_, flow_direction, density, + particle_spacing, initial_condition, + extrude_geometry, open_boundary_layers; + boundary_type=boundary_type_) - # Remove particles outside the boundary zone. - # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. - ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) - - return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), - typeof(zone_width), - typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, - flow_direction_) + return new{typeof(boundary_type_), ndims(ic), typeof(ic), typeof(spanning_set_), + typeof(zone_origin), typeof(zone_width), typeof(flow_direction), + typeof(plane_normal_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction, plane_normal_) end end -@doc raw""" - OutFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) - -Outflow boundary zone for [`OpenBoundarySPHSystem`](@ref). - -The specified plane (line in 2D or rectangle in 3D) will be extruded in downstream -direction (the direction in `flow_direction`) to create a box for the boundary zone. -There are three ways to specify the actual shape of the outflow: -1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then - be filled with outflow particles (default). -2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D, - which is then extruded in downstream direction to create the outflow particles. - - In 2D, the shape must be either an initial condition with 2D coordinates, which lies - on the line specified by `plane`, or an initial condition with 1D coordinates, which lies - on the line specified by `plane` when a y-coordinate of `0` is added. - - In 3D, the shape must be either an initial condition with 3D coordinates, which lies - in the rectangle specified by `plane`, or an initial condition with 2D coordinates, - which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added. -3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D, - which will be used for the outflow particles. - -!!! note "Note" - Particles outside the boundary zone box will be removed. - -# Keywords -- `plane`: Tuple of points defining a part of the surface of the domain. - The points must either span a line in 2D or a rectangle in 3D. - This line or rectangle is then extruded in downstream direction to obtain - the boundary zone. - In 2D, pass two points ``(A, B)``, so that the interval ``[A, B]`` is - the outflow surface. - In 3D, pass three points ``(A, B, C)``, so that the rectangular outflow surface - is spanned by the vectors ``\widehat{AB}`` and ``\widehat{AC}``. - These two vectors must be orthogonal. -- `flow_direction`: Vector defining the flow direction. -- `open_boundary_layers`: Number of particle layers in downstream direction. -- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). -- `density`: Particle density (see [`InitialCondition`](@ref)). -- `initial_condition=nothing`: `InitialCondition` for the outflow particles. - Particles outside the boundary zone will be removed. - Do not use together with `extrude_geometry`. -- `extrude_geometry=nothing`: 1D shape in 2D or 2D shape in 3D, which lies on the plane - and is extruded downstream to obtain the outflow particles. - See point 2 above for more details. - -# Examples -```julia -# 2D -plane_points = ([0.0, 0.0], [0.0, 1.0]) -flow_direction = [1.0, 0.0] +@inline Base.ndims(::BoundaryZone{F, NDIMS}) where {F, NDIMS} = NDIMS + +function set_up_boundary_zone(plane, plane_normal, flow_direction, density, + particle_spacing, initial_condition, extrude_geometry, + open_boundary_layers; boundary_type) + if boundary_type isa InFlow + extrude_direction = -flow_direction + elseif boundary_type isa OutFlow + extrude_direction = flow_direction + elseif boundary_type isa BidirectionalFlow + # `plane_normal` is always pointing in the fluid domain + extrude_direction = -plane_normal + end -outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, - open_boundary_layers=4) + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=extrude_direction, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, + density, + direction=extrude_direction, + n_extrude=open_boundary_layers) + end -# 3D -plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) -flow_direction = [0.0, 0.0, 1.0] + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) -outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, - open_boundary_layers=4) + zone_width = open_boundary_layers * initial_condition.particle_spacing -# 3D particles sampled as cylinder -circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) -outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, - extrude_geometry=circle, open_boundary_layers=4) -``` + # First vector of `spanning_vectors` is normal to the boundary plane. + dot_plane_normal = dot(normalize(spanning_set[:, 1]), plane_normal) -!!! warning "Experimental Implementation" - This is an experimental feature and may change in any future releases. -""" -struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} - initial_condition :: IC - spanning_set :: S - zone_origin :: ZO - zone_width :: ZW - flow_direction :: FD - - function OutFlow(; plane, flow_direction, density, particle_spacing, - initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer) - if open_boundary_layers <= 0 - throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) - end - - # Unit vector pointing in downstream direction - flow_direction_ = normalize(SVector(flow_direction...)) - - # Sample particles in boundary zone - if isnothing(initial_condition) && isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, - density, - direction=flow_direction_, - n_extrude=open_boundary_layers) - elseif !isnothing(extrude_geometry) - initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; - particle_spacing, density, - direction=flow_direction_, - n_extrude=open_boundary_layers) - end + if !isapprox(abs(dot_plane_normal), 1.0, atol=1e-7) + throw(ArgumentError("`plane_normal` is not normal to the boundary plane")) + end - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) + if boundary_type isa InFlow + # First vector of `spanning_vectors` is normal to the boundary plane. + dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction) - zone_width = open_boundary_layers * initial_condition.particle_spacing + # The vector must point in upstream direction for an inflow boundary. + # Flip the normal vector to point in the opposite direction of `flow_direction` + spanning_set[:, 1] .*= -sign(dot_flow) - # Vectors spanning the boundary zone/box - spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + elseif boundary_type isa OutFlow + # First vector of `spanning_vectors` is normal to the boundary plane. + dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction) - # First vector of `spanning_vectors` is normal to the outflow plane. - # The normal vector must point in downstream direction for an outflow boundary. - dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + # The vector must point in downstream direction for an outflow boundary. + # Flip the normal vector to point in `flow_direction` + spanning_set[:, 1] .*= sign(dot_flow) - if !isapprox(abs(dot_), 1.0, atol=1e-7) - throw(ArgumentError("`flow_direction` is not normal to outflow plane")) - else - # Flip the normal vector to point in `flow_direction` - spanning_set[:, 1] .*= sign(dot_) - end + elseif boundary_type isa BidirectionalFlow + # Flip the normal vector to point opposite to fluid domain + spanning_set[:, 1] .*= -sign(dot_plane_normal) + end - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) - # Remove particles outside the boundary zone. - # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. - ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + # Remove particles outside the boundary zone. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) - return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), - typeof(zone_width), - typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, - flow_direction_) - end + return ic, spanning_set_, zone_origin, zone_width end -@inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS - function calculate_spanning_vectors(plane, zone_width) return spanning_vectors(Tuple(plane), zone_width), SVector(plane[1]...) end @@ -289,8 +220,9 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) edge1 = plane_points[2] - plane_points[1] edge2 = plane_points[3] - plane_points[1] - if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) - throw(ArgumentError("the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal")) + # Check if the edges are linearly dependent (to avoid degenerate planes) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + throw(ArgumentError("The vectors `AB` and `AC` must not be collinear")) end # Calculate normal vector of plane @@ -299,7 +231,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) return hcat(c, edge1, edge2) end -@inline function is_in_boundary_zone(boundary_zone::Union{InFlow, OutFlow}, particle_coords) +@inline function is_in_boundary_zone(boundary_zone::BoundaryZone, particle_coords) (; zone_origin, spanning_set) = boundary_zone particle_position = particle_coords - zone_origin @@ -324,7 +256,7 @@ end end function remove_outside_particles(initial_condition, spanning_set, zone_origin) - (; coordinates, density, particle_spacing) = initial_condition + (; coordinates, density, pressure, particle_spacing) = initial_condition in_zone = fill(true, nparticles(initial_condition)) @@ -335,6 +267,6 @@ function remove_outside_particles(initial_condition, spanning_set, zone_origin) in_zone[particle] = is_in_boundary_zone(spanning_set, particle_position) end - return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), - particle_spacing) + return InitialCondition(; coordinates=coordinates[:, in_zone], density=density[in_zone], + pressure=pressure[in_zone], particle_spacing) end diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 38fdd477a..0ac59d6f0 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -6,15 +6,31 @@ This model uses the characteristic variables to propagate the appropriate values to the outlet or inlet and have been proposed by Lastiwka et al. (2009). For more information about the method see [description below](@ref method_of_characteristics). """ -struct BoundaryModelLastiwka end +struct BoundaryModelLastiwka + extrapolate_reference_values::Bool + function BoundaryModelLastiwka(; extrapolate_reference_values::Bool=true) + return new{}(extrapolate_reference_values) + end +end -@inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka, - v, u, v_ode, u_ode, semi, t) - (; density, pressure, cache, flow_direction, +# Called from update callback via update_open_boundary_eachstep! +@inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka, + v, u, v_ode, u_ode, semi, t) + (; density, pressure, cache, boundary_zone, reference_velocity, reference_pressure, reference_density) = system + (; flow_direction) = boundary_zone sound_speed = system_sound_speed(system.fluid_system) + if boundary_model.extrapolate_reference_values + (; prescribed_pressure, prescribed_velocity, prescribed_density) = cache + + @trixi_timeit timer() "extrapolate and correct values" begin + extrapolate_values!(system, v_ode, u_ode, semi, t; prescribed_pressure, + prescribed_velocity, prescribed_density) + end + end + # Update quantities based on the characteristic variables @threaded system for particle in each_moving_particle(system) particle_position = current_coords(u, system, particle) @@ -23,16 +39,16 @@ struct BoundaryModelLastiwka end J2 = cache.characteristics[2, particle] J3 = cache.characteristics[3, particle] - rho_ref = reference_value(reference_density, density[particle], system, particle, + rho_ref = reference_value(reference_density, density[particle], particle_position, t) density[particle] = rho_ref + ((-J1 + 0.5 * (J2 + J3)) / sound_speed^2) - p_ref = reference_value(reference_pressure, pressure[particle], system, particle, + p_ref = reference_value(reference_pressure, pressure[particle], particle_position, t) pressure[particle] = p_ref + 0.5 * (J2 + J3) v_current = current_velocity(v, system, particle) - v_ref = reference_value(reference_velocity, v_current, system, particle, + v_ref = reference_value(reference_velocity, v_current, particle_position, t) rho = density[particle] v_ = v_ref + ((J2 - J3) / (2 * sound_speed * rho)) * flow_direction @@ -45,10 +61,13 @@ struct BoundaryModelLastiwka end return system end +# Called from semidiscretization function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t) @trixi_timeit timer() "evaluate characteristics" begin evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end + + return system end # ==== Characteristics @@ -98,9 +117,16 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end end - characteristics[1, particle] = avg_J1 / counter - characteristics[2, particle] = avg_J2 / counter - characteristics[3, particle] = avg_J3 / counter + # To prevent NANs here if the boundary has not been in contact before. + if counter > 0 + characteristics[1, particle] = avg_J1 / counter + characteristics[2, particle] = avg_J2 / counter + characteristics[3, particle] = avg_J3 / counter + else + characteristics[1, particle] = 0.0 + characteristics[2, particle] = 0.0 + characteristics[3, particle] = 0.0 + end else characteristics[1, particle] /= volume[particle] characteristics[2, particle] /= volume[particle] @@ -114,8 +140,9 @@ end function evaluate_characteristics!(system, neighbor_system::FluidSystem, v, u, v_ode, u_ode, semi, t) - (; volume, cache, flow_direction, density, pressure, + (; volume, cache, boundary_zone, density, pressure, reference_velocity, reference_pressure, reference_density) = system + (; flow_direction) = boundary_zone (; characteristics) = cache v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) @@ -134,15 +161,16 @@ function evaluate_characteristics!(system, neighbor_system::FluidSystem, # Determine current and prescribed quantities rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_ref = reference_value(reference_density, density, system, particle, + rho_ref = reference_value(reference_density, density[particle], neighbor_position, t) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) - p_ref = reference_value(reference_pressure, pressure, system, particle, + p_ref = reference_value(reference_pressure, pressure[particle], neighbor_position, t) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - v_neighbor_ref = reference_value(reference_velocity, v, system, particle, + v_particle = current_velocity(v, system, particle) + v_neighbor_ref = reference_value(reference_velocity, v_particle, neighbor_position, t) # Determine characteristic variables @@ -162,15 +190,15 @@ function evaluate_characteristics!(system, neighbor_system::FluidSystem, return system end -@inline function prescribe_conditions!(characteristics, particle, ::OutFlow) +@inline function prescribe_conditions!(characteristics, particle, ::BoundaryZone{OutFlow}) # J3 is prescribed (i.e. determined from the exterior of the domain). - # J1 and J2 is transimtted from the domain interior. + # J1 and J2 is transmitted from the domain interior. characteristics[3, particle] = zero(eltype(characteristics)) return characteristics end -@inline function prescribe_conditions!(characteristics, particle, ::InFlow) +@inline function prescribe_conditions!(characteristics, particle, ::BoundaryZone{InFlow}) # Allow only J3 to propagate upstream to the boundary characteristics[1, particle] = zero(eltype(characteristics)) characteristics[2, particle] = zero(eltype(characteristics)) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl new file mode 100644 index 000000000..4279fa573 --- /dev/null +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -0,0 +1,198 @@ +# Two ways of providing the information to an open boundary are considered: +# - physical quantities are either assigned a priori +# - or extrapolated from the fluid domain to the buffer zones (inflow and outflow) using ghost nodes + +# The position of the ghost nodes is obtained by mirroring the boundary particles +# into the fluid along a direction that is normal to the open boundary. + +# In order to calculate fluid quantities at the ghost nodes, a standard particle interpolation +# would not be consistent due to the proximity of these points to an open boundary, +# which translates into the lack of a full kernel support. +struct BoundaryModelTafuni end + +function update_quantities!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) + @trixi_timeit timer() "extrapolate and correct values" begin + extrapolate_values!(system, v_ode, u_ode, semi, t; system.cache...) + end +end + +update_final!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system + +function extrapolate_values!(system, v_ode, u_ode, semi, t; prescribed_density=false, + prescribed_pressure=false, prescribed_velocity=false) + (; pressure, density, fluid_system, boundary_zone, reference_density, + reference_velocity, reference_pressure) = system + (; plane_normal) = boundary_zone + + state_equation = system_state_equation(system.fluid_system) + + # Static indices to avoid allocations + two_to_end = SVector{ndims(system)}(2:(ndims(system) + 1)) + + v_open_boundary = wrap_v(v_ode, system, semi) + v_fluid = wrap_v(v_ode, fluid_system, semi) + u_open_boundary = wrap_u(u_ode, system, semi) + u_fluid = wrap_u(u_ode, fluid_system, semi) + + # Use the fluid-fluid nhs, since the boundary particles are mirrored into the fluid domain + neighborhood_search = get_neighborhood_search(fluid_system, fluid_system, semi) + + @threaded system for particle in each_moving_particle(system) + particle_coords = current_coords(u_open_boundary, system, particle) + ghost_node_position = mirror_position(particle_coords, boundary_zone) + + # Set zero + correction_matrix = zero(SMatrix{ndims(system) + 1, ndims(system) + 1, + eltype(system)}) + interpolated_pressure = zero(SVector{ndims(system) + 1, eltype(system)}) + + interpolated_density = zero(SVector{ndims(system) + 1, eltype(system)}) + + interpolated_velocity = zero(SMatrix{ndims(system), ndims(system) + 1, + eltype(system)}) + + for neighbor in PointNeighbors.eachneighbor(ghost_node_position, + neighborhood_search) + neighbor_coords = current_coords(u_fluid, fluid_system, neighbor) + pos_diff = ghost_node_position - neighbor_coords + distance2 = dot(pos_diff, pos_diff) + + if distance2 <= neighborhood_search.search_radius^2 + distance = sqrt(distance2) + + m_b = hydrodynamic_mass(fluid_system, neighbor) + rho_b = particle_density(v_fluid, fluid_system, neighbor) + pressure_b = particle_pressure(v_fluid, fluid_system, neighbor) + v_b = current_velocity(v_fluid, fluid_system, neighbor) + + # Project `v_b` to the normal direction of the boundary zone + # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: + # "Because flow from the inlet interface occurs perpendicular to the boundary, + # only this component of interpolated velocity is kept [...]" + # TODO: Investigate further because it doesn't really make a difference. + #v_b = dot(v_b, plane_normal) * plane_normal + + kernel_value = smoothing_kernel(fluid_system, distance) + grad_kernel = smoothing_kernel_grad(fluid_system, pos_diff, distance) + + L, R = correction_arrays(kernel_value, grad_kernel, pos_diff, rho_b, m_b) + + correction_matrix += L + + !(prescribed_pressure) && (interpolated_pressure += pressure_b * R) + + !(prescribed_density) && (interpolated_density += rho_b * R) + + !(prescribed_velocity) && (interpolated_velocity += v_b * R') + end + end + + if abs(det(correction_matrix)) < 1e-9 + L_inv = I + else + L_inv = inv(correction_matrix) + end + + pos_diff = particle_coords - ghost_node_position + + # In Negi et al. (2020) https://doi.org/10.1016/j.cma.2020.113119, + # they have modified the equation for extrapolation to + # + # f_0 = f_k - (r_0 - r_k) ⋅ ∇f_k + # + # in order to get zero gradient at the outlet interface. + if prescribed_pressure + pressure[particle] = reference_value(reference_pressure, pressure[particle], + particle_coords, t) + else + f_p = L_inv * interpolated_pressure + df_p = f_p[two_to_end] + + pressure[particle] = f_p[1] + dot(pos_diff, df_p) + end + + # TODO: Check if the following is better + # Unlike Tafuni et al. (2018), we calculate the density using the inverse state-equation. + # if prescribed_density + # density[particle] = reference_value(reference_density, density[particle], + # particle_coords, t) + # else + # inverse_state_equation!(density, state_equation, pressure, particle) + # end + if prescribed_density + density[particle] = reference_value(reference_density, density[particle], + particle_coords, t) + else + f_p = L_inv * interpolated_density + df_p = f_p[two_to_end] + + density[particle] = f_p[1] + dot(pos_diff, df_p) + end + + if prescribed_velocity + v_particle = current_velocity(v_open_boundary, system, particle) + v_ref = reference_value(reference_velocity, v_particle, particle_coords, t) + @inbounds for dim in 1:ndims(system) + v_open_boundary[dim, particle] = v_ref[dim] + end + else + @inbounds for dim in 1:ndims(system) + f_v = L_inv * interpolated_velocity[dim, :] + df_v = f_v[two_to_end] + + v_open_boundary[dim, particle] = f_v[1] + dot(pos_diff, df_v) + end + end + end + + return system +end + +function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{3}, rho_b, m_b) + x_ab = pos_diff[1] + y_ab = pos_diff[2] + z_ab = pos_diff[3] + + grad_W_ab_x = grad_W_ab[1] + grad_W_ab_y = grad_W_ab[2] + grad_W_ab_z = grad_W_ab[3] + + V_b = m_b / rho_b + + M_ = @SMatrix [W_ab W_ab*x_ab W_ab*y_ab W_ab*z_ab; + grad_W_ab_x grad_W_ab_x*x_ab grad_W_ab_x*y_ab grad_W_ab_x*z_ab; + grad_W_ab_y grad_W_ab_y*x_ab grad_W_ab_y*y_ab grad_W_ab_y*z_ab; + grad_W_ab_z grad_W_ab_z*x_ab grad_W_ab_z*y_ab grad_W_ab_z*z_ab] + + L = V_b * M_ + + R = V_b * SVector(W_ab, grad_W_ab_x, grad_W_ab_y, grad_W_ab_z) + + return L, R +end + +function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{2}, rho_b, m_b) + x_ab = pos_diff[1] + y_ab = pos_diff[2] + + grad_W_ab_x = grad_W_ab[1] + grad_W_ab_y = grad_W_ab[2] + + V_b = m_b / rho_b + + M_ = @SMatrix [W_ab W_ab*x_ab W_ab*y_ab; + grad_W_ab_x grad_W_ab_x*x_ab grad_W_ab_x*y_ab; + grad_W_ab_y grad_W_ab_y*x_ab grad_W_ab_y*y_ab] + L = V_b * M_ + + R = V_b * SVector(W_ab, grad_W_ab_x, grad_W_ab_y) + + return L, R +end + +function mirror_position(particle_coords, boundary_zone) + particle_position = particle_coords - boundary_zone.zone_origin + dist = dot(particle_position, boundary_zone.plane_normal) + + return particle_coords - 2 * dist * boundary_zone.plane_normal +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 56e7746fc..986c83f4b 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -1,5 +1,5 @@ @doc raw""" - OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; + OpenBoundarySPHSystem(boundary_zone::BoundaryZone; fluid_system::FluidSystem, buffer_size::Integer, boundary_model, reference_velocity=nothing, @@ -9,7 +9,7 @@ Open boundary system for in- and outflow particles. # Arguments -- `boundary_zone`: Use [`InFlow`](@ref) for an inflow and [`OutFlow`](@ref) for an outflow boundary. +- `boundary_zone`: See [`BoundaryZone`](@ref). # Keywords - `fluid_system`: The corresponding fluid system @@ -39,7 +39,6 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] boundary_zone :: BZ - flow_direction :: SVector{NDIMS, ELTYPE} reference_velocity :: RV reference_pressure :: RP reference_density :: RD @@ -47,7 +46,7 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, update_callback_used :: Ref{Bool} cache :: C - function OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; + function OpenBoundarySPHSystem(boundary_zone::BoundaryZone; fluid_system::FluidSystem, buffer_size::Integer, boundary_model, reference_velocity=nothing, @@ -77,6 +76,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "an array where the ``i``-th column holds the velocity of particle ``i`` " * "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity")) else + if reference_velocity isa Function + test_result = reference_velocity(zeros(NDIMS), 0.0) + if length(test_result) != NDIMS + throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS")) + end + end reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) end @@ -86,6 +91,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "each particle's coordinates and time to its pressure, " * "a vector holding the pressure of each particle, or a scalar")) else + if reference_pressure isa Function + test_result = reference_pressure(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_pressure` function must be a scalar function")) + end + end reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) end @@ -95,12 +106,18 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "each particle's coordinates and time to its density, " * "a vector holding the density of each particle, or a scalar")) else + if reference_density isa Function + test_result = reference_density(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_density` function must be a scalar function")) + end + end reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) end - flow_direction_ = boundary_zone.flow_direction - - cache = create_cache_open_boundary(boundary_model, initial_condition) + cache = create_cache_open_boundary(boundary_model, initial_condition, + reference_density, reference_velocity, + reference_pressure) return new{typeof(boundary_model), typeof(boundary_zone), NDIMS, ELTYPE, typeof(initial_condition), typeof(fluid_system), typeof(mass), @@ -108,14 +125,26 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, typeof(reference_density_), typeof(buffer), typeof(cache)}(initial_condition, fluid_system, boundary_model, mass, density, volume, pressure, boundary_zone, - flow_direction_, reference_velocity_, reference_pressure_, + reference_velocity_, reference_pressure_, reference_density_, buffer, false, cache) end end -function create_cache_open_boundary(boundary_model, initial_condition) +function create_cache_open_boundary(boundary_model, initial_condition, + reference_density, reference_velocity, + reference_pressure) ELTYPE = eltype(initial_condition) + prescribed_pressure = isnothing(reference_pressure) ? false : true + prescribed_velocity = isnothing(reference_velocity) ? false : true + prescribed_density = isnothing(reference_density) ? false : true + + if boundary_model isa BoundaryModelTafuni + return (; prescribed_pressure=prescribed_pressure, + prescribed_density=prescribed_density, + prescribed_velocity=prescribed_velocity) + end + characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) @@ -125,12 +154,13 @@ end timer_name(::OpenBoundarySPHSystem) = "open_boundary" vtkname(system::OpenBoundarySPHSystem) = "open_boundary" +boundary_type_name(::BoundaryZone{ZT}) where {ZT} = string(nameof(ZT)) function Base.show(io::IO, system::OpenBoundarySPHSystem) @nospecialize system # reduce precompilation time print(io, "OpenBoundarySPHSystem{", ndims(system), "}(") - print(io, type2string(system.boundary_zone)) + print(io, boundary_type_name(system.boundary_zone)) print(io, ") with ", nparticles(system), " particles") end @@ -145,8 +175,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem) summary_line(io, "#buffer_particles", system.buffer.buffer_size) summary_line(io, "fluid system", type2string(system.fluid_system)) summary_line(io, "boundary model", type2string(system.boundary_model)) - summary_line(io, "boundary", type2string(system.boundary_zone)) - summary_line(io, "flow direction", system.flow_direction) + summary_line(io, "boundary type", boundary_type_name(system.boundary_zone)) summary_line(io, "prescribed velocity", type2string(system.reference_velocity)) summary_line(io, "prescribed pressure", type2string(system.reference_pressure)) summary_line(io, "prescribed density", type2string(system.reference_density)) @@ -173,6 +202,8 @@ end return system.pressure[particle] end +@inline system_smoothing_length(system::OpenBoundarySPHSystem) = system_smoothing_length(system.fluid_system) + function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi, t; update_from_callback=false) if !update_from_callback && !(system.update_callback_used[]) @@ -190,10 +221,12 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ # Update density, pressure and velocity based on the characteristic variables. # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971 - @trixi_timeit timer() "update quantities" update_quantities!(system, - system.boundary_model, - v, u, v_ode, u_ode, - semi, t) + @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system, + system.boundary_model, + v, u, + v_ode, + u_ode, + semi, t) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) @@ -238,7 +271,7 @@ end # Outflow particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::OutFlow, particle, v, u, + boundary_zone::BoundaryZone{OutFlow}, particle, v, u, v_fluid, u_fluid) deactivate_particle!(system, particle, u) @@ -247,7 +280,7 @@ end # Inflow particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::InFlow, particle, v, u, + boundary_zone::BoundaryZone{InFlow}, particle, v, u, v_fluid, u_fluid) (; spanning_set) = boundary_zone @@ -262,6 +295,31 @@ end return system end +# Buffer particle is outside the boundary zone +@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, + boundary_zone::BoundaryZone{BidirectionalFlow}, particle, + v, u, v_fluid, u_fluid) + particle_position = current_coords(u, system, particle) - boundary_zone.zone_origin + + # Check if particle is in- or outside the fluid domain. + # `plane_normal` is always pointing in the fluid domain. + if signbit(dot(particle_position, boundary_zone.plane_normal)) + deactivate_particle!(system, particle, u) + + return system + end + + # Activate a new particle in simulation domain + transfer_particle!(fluid_system, system, particle, v_fluid, u_fluid, v, u) + + # Reset position of boundary particle + for dim in 1:ndims(system) + u[dim, particle] += boundary_zone.spanning_set[1][dim] + end + + return system +end + # Fluid particle is in boundary zone @inline function convert_particle!(fluid_system::FluidSystem, system, boundary_zone, particle, v, u, v_fluid, u_fluid) @@ -341,17 +399,23 @@ function wrap_reference_function(constant_vector_, ::Val{NDIMS}) where {NDIMS} return constant_vector(coords, t) = SVector{NDIMS}(constant_vector_) end -function reference_value(value::Function, quantity, system, particle, position, t) +function reference_value(value::Function, quantity, position, t) return value(position, t) end # This method is used when extrapolating quantities from the domain # instead of using the method of characteristics -reference_value(value::Nothing, quantity, system, particle, position, t) = quantity +reference_value(value::Nothing, quantity, position, t) = quantity + +function check_reference_values!(boundary_model, reference_density, reference_pressure, + reference_velocity) + return boundary_model +end function check_reference_values!(boundary_model::BoundaryModelLastiwka, reference_density, reference_pressure, reference_velocity) - # TODO: Extrapolate the reference values from the domain + boundary_model.extrapolate_reference_values && return boundary_model + if any(isnothing.([reference_density, reference_pressure, reference_velocity])) throw(ArgumentError("for `BoundaryModelLastiwka` all reference values must be specified")) end diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index a3f47f611..b8a836a87 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -21,24 +21,25 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, movement :: M ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) adhesion_coefficient :: ELTYPE + color :: Int64 cache :: CA buffer :: Nothing # This constructor is necessary for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. function BoundarySPHSystem(initial_condition, coordinates, boundary_model, movement, - ismoving, adhesion_coefficient, cache, buffer) + ismoving, adhesion_coefficient, cache, buffer, color) ELTYPE = eltype(coordinates) new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition), typeof(coordinates), typeof(movement), typeof(ismoving), typeof(cache)}(initial_condition, coordinates, boundary_model, movement, - ismoving, adhesion_coefficient, cache, buffer) + ismoving, adhesion_coefficient, color, cache, buffer) end end function BoundarySPHSystem(initial_condition, model; movement=nothing, - adhesion_coefficient=0.0) + adhesion_coefficient=0.0, color_value=0) coordinates = copy(initial_condition.coordinates) ismoving = Ref(!isnothing(movement)) @@ -54,7 +55,36 @@ function BoundarySPHSystem(initial_condition, model; movement=nothing, # Because of dispatches boundary model needs to be first! return BoundarySPHSystem(initial_condition, coordinates, model, movement, - ismoving, adhesion_coefficient, cache, nothing) + ismoving, adhesion_coefficient, cache, nothing, color_value) +end + +function Base.show(io::IO, system::BoundarySPHSystem) + @nospecialize system # reduce precompilation time + + print(io, "BoundarySPHSystem{", ndims(system), "}(") + print(io, system.boundary_model) + print(io, ", ", system.movement) + print(io, ", ", system.adhesion_coefficient) + print(io, ", ", system.color) + print(io, ") with ", nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "BoundarySPHSystem{$(ndims(system))}") + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "boundary model", system.boundary_model) + summary_line(io, "movement function", + isnothing(system.movement) ? "nothing" : + string(system.movement.movement_function)) + summary_line(io, "adhesion coefficient", system.adhesion_coefficient) + summary_line(io, "color", system.color) + summary_footer(io) + end end """ @@ -162,33 +192,6 @@ function create_cache_boundary(::BoundaryMovement, initial_condition) return (; velocity, acceleration, initial_coordinates) end -function Base.show(io::IO, system::BoundarySPHSystem) - @nospecialize system # reduce precompilation time - - print(io, "BoundarySPHSystem{", ndims(system), "}(") - print(io, system.boundary_model) - print(io, ", ", system.movement) - print(io, ", ", system.adhesion_coefficient) - print(io, ") with ", nparticles(system), " particles") -end - -function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) - @nospecialize system # reduce precompilation time - - if get(io, :compact, false) - show(io, system) - else - summary_header(io, "BoundarySPHSystem{$(ndims(system))}") - summary_line(io, "#particles", nparticles(system)) - summary_line(io, "boundary model", system.boundary_model) - summary_line(io, "movement function", - isnothing(system.movement) ? "nothing" : - string(system.movement.movement_function)) - summary_line(io, "adhesion coefficient", system.adhesion_coefficient) - summary_footer(io) - end -end - timer_name(::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary" @inline function Base.eltype(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) @@ -395,6 +398,35 @@ end # viscosity model has to be used. @inline viscosity_model(system::BoundarySPHSystem, neighbor_system::FluidSystem) = neighbor_system.viscosity -function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem) +function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem, semi) return Inf end + +function initialize!(system::BoundarySPHSystem, neighborhood_search) + initialize_colorfield!(system, system.boundary_model, neighborhood_search) + return system +end + +function initialize_colorfield!(system, boundary_model, neighborhood_search) + return system +end + +function initialize_colorfield!(system, ::BoundaryModelDummyParticles, neighborhood_search) + system_coords = system.coordinates + (; smoothing_kernel, smoothing_length) = system.boundary_model + + foreach_point_neighbor(system, system, + system_coords, system_coords, + neighborhood_search, + points=eachparticle(system)) do particle, neighbor, pos_diff, + distance + system.boundary_model.cache.colorfield_bnd[particle] += system.initial_condition.mass[particle] / + system.initial_condition.density[particle] * + system.color * + kernel(smoothing_kernel, + distance, + smoothing_length) + system.boundary_model.cache.neighbor_count[particle] += 1 + end + return system +end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index be8d689f5..dfdc77edb 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -8,6 +8,9 @@ function interact!(dv, v_particle_system, u_particle_system, system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + surface_tension_a = surface_tension_model(particle_system) + surface_tension_b = surface_tension_model(neighbor_system) + # Loop over all pairs of particles and neighbors within the kernel cutoff. foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, @@ -35,8 +38,17 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel) + + dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, + particle, neighbor, pos_diff, distance) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension[i] + + dv_adhesion[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index cf73f5876..5f5db59eb 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -42,12 +42,14 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} + PF, ST, SRFT, SRFN, B, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K smoothing_length :: ELTYPE + number_density :: Int64 + color :: Int64 sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE @@ -55,6 +57,8 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, correction :: Nothing pressure_acceleration_formulation :: PF source_terms :: ST + surface_tension :: SRFT + surface_normal_method :: SRFN buffer :: B cache :: C @@ -65,7 +69,9 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - source_terms=nothing, buffer_size=nothing) + source_terms=nothing, surface_tension=nothing, + surface_normal_method=nothing, buffer_size=nothing, + reference_particle_spacing=0.0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -75,6 +81,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, ELTYPE = eltype(initial_condition) mass = copy(initial_condition.mass) + n_particles = length(initial_condition.mass) if ndims(smoothing_kernel) != NDIMS throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) @@ -85,6 +92,22 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + if surface_tension !== nothing && surface_normal_method === nothing + surface_normal_method = ColorfieldSurfaceNormal(smoothing_kernel, + smoothing_length) + end + + if surface_normal_method !== nothing && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + end + + number_density_ = 0 + if reference_particle_spacing > 0.0 + number_density_ = number_density(Val(NDIMS), reference_particle_spacing, + compact_support(smoothing_kernel, + smoothing_length)) + end + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, @@ -93,14 +116,23 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = create_cache_density(initial_condition, density_calculator) + cache = (; + create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, + n_particles)..., + create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, + n_particles)..., + cache...) new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(pressure_acceleration), typeof(source_terms), + typeof(surface_tension), typeof(surface_normal_method), typeof(buffer), typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, + smoothing_length, number_density_, color_value, sound_speed, + viscosity, nu_edac, acceleration_, nothing, pressure_acceleration, source_terms, + surface_tension, surface_normal_method, buffer, cache) end end @@ -113,6 +145,8 @@ function Base.show(io::IO, system::EntropicallyDampedSPHSystem) print(io, ", ", system.viscosity) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.acceleration) + print(io, ", ", system.surface_tension) + print(io, ", ", system.surface_normal_method) print(io, ") with ", nparticles(system), " particles") end @@ -135,6 +169,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "acceleration", system.acceleration) + summary_line(io, "surface tension", system.surface_tension) + summary_line(io, "surface normal method", system.surface_normal_method) summary_footer(io) end end @@ -160,6 +196,8 @@ end return v[end, particle] end +@inline system_state_equation(system::EntropicallyDampedSPHSystem) = nothing + # WARNING! # These functions are intended to be used internally to set the pressure # of newly activated particles in a callback. @@ -179,6 +217,21 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, compute_density!(system, u, u_ode, semi, system.density_calculator) end +function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) + compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, + t) + compute_surface_delta_function!(system, system.surface_tension) +end + +function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + (; surface_tension) = system + + # Surface normal of neighbor and boundary needs to have been calculated already + compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) +end + function write_v0!(v0, system::EntropicallyDampedSPHSystem, density_calculator) for particle in eachparticle(system) v0[end, particle] = system.initial_condition.pressure[particle] diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 9e2289a8c..a33d45313 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -42,6 +42,10 @@ write_v0!(v0, system, density_calculator) = v0 @inline viscosity_model(system::FluidSystem, neighbor_system::FluidSystem) = neighbor_system.viscosity @inline viscosity_model(system::FluidSystem, neighbor_system::BoundarySystem) = neighbor_system.boundary_model.viscosity +@inline system_state_equation(system::FluidSystem) = system.state_equation + +@inline system_smoothing_length(system::FluidSystem) = system.smoothing_length + function compute_density!(system, u, u_ode, semi, ::ContinuityDensity) # No density update with `ContinuityDensity` return system @@ -54,16 +58,22 @@ function compute_density!(system, u, u_ode, semi, ::SummationDensity) summation_density!(system, semi, u, u_ode, density) end -function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem) - (; smoothing_length, viscosity, acceleration) = system +function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem, semi) + (; smoothing_length, viscosity, acceleration, surface_tension) = system - dt_viscosity = 0.125 * smoothing_length^2 / kinematic_viscosity(system, viscosity) + if !(system.viscosity isa Nothing) + dt_viscosity = 0.125 * smoothing_length^2 / kinematic_viscosity(system, viscosity) + else + dt_viscosity = 0.125 * smoothing_length^2 + end # TODO Adami et al. (2012) just use the gravity here, but Antuono et al. (2012) # are using a per-particle acceleration. Is that supposed to be the previous RHS? + # Morris (2000) also uses the acceleration and cites Monaghan (1992) dt_acceleration = 0.25 * sqrt(smoothing_length / norm(acceleration)) # TODO Everyone seems to be doing this differently. + # Morris (2000) uses the additional condition CFL < 0.25. # Sun et al. (2017) only use h / c (because c depends on v_max as c >= 10 v_max). # Adami et al. (2012) use h / (c + v_max) with a fixed CFL of 0.25. # Antuono et al. (2012) use h / (c + v_max + h * pi_max), where pi is the viscosity coefficient. @@ -72,11 +82,37 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem) # See docstring of the callback for the references. dt_sound_speed = cfl_number * smoothing_length / system_sound_speed(system) - return min(dt_viscosity, dt_acceleration, dt_sound_speed) + # Eq. 28 in Morris (2000) + dt_surface_tension = max(dt_viscosity, dt_acceleration, dt_sound_speed) + if surface_tension isa SurfaceTensionMorris || + surface_tension isa SurfaceTensionMomentumMorris + v = wrap_v(v_ode, system, semi) + dt_surface_tension = sqrt(particle_density(v, system, 1) * smoothing_length^3 / + (2 * pi * surface_tension.surface_tension_coefficient)) + end + + return min(dt_viscosity, dt_acceleration, dt_sound_speed, dt_surface_tension) +end + +@inline function surface_tension_model(system::FluidSystem) + return system.surface_tension +end + +@inline function surface_tension_model(system) + return nothing +end + +@inline function surface_normal_method(system::FluidSystem) + return system.surface_normal_method +end + +@inline function surface_normal_method(system) + return nothing end include("pressure_acceleration.jl") include("viscosity.jl") include("surface_tension.jl") +include("surface_normal_sph.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl new file mode 100644 index 000000000..a55f47e05 --- /dev/null +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -0,0 +1,392 @@ +struct ColorfieldSurfaceNormal{ELTYPE, K} + smoothing_kernel::K + smoothing_length::ELTYPE + boundary_contact_threshold::ELTYPE +end + +function ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length; + boundary_contact_threshold=0.1) + return ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length, + boundary_contact_threshold) +end + +function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) + surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + neighbor_count = Array{ELTYPE, 1}(undef, nparticles) + colorfield = Array{ELTYPE, 1}(undef, nparticles) + return (; surface_normal, neighbor_count) +end + +@inline function surface_normal(particle_system::FluidSystem, particle) + (; cache) = particle_system + return extract_svector(cache.surface_normal, particle_system, particle) +end + +function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, + u_neighbor_system, semi, surfn, nsurfn) + # Normal not needed + return system +end + +# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" +# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_normal!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn, + ::ColorfieldSurfaceNormal) + (; cache) = system + (; smoothing_kernel, smoothing_length) = surfn + + system_coords = current_coordinates(u_system, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + + if smoothing_length != system.smoothing_length || + smoothing_kernel !== system.smoothing_kernel + # TODO: this is really slow but there is no way to easily implement multiple search radia + search_radius = compact_support(smoothing_kernel, smoothing_length) + nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius, + nparticles(system)) + PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords) + end + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + density_neighbor = particle_density(v_neighbor_system, + neighbor_system, neighbor) + grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) + for i in 1:ndims(system) + cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + end + + cache.neighbor_count[particle] += 1 + end + + return system +end + +# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" +# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn, nsurfn) + (; cache) = system + (; colorfield, colorfield_bnd) = neighbor_system.boundary_model.cache + (; smoothing_kernel, smoothing_length, boundary_contact_threshold) = surfn + + system_coords = current_coordinates(u_system, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + + if smoothing_length != system.smoothing_length || + smoothing_kernel !== system.smoothing_kernel + # TODO: this is really slow but there is no way to easily implement multiple search radia + search_radius = compact_support(smoothing_kernel, smoothing_length) + nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius, + nparticles(system)) + nhs_bnd = PointNeighbors.copy_neighborhood_search(nhs_bnd, search_radius, + nparticles(neighbor_system)) + PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords) + PointNeighbors.initialize!(nhs_bnd, neighbor_system_coords, neighbor_system_coords) + end + + # First we need to calculate the smoothed colorfield values of the boundary + # TODO: move colorfield to extra step + # TODO: this is only correct for a single fluid + + # Reset to the constant boundary interpolated color values + colorfield .= colorfield_bnd + + # Accumulate fluid neighbors + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + colorfield[neighbor] += hydrodynamic_mass(system, particle) / + particle_density(v, system, particle) * system.color * + kernel(smoothing_kernel, distance, smoothing_length) + end + + maximum_colorfield = maximum(colorfield) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + # we assume that we are in contact with the boundary if the color of the boundary particle is larger than the threshold + if colorfield[neighbor] / maximum_colorfield > boundary_contact_threshold + m_b = hydrodynamic_mass(system, particle) + density_neighbor = particle_density(v, system, particle) + grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length) + for i in 1:ndims(system) + cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + end + cache.neighbor_count[particle] += 1 + end + end + + return system +end + +function remove_invalid_normals!(system, surface_tension) + # Normal not needed + return system +end + +function remove_invalid_normals!(system::FluidSystem, surface_tension::SurfaceTensionAkinci) + (; cache) = system + + # We remove invalid normals (too few neighbors) to reduce the impact of underdefined normals + for particle in each_moving_particle(system) + # A corner has that many neighbors assuming a regular 2 * r distribution and a compact_support of 4r + if cache.neighbor_count[particle] < 2^ndims(system) + 1 + cache.surface_normal[1:ndims(system), particle] .= 0 + end + end + + return system +end + +# see Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function remove_invalid_normals!(system::FluidSystem, + surface_tension::Union{SurfaceTensionMorris, + SurfaceTensionMomentumMorris}) + (; cache, smoothing_length, smoothing_kernel, number_density) = system + (; free_surface_threshold) = surface_tension + + # println("compact_support ", compact_support(smoothing_kernel, smoothing_length)) + + # TODO: make settable + # We remove invalid normals i.e. they have a small norm (eq. 20) + # normal_condition2 = (0.01 / compact_support(smoothing_kernel, smoothing_length))^2 + normal_condition2 = (0.1 / compact_support(smoothing_kernel, smoothing_length))^2 + + for particle in each_moving_particle(system) + + # TODO: make selectable + # heuristic condition if there is no gas phase to find the free surface + # if free_surface_threshold * number_density < cache.neighbor_count[particle] + # cache.surface_normal[1:ndims(system), particle] .= 0 + # continue + # end + + particle_surface_normal = cache.surface_normal[1:ndims(system), particle] + norm2 = dot(particle_surface_normal, particle_surface_normal) + + # println(norm2, " > ", normal_condition2) + # see eq. 21 + if norm2 > normal_condition2 + cache.surface_normal[1:ndims(system), particle] = particle_surface_normal / + sqrt(norm2) + else + cache.surface_normal[1:ndims(system), particle] .= 0 + end + end + + # println("after removable: ") + # println(cache.surface_normal) + + return system +end + +function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + return system +end + +function compute_surface_normal!(system::FluidSystem, + surface_normal_method_::ColorfieldSurfaceNormal, + v, u, v_ode, u_ode, semi, t) + (; cache, surface_tension) = system + + # Reset surface normal + set_zero!(cache.surface_normal) + set_zero!(cache.neighbor_count) + + # TODO: if color values are set only different systems need to be called + # TODO: what to do if there is no gas phase? -> config values + @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + calc_normal!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method_, + surface_normal_method(neighbor_system)) + end + remove_invalid_normals!(system, surface_tension) + + return system +end + +function calc_curvature!(system, neighbor_system, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn, nsurfn) +end + +# Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_curvature!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, + surfn::ColorfieldSurfaceNormal, nsurfn::ColorfieldSurfaceNormal) + (; cache) = system + (; smoothing_kernel, smoothing_length) = surfn + (; curvature) = cache + + system_coords = current_coordinates(u_system, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + correction_factor = fill(eps(eltype(system)), n_moving_particles(system)) + + if smoothing_length != system.smoothing_length || + smoothing_kernel !== system.smoothing_kernel + # TODO: this is really slow but there is no way to easily implement multiple search radia + search_radius = compact_support(smoothing_kernel, smoothing_length) + nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius, + nparticles(system)) + PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords) + end + + no_valid_neighbors = 0 + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + n_a = surface_normal(system, particle) + n_b = surface_normal(neighbor_system, neighbor) + grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) + v_b = m_b / rho_b + + # eq. 22 we can test against eps() here since the surface normals that are invalid have been reset + if dot(n_a, n_a) > eps() && dot(n_b, n_b) > eps() + # for i in 1:ndims(system) + curvature[particle] += v_b * dot((n_b .- n_a), grad_kernel) + # end + # eq. 24 + correction_factor[particle] += v_b * kernel(smoothing_kernel, distance, + smoothing_length) + # prevent NaNs from systems that are entirely skipped + no_valid_neighbors += 1 + end + end + + # eq. 23 + if no_valid_neighbors > 0 + for i in 1:n_moving_particles(system) + curvature[i] /= correction_factor[i] + end + end + + # println("after curvature") + # println("surf_norm ", cache.surface_normal) + # println("curv ", cache.curvature) + # println("C ", correction_factor) + + return system +end + +function compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + return system +end + +function compute_curvature!(system::FluidSystem, surface_tension::SurfaceTensionMorris, v, + u, v_ode, u_ode, semi, t) + (; cache, surface_tension) = system + + # Reset surface curvature + set_zero!(cache.curvature) + + @trixi_timeit timer() "compute surface curvature" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + calc_curvature!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method(system), + surface_normal_method(neighbor_system)) + end + return system +end + +function calc_stress_tensors!(system::FluidSystem, neighbor_system::FluidSystem, u_system, + v, + v_neighbor_system, u_neighbor_system, semi, + surfn::ColorfieldSurfaceNormal, + nsurfn::ColorfieldSurfaceNormal) + (; cache) = system + (; smoothing_kernel, smoothing_length) = surfn + (; stress_tensor, delta_s) = cache + + neighbor_cache = neighbor_system.cache + neighbor_delta_s = neighbor_cache.delta_s + + NDIMS = ndims(system) + max_delta_s = max(maximum(delta_s), maximum(neighbor_delta_s)) + + for particle in each_moving_particle(system) + normal = surface_normal(system, particle) + delta_s_particle = delta_s[particle] + if delta_s_particle > eps() + for i in 1:NDIMS + for j in 1:NDIMS + delta_ij = (i == j) ? 1.0 : 0.0 + stress_tensor[i, j, particle] = delta_s_particle * + (delta_ij - normal[i] * normal[j]) - + delta_ij * max_delta_s + end + end + end + end + + return system +end + +function compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + return system +end + +# Section 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function compute_stress_tensors!(system::FluidSystem, ::SurfaceTensionMomentumMorris, + v, u, v_ode, u_ode, semi, t) + (; cache) = system + (; delta_s, stress_tensor) = cache + + # Reset surface stress_tensor + set_zero!(stress_tensor) + + max_delta_s = maximum(delta_s) + NDIMS = ndims(system) + + @trixi_timeit timer() "compute surface stress tensor" for particle in each_moving_particle(system) + normal = surface_normal(system, particle) + delta_s_particle = delta_s[particle] + if delta_s_particle > eps() + for i in 1:NDIMS + for j in 1:NDIMS + delta_ij = (i == j) ? 1.0 : 0.0 + stress_tensor[i, j, particle] = delta_s_particle * + (delta_ij - normal[i] * normal[j]) - + delta_ij * max_delta_s + end + end + end + end + + return system +end + +function compute_surface_delta_function!(system, surface_tension) + return system +end + +# eq. 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function compute_surface_delta_function!(system, ::SurfaceTensionMomentumMorris) + (; cache) = system + (; delta_s) = cache + + set_zero!(delta_s) + + for particle in each_moving_particle(system) + delta_s[particle] = norm(surface_normal(system, particle)) + end + return system +end diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 66140775e..35ed7084e 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -1,4 +1,5 @@ -abstract type AkinciTypeSurfaceTension end +abstract type SurfaceTensionModel end +abstract type AkinciTypeSurfaceTension <: SurfaceTensionModel end @doc raw""" CohesionForceAkinci(surface_tension_coefficient=1.0) @@ -38,6 +39,45 @@ struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension end end +struct SurfaceTensionMorris{ELTYPE} <: SurfaceTensionModel + surface_tension_coefficient::ELTYPE + free_surface_threshold::ELTYPE + + function SurfaceTensionMorris(; surface_tension_coefficient=1.0, + free_surface_threshold=0.75) + new{typeof(surface_tension_coefficient)}(surface_tension_coefficient, + free_surface_threshold) + end +end + +function create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_surface_tension(::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles) + curvature = Array{ELTYPE, 1}(undef, nparticles) + return (; curvature) +end + +struct SurfaceTensionMomentumMorris{ELTYPE} <: SurfaceTensionModel + surface_tension_coefficient::ELTYPE + free_surface_threshold::ELTYPE + + function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0, + free_surface_threshold=0.75) + new{typeof(surface_tension_coefficient)}(surface_tension_coefficient, + free_surface_threshold) + end +end + +function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, + nparticles) + # Allocate stress tensor for each particle: NDIMS x NDIMS x nparticles + delta_s = Array{ELTYPE, 1}(undef, nparticles) + stress_tensor = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, nparticles) + return (; stress_tensor, delta_s) +end + # 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 @@ -89,47 +129,18 @@ end return adhesion_force end -# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" -# Note: most of the time this only leads to an approximation of the surface normal - -function calc_normal_akinci!(system, neighbor_system::FluidSystem, - surface_tension::SurfaceTensionAkinci, u_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search) - (; smoothing_length, cache) = system - - system_coords = current_coordinates(u_system, system) - neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, - system_coords, neighbor_system_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - density_neighbor = particle_density(v_neighbor_system, - neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, - particle) - for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * - grad_kernel[i] * smoothing_length - end - end - - return system -end - -function calc_normal_akinci!(system, neighbor_system, surface_tension, u_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search) - # Normal not needed - return system +# Skip +@inline function surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel) + return zero(pos_diff) end @inline function surface_tension_force(surface_tension_a::CohesionForceAkinci, surface_tension_b::CohesionForceAkinci, particle_system::FluidSystem, neighbor_system::FluidSystem, particle, neighbor, - pos_diff, distance) + pos_diff, distance, rho_a, rho_b, grad_kernel) (; smoothing_length) = particle_system # No cohesion with oneself distance < sqrt(eps()) && return zero(pos_diff) @@ -144,7 +155,7 @@ end surface_tension_b::SurfaceTensionAkinci, particle_system::FluidSystem, neighbor_system::FluidSystem, particle, neighbor, - pos_diff, distance) + pos_diff, distance, rho_a, rho_b, grad_kernel) (; smoothing_length, smoothing_kernel) = particle_system (; surface_tension_coefficient) = surface_tension_a @@ -152,20 +163,47 @@ end distance < sqrt(eps()) && return zero(pos_diff) m_b = hydrodynamic_mass(neighbor_system, neighbor) - n_a = surface_normal(surface_tension_a, particle_system, particle) - n_b = surface_normal(surface_tension_b, neighbor_system, neighbor) + n_a = surface_normal(particle_system, particle) + n_b = surface_normal(neighbor_system, neighbor) support_radius = compact_support(smoothing_kernel, smoothing_length) return cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) .- - (surface_tension_coefficient * (n_a - n_b)) + (surface_tension_coefficient * (n_a - n_b) * smoothing_length) end -# Skip -@inline function surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, particle, neighbor, - pos_diff, distance) - return zero(pos_diff) +@inline function surface_tension_force(surface_tension_a::SurfaceTensionMorris, + surface_tension_b::SurfaceTensionMorris, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel) + (; surface_tension_coefficient) = surface_tension_a + + # No surface tension with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + n_a = surface_normal(particle_system, particle) + curvature_a = curvature(particle_system, particle) + + return -surface_tension_coefficient / rho_a * curvature_a * n_a +end + +@inline function surface_tension_force(surface_tension_a::SurfaceTensionMomentumMorris, + surface_tension_b::SurfaceTensionMomentumMorris, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel) + (; surface_tension_coefficient) = surface_tension_a + + # No surface tension with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + S_a = particle_system.cache.stress_tensor[:, :, particle] + S_b = neighbor_system.cache.stress_tensor[:, :, neighbor] + + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + return surface_tension_coefficient * m_b * (S_a + S_b) / (rho_a * rho_b) * grad_kernel end @inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 09bde77f9..61fc72e6a 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -167,7 +167,7 @@ end # In: Reports on Progress in Physics (2005), pages 1703-1759. # [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01) function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan) - (; smoothing_length) = system + smoothing_length = system_smoothing_length(system) (; alpha) = viscosity sound_speed = system_sound_speed(system) diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c4ddc8b81..4eeebcc2c 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -208,7 +208,7 @@ end (; delta) = density_diffusion (; smoothing_length, state_equation) = particle_system - (; sound_speed) = state_equation + sound_speed_ = sound_speed(state_equation) volume_b = m_b / rho_b @@ -216,7 +216,7 @@ end particle_system, particle, neighbor) density_diffusion_term = dot(psi, grad_kernel) * volume_b - dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term + dv[end, particle] += delta * smoothing_length * sound_speed_ * density_diffusion_term end # Density diffusion `nothing` or interaction other than fluid-fluid diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 9e40eacfa..4330c96a7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -6,8 +6,9 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::WeaklyCompressibleSPHSystem, neighbor_system) - (; density_calculator, state_equation, correction, surface_tension) = particle_system - (; sound_speed) = state_equation + (; density_calculator, state_equation, correction) = particle_system + + sound_speed_ = sound_speed(state_equation) surface_tension_a = surface_tension_model(particle_system) surface_tension_b = surface_tension_model(neighbor_system) @@ -47,23 +48,35 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, neighbor_system, particle, neighbor) - dv_pressure = pressure_correction * + n = 4.0 + alpha = 0.01 + q = distance/particle_system.reference_particle_spacing + if q < 1.0 + Pi_ij = -alpha * (1 - q)^n / rho_mean + else + Pi_ij = 0.0 + end + + dv_pressure = pressure_correction * ( pressure_acceleration(particle_system, neighbor_system, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) + distance, grad_kernel, correction) -m_b * Pi_ij * grad_kernel) + + dv_viscosity_ = viscosity_correction * dv_viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + sound_speed_, m_a, m_b, rho_a, rho_b, grad_kernel) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, - particle, neighbor, pos_diff, distance) + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel) - dv_adhesion = adhesion_force(surface_tension, particle_system, neighbor_system, + dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) @inbounds for i in 1:ndims(particle_system) diff --git a/src/schemes/fluid/weakly_compressible_sph/state_equations.jl b/src/schemes/fluid/weakly_compressible_sph/state_equations.jl index ea0f7d2ef..0a47f1f4b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/state_equations.jl +++ b/src/schemes/fluid/weakly_compressible_sph/state_equations.jl @@ -49,3 +49,65 @@ function inverse_state_equation(state_equation::StateEquationCole, pressure) return reference_density * tmp^(1 / exponent) end + +@inline sound_speed(eos) = eos.sound_speed + +@doc raw""" + StateEquationIdealGas(; gas_constant, temperature, gamma) +Equation of state to describe the relationship between pressure and density +of a gas using the Ideal Gas Law. +# Keywords +- `gas_constant`: Specific gas constant (R) for the gas, typically in J/(kg*K). +- `temperature` : Absolute temperature of the gas in Kelvin. +- `gamma` : Heat-capacity ratio +This struct calculates the pressure of a gas from its density using the formula: +\[ P = \rho \cdot R \cdot T \] +where \( P \) is pressure, \( \rho \) is density, \( R \) is the gas constant, and \( T \) is temperature. +Note: +For basic WCSPH this boils down to the assumption of a linear pressure-density relationship. +""" +struct StateEquationIdealGas{ELTYPE} + gas_constant :: ELTYPE + temperature :: ELTYPE + gamma :: ELTYPE + + function StateEquationIdealGas(; gas_constant, temperature, gamma) + new{typeof(gas_constant)}(gas_constant, temperature, gamma) + end +end + +function (state_equation::StateEquationIdealGas)(density) + (; gas_constant, temperature) = state_equation + pressure = density * gas_constant * temperature + return pressure +end + +# This version is for simulations that include a temperature. +function (state_equation::StateEquationIdealGas)(density, internal_energy) + (; gamma) = state_equation + pressure = (gamma - 1.0) * density * internal_energy + return pressure +end + +function inverse_state_equation(state_equation::StateEquationIdealGas, pressure) + (; gas_constant, temperature) = state_equation + density = pressure / (gas_constant * temperature) + return density +end + +# This version is for simulations that include a temperature. +function inverse_state_equation(state_equation::StateEquationIdealGas, pressure, + internal_energy) + gamma = state_equation.gamma + + density = pressure / ((gamma - 1.0) * internal_energy) + return density +end + +@inline sound_speed(eos::StateEquationIdealGas) = sqrt(eos.gamma * eos.gas_constant * + eos.temperature) +# This version is for simulations that include a temperature. +@inline sound_speed(eos::StateEquationIdealGas, pressure, density) = sqrt(eos.gamma * + pressure / + (density * + eos.gas_constant)) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 482261fc1..846069547 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -45,7 +45,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} + V, DD, COR, PF, ST, B, SRFT, SRFN, C} <: + FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -53,6 +54,9 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, state_equation :: SE smoothing_kernel :: K smoothing_length :: ELTYPE + reference_particle_spacing :: ELTYPE + number_density :: Int64 + color :: Int64 acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V density_diffusion :: DD @@ -60,6 +64,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, pressure_acceleration_formulation :: PF source_terms :: ST surface_tension :: SRFT + surface_normal_method :: SRFN buffer :: B cache :: C end @@ -75,7 +80,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, - surface_tension=nothing) + surface_tension=nothing, surface_normal_method=nothing, + reference_particle_spacing=0.0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -103,6 +109,21 @@ function WeaklyCompressibleSPHSystem(initial_condition, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end + if surface_tension !== nothing && surface_normal_method === nothing + surface_normal_method = ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length) + end + + if surface_normal_method !== nothing && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + end + + number_density_ = 0 + if reference_particle_spacing > 0.0 + number_density_ = number_density(Val(NDIMS), reference_particle_spacing, + compact_support(smoothing_kernel, + smoothing_length)) + end + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, @@ -113,16 +134,21 @@ function WeaklyCompressibleSPHSystem(initial_condition, create_cache_wcsph(correction, initial_condition.density, NDIMS, n_particles)..., cache...) cache = (; - create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., + create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, + n_particles)..., + create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, + n_particles)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, + smoothing_kernel, smoothing_length, reference_particle_spacing, + number_density_, color_value, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, - source_terms, surface_tension, buffer, cache) + source_terms, surface_tension, surface_normal_method, + buffer, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) @@ -149,11 +175,6 @@ function create_cache_wcsph(::MixedKernelGradientCorrection, density, NDIMS, n_p return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end -function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) - surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) - return (; surface_normal) -end - function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time @@ -165,6 +186,10 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.viscosity) print(io, ", ", system.density_diffusion) print(io, ", ", system.surface_tension) + print(io, ", ", system.surface_normal_method) + if system.surface_normal_method isa ColorfieldSurfaceNormal + print(io, ", ", system.color) + end print(io, ", ", system.acceleration) print(io, ", ", system.source_terms) print(io, ") with ", nparticles(system), " particles") @@ -192,6 +217,10 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "viscosity", system.viscosity) summary_line(io, "density diffusion", system.density_diffusion) summary_line(io, "surface tension", system.surface_tension) + summary_line(io, "surface normal method", system.surface_normal_method) + if system.surface_normal_method isa ColorfieldSurfaceNormal + summary_line(io, "color", system.color) + end summary_line(io, "acceleration", system.acceleration) summary_line(io, "source terms", system.source_terms |> typeof |> nameof) summary_footer(io) @@ -214,7 +243,10 @@ end return system.pressure[particle] end -@inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed +@inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = sound_speed(system.state_equation) +@inline system_sound_speed(system::WeaklyCompressibleSPHSystem, pressure, density) = sound_speed(system.state_equation, + pressure, + density) function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) @@ -230,7 +262,7 @@ function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction, surface_tension) = system + (; density_calculator, correction, surface_normal_method, surface_tension) = system compute_correction_values!(system, correction, u, v_ode, u_ode, semi) @@ -240,11 +272,20 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) compute_pressure!(system, v) - compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - + compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + compute_surface_delta_function!(system, surface_tension) return system end +function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + (; surface_tension) = system + + # Surface normal of neighbor and boundary needs to have been calculated already + compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) +end + function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, correction, density_calculator) return system @@ -356,38 +397,7 @@ end extract_smatrix(system.cache.correction_matrix, system, particle) end -function compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - return system -end - -function compute_surface_normal!(system, surface_tension::SurfaceTensionAkinci, v, u, v_ode, - u_ode, semi, t) - (; cache) = system - - # Reset surface normal - set_zero!(cache.surface_normal) - - @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system - u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) - v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) - nhs = get_neighborhood_search(system, semi) - - calc_normal_akinci!(system, neighbor_system, surface_tension, u, v_neighbor_system, - u_neighbor_system, nhs) - end - return system -end - -@inline function surface_normal(::SurfaceTensionAkinci, particle_system::FluidSystem, - particle) +@inline function curvature(particle_system::FluidSystem, particle) (; cache) = particle_system - return extract_svector(cache.surface_normal, particle_system, particle) -end - -@inline function surface_tension_model(system::WeaklyCompressibleSPHSystem) - return system.surface_tension -end - -@inline function surface_tension_model(system) - return nothing + return cache.curvature[particle] end diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4da0f7783..0f5754ce9 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -157,7 +157,7 @@ end function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::TotalLagrangianSPHSystem, - neighbor_system::BoundarySPHSystem) + neighbor_system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}) # TODO continuity equation? return dv end diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index c666ec3d7..0846bbb25 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -177,8 +177,8 @@ function sample_plane(plane_points::NTuple{3}, particle_spacing; tlsph=nothing) edge2 = point3_ - point1_ # Check if the points are collinear - if norm(cross(edge1, edge2)) == 0 - throw(ArgumentError("the points must not be collinear")) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + throw(ArgumentError("The vectors `AB` and `AC` must not be collinear")) end # Determine the number of points along each edge diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 63d7628f4..f3e6f522e 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -247,6 +247,48 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["pressure"] = [particle_pressure(v, system, particle) for particle in active_particles(system)] + if system.surface_normal_method !== nothing + vtk["surf_normal"] = [surface_normal(system, particle) + for particle in eachparticle(system)] + vtk["neighbor_count"] = system.cache.neighbor_count + vtk["color"] = system.color + end + + if system.surface_tension isa SurfaceTensionMorris || + system.surface_tension isa SurfaceTensionMomentumMorris + surft = zeros((ndims(system), n_moving_particles(system))) + system_coords = current_coordinates(u, system) + + surface_tension_a = surface_tension_model(system) + surface_tension_b = surface_tension_model(system) + nhs = create_neighborhood_search(nothing, system, system) + + foreach_point_neighbor(system, system, + system_coords, system_coords, + nhs) do particle, neighbor, pos_diff, + distance + rho_a = particle_density(v, system, particle) + rho_b = particle_density(v, system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + surft[1:ndims(system), particle] .+= surface_tension_force(surface_tension_a, + surface_tension_b, + system, system, + particle, neighbor, + pos_diff, distance, + rho_a, rho_b, + grad_kernel) + end + vtk["surface_tension"] = surft + + if system.surface_tension isa SurfaceTensionMorris + vtk["curvature"] = system.cache.curvature + end + if system.surface_tension isa SurfaceTensionMomentumMorris + vtk["surface_stress_tensor"] = system.cache.stress_tensor + end + end + if write_meta_data vtk["acceleration"] = system.acceleration vtk["viscosity"] = type2string(system.viscosity) @@ -256,17 +298,31 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["density_calculator"] = type2string(system.density_calculator) if system isa WeaklyCompressibleSPHSystem + vtk["solver"] = "WCSPH" + vtk["correction_method"] = type2string(system.correction) + vtk["number_density"] = system.number_density if system.correction isa AkinciFreeSurfaceCorrection vtk["correction_rho0"] = system.correction.rho0 end vtk["state_equation"] = type2string(system.state_equation) - vtk["state_equation_rho0"] = system.state_equation.reference_density - vtk["state_equation_pa"] = system.state_equation.background_pressure - vtk["state_equation_c"] = system.state_equation.sound_speed - vtk["state_equation_exponent"] = system.state_equation.exponent + if system.state_equation isa StateEquationCole + vtk["state_equation_rho0"] = system.state_equation.reference_density + vtk["state_equation_pa"] = system.state_equation.background_pressure + vtk["state_equation_c"] = system.state_equation.sound_speed + vtk["state_equation_exponent"] = system.state_equation.exponent + end - vtk["solver"] = "WCSPH" + if system.state_equation isa StateEquationIdealGas + vtk["state_equation_gamma"] = system.state_equation.gamma + vtk["state_equation_gas_c"] = system.state_equation.gas_constant + vtk["state_equation_temperature"] = system.state_equation.temperature + end + + if system isa StateEquationIdealGas + vtk["state_equation_temperature"] = system.state_equation.temperature + vtk["state_equation_gas_constant"] = system.state_equation.gas_constant + end else vtk["solver"] = "EDAC" vtk["sound_speed"] = system.sound_speed @@ -327,14 +383,15 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data for particle in active_particles(system)] vtk["pressure"] = [particle_pressure(v, system, particle) for particle in active_particles(system)] + vtk["colorfield"] = system.cache.colorfield if write_meta_data - vtk["boundary_zone"] = type2string(system.boundary_zone) + vtk["boundary_zone"] = type2string(first(typeof(system.boundary_zone).parameters)) vtk["width"] = round(system.boundary_zone.zone_width, digits=3) - vtk["flow_direction"] = system.flow_direction vtk["velocity_function"] = type2string(system.reference_velocity) vtk["pressure_function"] = type2string(system.reference_pressure) vtk["density_function"] = type2string(system.reference_density) + vtk["color"] = system.color end return vtk @@ -361,16 +418,19 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system; write_meta_data=true) if write_meta_data vtk["boundary_model"] = "BoundaryModelDummyParticles" - vtk["smoothing_kernel"] = type2string(system.boundary_model.smoothing_kernel) - vtk["smoothing_length"] = system.boundary_model.smoothing_length - vtk["density_calculator"] = type2string(system.boundary_model.density_calculator) - vtk["state_equation"] = type2string(system.boundary_model.state_equation) + vtk["smoothing_kernel"] = type2string(model.smoothing_kernel) + vtk["smoothing_length"] = model.smoothing_length + vtk["density_calculator"] = type2string(model.density_calculator) + vtk["state_equation"] = type2string(model.state_equation) vtk["viscosity_model"] = type2string(model.viscosity) end vtk["hydrodynamic_density"] = [particle_density(v, system, particle) for particle in eachparticle(system)] vtk["pressure"] = model.pressure + vtk["colorfield_bnd"] = model.cache.colorfield_bnd + vtk["colorfield"] = model.cache.colorfield + vtk["neighbor_count"] = model.cache.neighbor_count if model.viscosity isa ViscosityAdami vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :) diff --git a/test/examples/examples.jl b/test/examples/examples.jl index c962b5912..256fdc29d 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -191,6 +191,18 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/dam_break_2phase_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2phase_2d.jl"), + tspan=(0.0, 0.05)) [ + 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_3d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", @@ -226,6 +238,14 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/pipe_flow_3d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), + joinpath(examples_dir(), "fluid", + "pipe_flow_3d.jl")) + @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", diff --git a/test/general/buffer.jl b/test/general/buffer.jl index 841945784..a188d5fce 100644 --- a/test/general/buffer.jl +++ b/test/general/buffer.jl @@ -2,13 +2,14 @@ # Mock fluid system struct FluidSystemMock3 <: TrixiParticles.FluidSystem{2, Nothing} end - inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2, - open_boundary_layers=2, density=1.0, flow_direction=[1.0, 0.0]) - system = OpenBoundarySPHSystem(inflow; fluid_system=FluidSystemMock3(), + zone = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2, + open_boundary_layers=2, density=1.0, plane_normal=[1.0, 0.0], + boundary_type=:inflow) + system = OpenBoundarySPHSystem(zone; fluid_system=FluidSystemMock3(), reference_density=0.0, reference_pressure=0.0, reference_velocity=[0, 0], boundary_model=BoundaryModelLastiwka(), buffer_size=0) - system_buffer = OpenBoundarySPHSystem(inflow; buffer_size=5, + system_buffer = OpenBoundarySPHSystem(zone; buffer_size=5, reference_density=0.0, reference_pressure=0.0, reference_velocity=[0, 0], boundary_model=BoundaryModelLastiwka(), diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index 7a34a4bca..0bbc9a9f5 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -18,22 +18,23 @@ ] @testset verbose=true "Flow Direction $j" for j in eachindex(flow_directions) - inflow = InFlow(; plane=(point_1, point_2), particle_spacing, - flow_direction=flow_directions[j], density=1.0, - open_boundary_layers) - outflow = OutFlow(; plane=(point_1, point_2), particle_spacing, - flow_direction=flow_directions[j], density=1.0, - open_boundary_layers) + inflow = BoundaryZone(; plane=(point_1, point_2), particle_spacing, + plane_normal=flow_directions[j], density=1.0, + open_boundary_layers, boundary_type=:inflow) + outflow = BoundaryZone(; plane=(point_1, point_2), particle_spacing, + plane_normal=-flow_directions[j], density=1.0, + open_boundary_layers, boundary_type=:outflow) boundary_zones = [ inflow, outflow, ] - @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing - sign_ = (boundary_zone isa InFlow) ? -1 : 1 + sign_ = (first(typeof(boundary_zone).parameters) === + TrixiParticles.InFlow) ? -1 : 1 @test plane_points_1[i] == boundary_zone.zone_origin @test plane_points_2[i] - boundary_zone.zone_origin == @@ -78,22 +79,24 @@ ] @testset verbose=true "Flow Direction $j" for j in eachindex(flow_directions) - inflow = InFlow(; plane=(point_1, point_2, point_3), particle_spacing, - flow_direction=flow_directions[j], density=1.0, - open_boundary_layers) - outflow = OutFlow(; plane=(point_1, point_2, point_3), particle_spacing, - flow_direction=flow_directions[j], density=1.0, - open_boundary_layers) + inflow = BoundaryZone(; plane=(point_1, point_2, point_3), particle_spacing, + plane_normal=flow_directions[j], density=1.0, + open_boundary_layers, boundary_type=:inflow) + outflow = BoundaryZone(; plane=(point_1, point_2, point_3), + particle_spacing, + plane_normal=-flow_directions[j], density=1.0, + open_boundary_layers, boundary_type=:outflow) boundary_zones = [ inflow, outflow, ] - @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing - sign_ = (boundary_zone isa InFlow) ? -1 : 1 + sign_ = (first(typeof(boundary_zone).parameters) === + TrixiParticles.InFlow) ? -1 : 1 @test plane_points_1[i] == boundary_zone.zone_origin @test plane_points_2[i] - boundary_zone.zone_origin == @@ -115,18 +118,22 @@ flow_direction = normalize([-plane_size[2], plane_size[1]]) - inflow = InFlow(; plane=plane_points, particle_spacing=0.1, - flow_direction, density=1.0, open_boundary_layers=4) - outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, - flow_direction, density=1.0, open_boundary_layers=4) + inflow = BoundaryZone(; plane=plane_points, particle_spacing=0.1, + plane_normal=flow_direction, density=1.0, + open_boundary_layers=4, boundary_type=:inflow) + outflow = BoundaryZone(; plane=plane_points, particle_spacing=0.1, + plane_normal=-flow_direction, density=1.0, + open_boundary_layers=4, boundary_type=:outflow) boundary_zones = [ inflow, outflow, ] - @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones - perturb_ = boundary_zone isa InFlow ? sqrt(eps()) : -sqrt(eps()) + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones + perturb_ = first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow ? + sqrt(eps()) : + -sqrt(eps()) point1 = plane_points[1] point2 = plane_points[2] @@ -155,18 +162,21 @@ flow_direction = normalize(cross(point2 - point1, point3 - point1)) - inflow = InFlow(; plane=[point1, point2, point3], particle_spacing=0.1, - flow_direction, density=1.0, open_boundary_layers=4) - outflow = OutFlow(; plane=[point1, point2, point3], particle_spacing=0.1, - flow_direction, density=1.0, open_boundary_layers=4) + inflow = BoundaryZone(; plane=[point1, point2, point3], particle_spacing=0.1, + plane_normal=flow_direction, density=1.0, + open_boundary_layers=4, boundary_type=:inflow) + outflow = BoundaryZone(; plane=[point1, point2, point3], particle_spacing=0.1, + plane_normal=-flow_direction, density=1.0, + open_boundary_layers=4, boundary_type=:outflow) boundary_zones = [ inflow, outflow, ] - @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones - perturb_ = boundary_zone isa InFlow ? eps() : -eps() + @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones + perturb_ = first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow ? + eps() : -eps() point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin query_points = Dict( @@ -187,35 +197,43 @@ end @testset verbose=true "Illegal Inputs" begin - no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] + no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.4, 0.9, -0.15]] flow_direction = [0.0, 0.0, 1.0] - error_str = "the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal" - - @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane, - particle_spacing=0.1, - flow_direction, density=1.0, - open_boundary_layers=2) - @test_throws ArgumentError(error_str) OutFlow(; plane=no_rectangular_plane, - particle_spacing=0.1, - flow_direction, density=1.0, - open_boundary_layers=2) + error_str = "The vectors `AB` and `AC` must not be collinear" + + @test_throws ArgumentError(error_str) BoundaryZone(; plane=no_rectangular_plane, + particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + open_boundary_layers=2, + boundary_type=:inflow) + @test_throws ArgumentError(error_str) BoundaryZone(; plane=no_rectangular_plane, + particle_spacing=0.1, + plane_normal=-flow_direction, + density=1.0, + open_boundary_layers=2, + boundary_type=:outflow) rectangular_plane = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] flow_direction = [0.0, 1.0, 0.0] - error_str = "`flow_direction` is not normal to inflow plane" + error_str = "`plane_normal` is not normal to the boundary plane" - @test_throws ArgumentError(error_str) InFlow(; plane=rectangular_plane, - particle_spacing=0.1, - flow_direction, density=1.0, - open_boundary_layers=2) + @test_throws ArgumentError(error_str) BoundaryZone(; plane=rectangular_plane, + particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + open_boundary_layers=2, + boundary_type=:inflow) - error_str = "`flow_direction` is not normal to outflow plane" + error_str = "`plane_normal` is not normal to the boundary plane" - @test_throws ArgumentError(error_str) OutFlow(; plane=rectangular_plane, - particle_spacing=0.1, - flow_direction, density=1.0, - open_boundary_layers=2) + @test_throws ArgumentError(error_str) BoundaryZone(; plane=rectangular_plane, + particle_spacing=0.1, + plane_normal=-flow_direction, + density=1.0, + open_boundary_layers=2, + boundary_type=:outflow) end end diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 3760f73a8..bd45182ba 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -34,18 +34,21 @@ @testset "Flow Direction $j" for j in eachindex(flow_directions) flow_direction = flow_directions[j] - inflow = InFlow(; plane=plane_points, particle_spacing, density, - flow_direction, open_boundary_layers) - outflow = OutFlow(; plane=plane_points, particle_spacing, density, - flow_direction, open_boundary_layers) + inflow = BoundaryZone(; plane=plane_points, particle_spacing, density, + plane_normal=flow_direction, open_boundary_layers, + boundary_type=:inflow) + outflow = BoundaryZone(; plane=plane_points, particle_spacing, density, + plane_normal=-flow_direction, open_boundary_layers, + boundary_type=:outflow) boundary_zones = [ inflow, outflow, ] - @testset "`$(nameof(typeof(boundary_zone)))`" for boundary_zone in boundary_zones - sign_ = (boundary_zone isa InFlow) ? 1 : -1 + @testset "`$(TrixiParticles.boundary_type_name(boundary_zone))`" for boundary_zone in boundary_zones + sign_ = (first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow) ? + 1 : -1 fluid = extrude_geometry(plane_points; particle_spacing, n_extrude=4, density, pressure, direction=(sign_ * flow_direction)) @@ -96,13 +99,13 @@ v, u, v0_ode, u0_ode, semi, t1) evaluated_vars1 = boundary_system.cache.characteristics - if boundary_zone isa InFlow + if first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow @test all(isapprox.(evaluated_vars1[1, :], 0.0)) @test all(isapprox.(evaluated_vars1[2, :], 0.0)) @test all(isapprox.(evaluated_vars1[3, 1:n_influenced], J3(t1))) @test all(isapprox.(evaluated_vars1[3, (n_influenced + 1):end], 0.0)) - elseif boundary_zone isa OutFlow + elseif first(typeof(boundary_zone).parameters) === TrixiParticles.OutFlow @test all(isapprox.(evaluated_vars1[1, 1:n_influenced], J1(t1))) @test all(isapprox.(evaluated_vars1[2, 1:n_influenced], J2(t1))) @test all(isapprox.(evaluated_vars1[1:2, (n_influenced + 1):end], 0.0)) @@ -116,13 +119,13 @@ v, u, v0_ode, u0_ode, semi, t2) evaluated_vars2 = boundary_system.cache.characteristics - if boundary_zone isa InFlow + if first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow @test all(isapprox.(evaluated_vars2[1, :], 0.0)) @test all(isapprox.(evaluated_vars2[2, :], 0.0)) @test all(isapprox.(evaluated_vars2[3, 1:n_influenced], J3(t2))) @test all(isapprox.(evaluated_vars2[3, (n_influenced + 1):end], J3(t1))) - elseif boundary_zone isa OutFlow + elseif first(typeof(boundary_zone).parameters) === TrixiParticles.OutFlow @test all(isapprox.(evaluated_vars2[1, 1:n_influenced], J1(t2))) @test all(isapprox.(evaluated_vars2[2, 1:n_influenced], J2(t2))) @test all(isapprox.(evaluated_vars2[1, (n_influenced + 1):end], J1(t1))) diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index 639fa6429..664325a1a 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -1,5 +1,6 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") +include("surface_normal_sph.jl") include("surface_tension.jl") include("viscosity.jl") diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl new file mode 100644 index 000000000..a9469aa25 --- /dev/null +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -0,0 +1,152 @@ +function create_fluid_system(coordinates, velocity, mass, density, particle_spacing; + buffer_size=0, NDIMS=2, smoothing_length=1.0) + smoothing_kernel = SchoenbergCubicSplineKernel{NDIMS}() + sound_speed = 20.0 + fluid_density = 1000.0 + exponent = 1 + clip_negative_pressure = false + density_calculator = SummationDensity() + surface_normal_method = ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length) + reference_particle_spacing = particle_spacing + tspan = (0.0, 0.01) + + fluid = InitialCondition(coordinates=coordinates, velocity=velocity, mass=mass, + density=density, particle_spacing=particle_spacing) + + state_equation = StateEquationCole(sound_speed=sound_speed, + reference_density=fluid_density, + exponent=exponent, + clip_negative_pressure=clip_negative_pressure) + + system = WeaklyCompressibleSPHSystem(fluid, density_calculator, state_equation, + smoothing_kernel, smoothing_length; + surface_normal_method=surface_normal_method, + reference_particle_spacing=reference_particle_spacing, + buffer_size=buffer_size) + + semi = Semidiscretization(system) + ode = semidiscretize(semi, tspan) + + TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0) + + return system, semi, ode +end + +function compute_and_test_surface_normals(system, semi, ode; NDIMS=2) + # Set values within the function if they are not changed + surface_tension = SurfaceTensionAkinci() + + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, system, semi) + u = TrixiParticles.wrap_u(u0_ode, system, semi) + + # Compute the surface normals + TrixiParticles.compute_surface_normal!(system, surface_tension, v, u, v0_ode, u0_ode, + semi, 0.0) + + # After computation, check that surface normals have been computed + @test all(isfinite.(system.cache.surface_normal)) + @test all(isfinite.(system.cache.neighbor_count)) + @test size(system.cache.surface_normal, 1) == NDIMS + + # Use actual neighbor counts instead of random values + nparticles = size(u, 2) + + # Ensure the threshold is reasonable + threshold = 2^ndims(system) + 1 + + # Test the surface normals based on neighbor counts + for i in 1:nparticles + if system.cache.neighbor_count[i] < threshold + @test all(system.cache.surface_normal[:, i] .== 0.0) + else + # For the linear arrangement, surface normals may still be zero + # Adjust the test to account for this possibility + @test true + end + end +end + +@testset "Surface Normal Computation" begin + # Test case 1: Simple linear arrangement of particles + nparticles = 5 + particle_spacing = 1.0 + NDIMS = 2 + + coordinates = zeros(NDIMS, nparticles) + for i in 1:nparticles + coordinates[:, i] = [i * particle_spacing, 0.0] + end + velocity = zeros(NDIMS, nparticles) + mass = fill(1.0, nparticles) + fluid_density = 1000.0 + density = fill(fluid_density, nparticles) + + system, semi, ode = create_fluid_system(coordinates, velocity, mass, density, + particle_spacing; + buffer_size=0, NDIMS=NDIMS) + + compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS) +end + +@testset "Sphere Surface Normals" begin + # Test case 2: Particles arranged in a circle + particle_spacing = 0.25 + radius = 1.0 + center = (0.0, 0.0) + NDIMS = 2 + + # Create a SphereShape (which is a circle in 2D) + sphere_ic = SphereShape(particle_spacing, radius, center, 1000.0) + + coordinates = sphere_ic.coordinates + velocity = zeros(NDIMS, size(coordinates, 2)) + mass = sphere_ic.mass + density = sphere_ic.density + + # To get some what accurate normals we increase the smoothing length unrealistically + system, semi, ode = create_fluid_system(coordinates, velocity, mass, density, + particle_spacing; + buffer_size=0, NDIMS=NDIMS, + smoothing_length=3.0 * particle_spacing) + + compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS) + + nparticles = size(coordinates, 2) + expected_normals = zeros(NDIMS, nparticles) + surface_particles = Int[] + + # Compute expected normals and identify surface particles + for i in 1:nparticles + pos = coordinates[:, i] + r = pos - SVector(center...) + norm_r = norm(r) + + if abs(norm_r - radius) < particle_spacing + expected_normals[:, i] = -r / norm_r + + push!(surface_particles, i) + else + expected_normals[:, i] .= 0.0 + end + end + + # Normalize computed normals + computed_normals = copy(system.cache.surface_normal) + for i in surface_particles + norm_computed = norm(computed_normals[:, i]) + if norm_computed > 0 + computed_normals[:, i] /= norm_computed + end + end + + # Compare computed normals to expected normals for surface particles + for i in surface_particles + @test isapprox(computed_normals[:, i], expected_normals[:, i], atol=0.05) + end + + # Optionally, check that normals for interior particles are zero + # for i in setdiff(1:nparticles, surface_particles) + # @test isapprox(norm(system.cache.surface_normal[:, i]), 0.0, atol=1e-4) + # end +end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 33b06384d..c32854c32 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -126,7 +126,7 @@ system = EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed) - show_compact = "EntropicallyDampedSPHSystem{2}(SummationDensity(), nothing, Val{:smoothing_kernel}(), [0.0, 0.0]) with 2 particles" + show_compact = "EntropicallyDampedSPHSystem{2}(SummationDensity(), nothing, Val{:smoothing_kernel}(), [0.0, 0.0], nothing, nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -138,6 +138,8 @@ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ │ acceleration: …………………………………………… [0.0, 0.0] │ + │ surface tension: …………………………………… nothing │ + │ surface normal method: …………………… nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box end diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index beb54b320..da39ed397 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -6,8 +6,9 @@ # Mock fluid system struct FluidSystemMock2 <: TrixiParticles.FluidSystem{2, Nothing} end - inflow = InFlow(; plane, particle_spacing=0.1, - flow_direction, density=1.0, open_boundary_layers=2) + inflow = BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, density=1.0, + open_boundary_layers=2, boundary_type=:inflow) error_str = "`reference_velocity` must be either a function mapping " * "each particle's coordinates and time to its velocity, " * @@ -52,8 +53,9 @@ reference_pressure=0) end @testset "`show`" begin - inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, - flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + inflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + plane_normal=(1.0, 0.0), density=1.0, + open_boundary_layers=4, boundary_type=:inflow) system = OpenBoundarySPHSystem(inflow; buffer_size=0, boundary_model=BoundaryModelLastiwka(), reference_density=0.0, @@ -71,7 +73,7 @@ │ #buffer_particles: ……………………………… 0 │ │ fluid system: …………………………………………… FluidSystemMock2 │ │ boundary model: ……………………………………… BoundaryModelLastiwka │ - │ boundary: ……………………………………………………… InFlow │ + │ boundary type: ………………………………………… InFlow │ │ flow direction: ……………………………………… [1.0, 0.0] │ │ prescribed velocity: ………………………… constant_vector │ │ prescribed pressure: ………………………… constant_scalar │ @@ -81,8 +83,9 @@ @test repr("text/plain", system) == show_box - outflow = OutFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, - flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + outflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + plane_normal=(1.0, 0.0), density=1.0, open_boundary_layers=4, + boundary_type=:outflow) system = OpenBoundarySPHSystem(outflow; buffer_size=0, boundary_model=BoundaryModelLastiwka(), reference_density=0.0, @@ -100,8 +103,8 @@ │ #buffer_particles: ……………………………… 0 │ │ fluid system: …………………………………………… FluidSystemMock2 │ │ boundary model: ……………………………………… BoundaryModelLastiwka │ - │ boundary: ……………………………………………………… OutFlow │ - │ flow direction: ……………………………………… [1.0, 0.0] │ + │ boundary type: ………………………………………… OutFlow │ + │ flow direction: ……………………………………… [-1.0, -0.0] │ │ prescribed velocity: ………………………… constant_vector │ │ prescribed pressure: ………………………… constant_scalar │ │ prescribed density: …………………………… constant_scalar │ diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 68f69aea8..a9ee98a48 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -193,7 +193,7 @@ smoothing_length, density_diffusion=density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -207,6 +207,7 @@ │ viscosity: …………………………………………………… nothing │ │ density diffusion: ……………………………… Val{:density_diffusion}() │ │ surface tension: …………………………………… nothing │ + │ surface normal method: …………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ source terms: …………………………………………… Nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""