diff --git a/NEWS.md b/NEWS.md index f1940cbc4..d4e2ccdd7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,14 +4,27 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Version 0.4 + +### API Changes + +- Combined transport velocity formulation (TVF) and particle shifting technique (PST) into + one unified framework. + The keyword argument `transport_velocity` now changed to `shifting_technique`. + The `ParticleShiftingCallback` has been removed. To use PST, use the `UpdateCallback` + instead, and pass `shifting_technique=ParticleShiftingTechnique()` to the system. + +- Renamed the keyword argument `tlsph` to `place_on_shell` for `ParticlePackingSystem`, + `sample_boundary`, `extrude_geometry`, `RectangularShape`, and `SphereShape`. + ## Version 0.3.1 ### Features -- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, +- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, which implement a simplified Smagorinsky-type sub-grid-scale viscosity. (#753) -- **Multithreaded Integration Array**: Introduced a new array type for CPU backends +- **Multithreaded Integration Array**: Introduced a new array type for CPU backends that enables multithreaded broadcasting, delivering speed-ups of up to 5× on systems with many threads when combined with thread pinning. (#722) @@ -21,17 +34,17 @@ used in the Julia ecosystem. Notable changes will be documented in this file for - **DXF file format support**: Import complex geometries using the DXF file format. (#821) - **Improved Plane interpolation**: Massively improved interpolation performance for planes (#763). - + ### GPU - Make PST GPU-compatible (#813). - + - Make open boundaries GPU-compatible (#773). - + - Make interpolation GPU-compatible (#812). ### Important Bugfixes - + - Fix validation setups (#801). - Calculate interpolated density instead of computed density when using interpolation (#808). diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 0db8d275e..d035fb6fc 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -45,64 +45,3 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). Modules = [TrixiParticles] Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")] ``` - -## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) -Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow. -To address these problems, [Adami (2013)](@cite Adami2013) modified the advection velocity and added an extra term to the momentum equation. -The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. [Ramachandran (2019)](@cite Ramachandran2019) applied the TVF -also for the [EDAC](@ref edac) scheme. - -The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by - -```math -\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a -``` - -and is obtained at every time-step ``\Delta t`` from - -```math -\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), -``` - -where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field. -The tilde in the second term of the right hand side indicates that the material derivative has an advection part. - -The discretized form of the last term is - -```math - -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, -``` - -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. -Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, -which means that there is a non-vanishing contribution only when particles are disordered. -That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions. -Suggested is a background pressure which is in the order of the reference pressure but can be chosen arbitrarily large when the time-step criterion is adjusted. - -The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is - -```math -\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, -``` - - where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified - advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. - -The discretized form of the momentum equation for a particle ``a`` reads as - -```math -\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. -``` - -Here, ``\tilde{p}_{ab}`` is the density-weighted pressure - -```math -\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, -``` - -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and are given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. - -```@autodocs -Modules = [TrixiParticles] -Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")] -``` diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 51355961f..62a1bc8d0 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -152,8 +152,74 @@ as explained in [Sun2018](@cite Sun2018) on page 29, right above Equation 9. The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation of PST is commonly referred to as ``\delta^+``-SPH. -The Particle Shifting Technique can be applied in form -of the [`ParticleShiftingCallback`](@ref). +To apply particle shifting, use the keyword argument `shifting_technique` in the constructor +of a system that supports it. + + +## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) + +An alternative formulation is the so-called Transport Velocity Formulation (TVF) +by [Adami (2013)](@cite Adami2013). +[Ramachandran (2019)](@cite Ramachandran2019) applied the TVF also for the [EDAC](@ref edac) +scheme. + +The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position +of the particle ``r_a`` from one time step to the next by +```math +\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a +``` +and is obtained at every time step ``\Delta t`` from +```math +\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), +``` +where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` +is a constant background pressure field. +The tilde in the second term of the right-hand side indicates that the material derivative +has an advection part. + +The discretized form of the last term is +```math + -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, +``` +where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. +Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, +the discretization is not 0th-order consistent for **non**-uniform particle distribution, +which means that there is a non-vanishing contribution only when particles are disordered. +That also means that ``p_{\text{background}}`` occurs as pre-factor to correct +the trajectory of a particle resulting in uniform pressure distributions. +Suggested is a background pressure which is in the order of the reference pressure, +but it can be chosen arbitrarily large when the time-step criterion is adjusted. + +The inviscid momentum equation with an additional convection term for a particle +moving with ``\tilde{v}`` is +```math +\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, +``` +where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence +of the modified advection velocity and can be interpreted as the convection of momentum +with the relative velocity ``\tilde{v}-v``. + +The discretized form of the momentum equation for a particle ``a`` reads as +```math +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. +``` +Here, ``\tilde{p}_{ab}`` is the density-weighted pressure +```math +\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, +``` +with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` +and ``b``, respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors +for particle ``a`` and ``b``, respectively, and are given, e.g., for particle ``a``, +as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. + +To apply the TVF, use the keyword argument `shifting_technique` in the constructor +of a system that supports it. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "shifting_techniques.jl")] +``` + ## [Tensile Instability Control](@id tic) diff --git a/examples/dem/collapsing_sand_pile_3d.jl b/examples/dem/collapsing_sand_pile_3d.jl index aa9ba7a56..f463ac241 100644 --- a/examples/dem/collapsing_sand_pile_3d.jl +++ b/examples/dem/collapsing_sand_pile_3d.jl @@ -55,7 +55,8 @@ min_coords_floor = (min_boundary[1] - boundary_thickness, floor_particles = RectangularShape(particle_spacing, (n_particles_floor_x, n_particles_floor_y, n_particles_floor_z), - min_coords_floor; density=boundary_density, tlsph=true) + min_coords_floor; density=boundary_density, + place_on_shell=true) boundary_particles = floor_particles # ========================================================================================== diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index b6582a170..1b8b15cda 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -65,7 +65,7 @@ if wcsph state_equation, smoothing_kernel, pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure)) + shifting_technique=TransportVelocityAdami(pressure)) else state_equation = nothing density_calculator = ContinuityDensity() @@ -73,7 +73,7 @@ else smoothing_length, density_calculator=density_calculator, sound_speed, viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure)) + shifting_technique=TransportVelocityAdami(pressure)) end # ========================================================================================== diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 850c9f2ff..648a98be1 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -60,7 +60,7 @@ smoothing_length = 1.2 * particle_spacing smoothing_kernel = SchoenbergQuarticSplineKernel{2}() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, viscosity=ViscosityAdami(; nu), - transport_velocity=TransportVelocityAdami(pressure), + shifting_technique=TransportVelocityAdami(pressure), acceleration=(acceleration_x, 0.0)) # ========================================================================================== diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 74770c318..9cd4bd95a 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -28,7 +28,7 @@ initial_fluid_size = tank_size initial_velocity = (1.0, 0.0) fluid_density = 1000.0 -sound_speed = initial_velocity[1] +sound_speed = 10 * initial_velocity[1] state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7) @@ -48,6 +48,7 @@ viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, + shifting_technique=nothing, pressure_acceleration=nothing) # ========================================================================================== diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 3544ac7a8..4449a5738 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -12,7 +12,7 @@ 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 = 4 @@ -21,7 +21,7 @@ boundary_layers = 4 # 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 = 8 +open_boundary_layers = 6 # ========================================================================================== # ==== Experiment Setup @@ -30,29 +30,31 @@ tspan = (0.0, 2.0) # Boundary geometry and initial fluid particle positions domain_size = (1.0, 0.4) -flow_direction = [1.0, 0.0] reynolds_number = 100 -const prescribed_velocity = 2.0 +const prescribed_velocity = (1.0, 0.0) +flow_direction = [1.0, 0.0] -boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, - domain_size[2]) +open_boundary_size = (particle_spacing * open_boundary_layers, domain_size[2]) 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 - -sound_speed = 20 * prescribed_velocity +sound_speed = 10 * maximum(abs.(prescribed_velocity)) -state_equation = nothing - -pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, - pressure=pressure, n_layers=boundary_layers, +pipe = RectangularTank(particle_spacing, domain_size, domain_size, fluid_density, + n_layers=boundary_layers, velocity=prescribed_velocity, faces=(false, false, true, true)) -# Shift pipe walls in negative x-direction for the inflow -pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers +min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0) +inlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_inlet, + faces=(false, false, true, true)) + +min_coords_outlet = (pipe.fluid_size[1], 0.0) +outlet = RectangularTank(particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_outlet, + faces=(false, false, true, true)) NDIMS = ndims(pipe.fluid) @@ -60,32 +62,44 @@ n_buffer_particles = 5 * pipe.n_particles_per_dimension[2]^(NDIMS - 1) # ========================================================================================== # ==== Fluid -wcsph = false +wcsph = true smoothing_length = 1.5 * particle_spacing smoothing_kernel = WendlandC2Kernel{NDIMS}() fluid_density_calculator = ContinuityDensity() -kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number +kinematic_viscosity = maximum(prescribed_velocity) * domain_size[2] / reynolds_number viscosity = ViscosityAdami(nu=kinematic_viscosity) fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, sound_speed, viscosity=viscosity, density_calculator=fluid_density_calculator, + shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) # Alternatively the WCSPH scheme can be used if wcsph 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) + exponent=1) + density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, state_equation, smoothing_kernel, + density_diffusion=density_diffusion, smoothing_length, viscosity=viscosity, + shifting_technique=ParticleShiftingTechnique(), + buffer_size=n_buffer_particles) +else + # Alternatively the EDAC scheme can be used + state_equation = nothing + + fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, + smoothing_length, sound_speed, + viscosity=viscosity, + density_calculator=fluid_density_calculator, + shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) end @@ -96,75 +110,71 @@ 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) + return SVector(prescribed_velocity) end -open_boundary_model = BoundaryModelLastiwka() +open_boundary_model = BoundaryModelTafuniMirroring(; mirror_method=ZerothOrderMirroring()) +reference_velocity_in = velocity_function2d +reference_pressure_in = nothing +reference_density_in = nothing boundary_type_in = InFlow() plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers, density=fluid_density, particle_spacing, - boundary_type=boundary_type_in) - -reference_velocity_in = velocity_function2d -reference_pressure_in = pressure -reference_density_in = fluid_density -open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles, - reference_density=reference_density_in, - reference_pressure=reference_pressure_in, - reference_velocity=reference_velocity_in) - + reference_density=reference_density_in, + reference_pressure=reference_pressure_in, + reference_velocity=reference_velocity_in, + initial_condition=inlet.fluid, boundary_type=boundary_type_in) + +reference_velocity_out = nothing +reference_pressure_out = nothing +reference_density_out = nothing boundary_type_out = OutFlow() -plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]]) +plane_out = ([pipe.fluid_size[1], 0.0], [pipe.fluid_size[1], domain_size[2]]) outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction), open_boundary_layers, density=fluid_density, particle_spacing, - boundary_type=boundary_type_out) - -reference_velocity_out = velocity_function2d -reference_pressure_out = pressure -reference_density_out = fluid_density -open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles, - reference_density=reference_density_out, - reference_pressure=reference_pressure_out, - reference_velocity=reference_velocity_out) + reference_density=reference_density_out, + reference_pressure=reference_pressure_out, + reference_velocity=reference_velocity_out, + initial_condition=outlet.fluid, boundary_type=boundary_type_out) + +open_boundary = OpenBoundarySPHSystem(inflow, outflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles) + # ========================================================================================== # ==== Boundary -viscosity_boundary = ViscosityAdami(nu=1e-4) -boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, +wall = union(pipe.boundary, inlet.boundary, outlet.boundary) +viscosity_boundary = viscosity +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, AdamiPressureExtrapolation(), state_equation=state_equation, viscosity=viscosity_boundary, smoothing_kernel, smoothing_length) -boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) +boundary_system = BoundarySPHSystem(wall, boundary_model) # ========================================================================================== # ==== Simulation -min_corner = minimum(pipe.boundary.coordinates .- particle_spacing, dims=2) -max_corner = maximum(pipe.boundary.coordinates .+ particle_spacing, dims=2) +min_corner = minimum(wall.coordinates .- particle_spacing, dims=2) +max_corner = maximum(wall.coordinates .+ particle_spacing, dims=2) nhs = GridNeighborhoodSearch{NDIMS}(; cell_list=FullGridCellList(; min_corner, max_corner), update_strategy=ParallelUpdate()) -semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, - boundary_system, neighborhood_search=nhs, +semi = Semidiscretization(fluid_system, open_boundary, boundary_system, + neighborhood_search=nhs, parallelization_backend=PolyesterBackend()) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") -particle_shifting = ParticleShiftingCallback() extra_callback = nothing -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), - particle_shifting, extra_callback) +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback) sol = solve(ode, RDPK3SpFSAL35(), abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl index 10035f430..595859aff 100644 --- a/examples/fluid/pipe_flow_3d.jl +++ b/examples/fluid/pipe_flow_3d.jl @@ -22,28 +22,24 @@ open_boundary_layers = 6 # ========================================================================================== # ==== Experiment Setup -tspan = (0.0, 2.0) - -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 +tspan = (0.0, 0.5) 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]) - +const prescribed_velocity = (1.0, 0.0, 0.0) flow_direction = [1.0, 0.0, 0.0] +open_boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) +min_coords_inlet = (-open_boundary_layers * particle_spacing, 0.0, 0.0) +min_coords_outlet = (-open_boundary_layers * particle_spacing, 0.0, 0.0) + # setup simulation trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - domain_size=domain_size, boundary_size=boundary_size, + domain_size=domain_size, open_boundary_size=open_boundary_size, flow_direction=flow_direction, faces=(false, false, true, true, true, true), - tspan=tspan, reference_velocity=velocity_function3d, - open_boundary_layers=open_boundary_layers, + tspan=tspan, prescribed_velocity=prescribed_velocity, + open_boundary_layers=open_boundary_layers, min_coords_inlet=min_coords_inlet, + min_coords_outlet=min_coords_outlet, plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0], [0.0, 0.0, domain_size[3]]), plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0], diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index e6e74fbc4..8a99c455c 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -86,13 +86,13 @@ if wcsph pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, viscosity=ViscosityAdami(; nu), - transport_velocity=TransportVelocityAdami(background_pressure)) + shifting_technique=TransportVelocityAdami(background_pressure)) else density_calculator = SummationDensity() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, density_calculator=density_calculator, - transport_velocity=TransportVelocityAdami(background_pressure), + shifting_technique=TransportVelocityAdami(background_pressure), viscosity=ViscosityAdami(; nu)) end diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 529fd5652..44adf7bc5 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -83,18 +83,18 @@ solid_particle_spacing = thickness / (n_particles_x - 1) n_particles_y = round(Int, length_beam / solid_particle_spacing) + 1 # The bottom layer is sampled separately below. Note that the `RectangularShape` puts the -# first particle half a particle spacing away from the boundary, which is correct for fluids, -# but not for solids. We therefore need to pass `tlsph=true`. +# first particle half a particle spacing away from the shell of the shape, which is +# correct for fluids, but not for solids. We therefore need to pass `place_on_shell=true`. # # The right end of the plate is 0.2 from the right end of the tank. plate_position = 0.6 - n_particles_x * solid_particle_spacing plate = RectangularShape(solid_particle_spacing, (n_particles_x, n_particles_y - 1), (plate_position, solid_particle_spacing), - density=solid_density, tlsph=true) + density=solid_density, place_on_shell=true) fixed_particles = RectangularShape(solid_particle_spacing, (n_particles_x, 1), (plate_position, 0.0), - density=solid_density, tlsph=true) + density=solid_density, place_on_shell=true) solid = union(plate, fixed_particles) diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl index f5b42ecaa..d4f8b206f 100644 --- a/examples/fsi/dam_break_plate_2d.jl +++ b/examples/fsi/dam_break_plate_2d.jl @@ -57,15 +57,15 @@ solid_particle_spacing = thickness / (n_particles_x - 1) n_particles_y = round(Int, length_beam / solid_particle_spacing) + 1 # The bottom layer is sampled separately below. Note that the `RectangularShape` puts the -# first particle half a particle spacing away from the boundary, which is correct for fluids, -# but not for solids. We therefore need to pass `tlsph=true`. +# first particle half a particle spacing away from the shell of the shape, which is +# correct for fluids, but not for solids. We therefore need to pass `place_on_shell=true`. plate = RectangularShape(solid_particle_spacing, (n_particles_x, n_particles_y - 1), (2initial_fluid_size[1], solid_particle_spacing), - density=solid_density, tlsph=true) + density=solid_density, place_on_shell=true) fixed_particles = RectangularShape(solid_particle_spacing, (n_particles_x, 1), (2initial_fluid_size[1], 0.0), - density=solid_density, tlsph=true) + density=solid_density, place_on_shell=true) solid = union(plate, fixed_particles) diff --git a/examples/preprocessing/packing_2d.jl b/examples/preprocessing/packing_2d.jl index fc8b33026..0975b6e2a 100644 --- a/examples/preprocessing/packing_2d.jl +++ b/examples/preprocessing/packing_2d.jl @@ -20,7 +20,7 @@ file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", filename * ". # ========================================================================================== # ==== Packing parameters -tlsph = false +place_on_shell = false # ========================================================================================== # ==== Resolution @@ -50,7 +50,7 @@ shape_sampled = ComplexShape(geometry; particle_spacing, density, # Returns `InitialCondition` boundary_sampled = sample_boundary(signed_distance_field; boundary_density=density, - boundary_thickness, tlsph=tlsph) + boundary_thickness, place_on_shell=place_on_shell) trixi2vtk(shape_sampled) trixi2vtk(boundary_sampled, filename="boundary") @@ -66,12 +66,13 @@ background_pressure = 1.0 smoothing_length = 0.8 * particle_spacing packing_system = ParticlePackingSystem(shape_sampled; smoothing_length=smoothing_length, - signed_distance_field, tlsph=tlsph, + signed_distance_field, place_on_shell=place_on_shell, background_pressure) boundary_system = ParticlePackingSystem(boundary_sampled; smoothing_length=smoothing_length, is_boundary=true, signed_distance_field, - tlsph=tlsph, boundary_compress_factor=0.8, + place_on_shell=place_on_shell, + boundary_compress_factor=0.8, background_pressure) # ========================================================================================== diff --git a/examples/preprocessing/packing_3d.jl b/examples/preprocessing/packing_3d.jl index cb15a255b..ede3433b3 100644 --- a/examples/preprocessing/packing_3d.jl +++ b/examples/preprocessing/packing_3d.jl @@ -24,5 +24,5 @@ boundary_thickness = 8 * particle_spacing trixi_include(joinpath(examples_dir(), "preprocessing", "packing_2d.jl"), density=1000.0, particle_spacing=particle_spacing, file=file, - boundary_thickness=boundary_thickness, tlsph=true, + boundary_thickness=boundary_thickness, place_on_shell=true, save_intervals=false) diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl index 8df52c28f..28d371634 100644 --- a/examples/solid/oscillating_beam_2d.jl +++ b/examples/solid/oscillating_beam_2d.jl @@ -38,7 +38,7 @@ fixed_particles = SphereShape(particle_spacing, clamp_radius + particle_spacing (0.0, elastic_beam.thickness / 2), material.density, cutout_min=(0.0, 0.0), cutout_max=(clamp_radius, elastic_beam.thickness), - tlsph=true) + place_on_shell=true) n_particles_clamp_x = round(Int, clamp_radius / particle_spacing) @@ -48,9 +48,9 @@ n_particles_per_dimension = (round(Int, elastic_beam.length / particle_spacing) # Note that the `RectangularShape` puts the first particle half a particle spacing away # from the boundary, which is correct for fluids, but not for solids. -# We therefore need to pass `tlsph=true`. +# We therefore need to pass `place_on_shell=true`. beam = RectangularShape(particle_spacing, n_particles_per_dimension, - (0.0, 0.0), density=material.density, tlsph=true) + (0.0, 0.0), density=material.density, place_on_shell=true) solid = union(beam, fixed_particles) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 987e20875..023ba1e10 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -65,10 +65,9 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - ParticleShiftingCallback + PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback export ContinuityDensity, SummationDensity -export PenaltyForceGanzenmueller, TransportVelocityAdami +export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel @@ -79,7 +78,8 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr DensityDiffusionAntuono export tensile_instability_control export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni, + PressureMirroring, PressureZeroing, BoundaryModelLastiwkaCharacteristics, + BoundaryModelTafuniMirroring, BernoulliPressureExtrapolation export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring export HertzContactModel, LinearContactModel diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 062a16e3a..e9eb048c3 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -32,4 +32,3 @@ include("post_process.jl") include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") -include("particle_shifting.jl") diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl deleted file mode 100644 index 921370f3b..000000000 --- a/src/callbacks/particle_shifting.jl +++ /dev/null @@ -1,148 +0,0 @@ -@doc raw""" - ParticleShiftingCallback() - -Callback to apply the Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017). -Following the original paper, the callback is applied in every time step and not -in every stage of a multi-stage time integration method to reduce the computational -cost and improve the stability of the scheme. - -See [Callbacks](@ref Callbacks) for more information on how to use this callback. -See [Particle Shifting Technique](@ref shifting) for more information on the method itself. - -!!! warning - The Particle Shifting Technique needs to be disabled close to the free surface - and therefore requires a free surface detection method. This is not yet implemented. - **This callback cannot be used in a free surface simulation.** -""" -function ParticleShiftingCallback() - # The first one is the `condition`, the second the `affect!` - return DiscreteCallback((particle_shifting_condition), particle_shifting!, - save_positions=(false, false)) -end - -# `condition` -function particle_shifting_condition(u, t, integrator) - return true -end - -# `affect!` -function particle_shifting!(integrator) - t = integrator.t - semi = integrator.p - v_ode, u_ode = integrator.u.x - dt = integrator.dt - # Internal cache vector, which is safe to use as temporary array - vu_cache = first(get_tmp_cache(integrator)) - - @trixi_timeit timer() "particle shifting callback" begin - # Update quantities that are stored in the systems. These quantities (e.g. pressure) - # still have the values from the last stage of the previous step if not updated here. - @trixi_timeit timer() "update systems and nhs" begin - # Don't create sub-timers here to avoid cluttering the timer output - @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) - end - - @trixi_timeit timer() "particle shifting" foreach_system(semi) do system - u = wrap_u(u_ode, system, semi) - v = wrap_v(v_ode, system, semi) - particle_shifting!(u, v, system, v_ode, u_ode, semi, vu_cache, dt) - end - end - - # Tell OrdinaryDiffEq that `u` has been modified - u_modified!(integrator, true) - - return integrator -end - -function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt) - return u -end - -function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, - vu_cache, dt) - # Wrap the cache vector to an NDIMS x NPARTICLES matrix. - # We need this buffer because we cannot safely update `u` while iterating over it. - _, u_cache = vu_cache.x - delta_r = wrap_u(u_cache, system, semi) - set_zero!(delta_r) - - # This has similar performance to `maximum(..., eachparticle(system))`, - # but is GPU-compatible. - v_max = maximum(x -> sqrt(dot(x, x)), - reinterpret(reshape, SVector{ndims(system), eltype(v)}, - current_velocity(v, system))) - - # TODO this needs to be adapted to multi-resolution. - # Section 3.2 explains what else needs to be changed. - dx = particle_spacing(system, 1) - Wdx = smoothing_kernel(system, dx, 1) - h = smoothing_length(system, 1) - - foreach_system(semi) do neighbor_system - u_neighbor = wrap_u(u_ode, neighbor_system, semi) - v_neighbor = wrap_v(v_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - rho_a = current_density(v, system, particle) - rho_b = current_density(v_neighbor, neighbor_system, neighbor) - - kernel = smoothing_kernel(system, distance, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - # According to p. 29 below Eq. 9 - R = 2 // 10 - n = 4 - - # Eq. 7 in Sun et al. (2017). - # According to the paper, CFL * Ma can be rewritten as Δt * v_max / h - # (see p. 29, right above Eq. 9), but this does not work when scaling h. - # When setting CFL * Ma = Δt * v_max / (2 * Δx), PST works as expected - # for both small and large smoothing length factors. - # We need to scale - # - quadratically with the smoothing length, - # - linearly with the particle spacing, - # - linearly with the time step. - # See https://github.com/trixi-framework/TrixiParticles.jl/pull/834. - delta_r_ = -dt * v_max * (2 * h)^2 / (2 * dx) * (1 + R * (kernel / Wdx)^n) * - m_b / (rho_a + rho_b) * grad_kernel - - # Write into the buffer - for i in eachindex(delta_r_) - @inbounds delta_r[i, particle] += delta_r_[i] - end - end - end - - # Add δ_r from the buffer to the current coordinates - @threaded semi for particle in eachparticle(system) - for i in axes(delta_r, 1) - @inbounds u[i, particle] += delta_r[i, particle] - end - end - - return u -end - -function Base.show(io::IO, cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - print(io, "ParticleShiftingCallback()") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - - if get(io, :compact, false) - show(io, cb) - else - summary_box(io, "ParticleShiftingCallback") - end -end diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 77496a46c..8eecb8380 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -74,22 +74,32 @@ function (update_callback!::UpdateCallback)(integrator) semi = integrator.p v_ode, u_ode = integrator.u.x - # Update quantities that are stored in the systems. These quantities (e.g. pressure) - # still have the values from the last stage of the previous step if not updated here. - update_systems_and_nhs(v_ode, u_ode, semi, t) - - # Update open boundaries first, since particles might be activated or deactivated - @trixi_timeit timer() "update open boundary" foreach_system(semi) do system - update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) - end - - @trixi_timeit timer() "update particle packing" foreach_system(semi) do system - update_particle_packing(system, v_ode, u_ode, semi, integrator) - end - - # This is only used by the particle packing system and should be removed in the future - @trixi_timeit timer() "update TVF" foreach_system(semi) do system - update_transport_velocity!(system, v_ode, semi) + @trixi_timeit timer() "update callback" begin + # Update quantities that are stored in the systems. These quantities (e.g. pressure) + # still have the values from the last stage of the previous step if not updated here. + @trixi_timeit timer() "update systems and nhs" begin + # Don't create sub-timers here to avoid cluttering the timer output + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) + end + + # Update open boundaries first, since particles might be activated or deactivated + @trixi_timeit timer() "update open boundary" foreach_system(semi) do system + update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) + end + + @trixi_timeit timer() "update particle packing" foreach_system(semi) do system + update_particle_packing(system, v_ode, u_ode, semi, integrator) + end + + # This is only used by the particle packing system and should be removed in the future + @trixi_timeit timer() "update TVF" foreach_system(semi) do system + update_transport_velocity!(system, v_ode, semi) + end + + @trixi_timeit timer() "particle shifting" foreach_system(semi) do system + particle_shifting_from_callback!(u_ode, shifting_technique(system), system, + v_ode, semi, integrator.dt) + end end # Tell OrdinaryDiffEq that `u` has been modified diff --git a/src/general/buffer.jl b/src/general/buffer.jl index 3e1c28e8c..214d9765b 100644 --- a/src/general/buffer.jl +++ b/src/general/buffer.jl @@ -36,6 +36,10 @@ function allocate_buffer(initial_condition, buffer::SystemBuffer) return union(initial_condition, buffer_ic) end +# By default, there is no buffer. +# Dispatch by system type to handle systems that provide a buffer. +@inline buffer(system) = nothing + @inline update_system_buffer!(buffer::Nothing, semi) = buffer # TODO `resize` allocates. Find a non-allocating version diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index 9b0ac40d7..6ca3fa16f 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -191,9 +191,10 @@ function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_syst x_range = range(min_corner[1], max_corner[1], length=n_points_per_dimension[1]) y_range = range(min_corner[2], max_corner[2], length=n_points_per_dimension[2]) - # Generate points within the plane. Use `tlsph=true` to generate points on the boundary + # Generate points within the plane. Use `place_on_shell=true` to generate points + # on the shell of the geometry. point_coords = rectangular_shape_coords(resolution, n_points_per_dimension, min_corner, - tlsph=true) + place_on_shell=true) results = interpolate_points(point_coords, semi, ref_system, v_ode, u_ode, smoothing_length=smoothing_length, diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a9c78b334..66370437a 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -485,7 +485,7 @@ end @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du @inline function add_velocity!(du, v, particle, system::FluidSystem) - # This is zero unless a transport velocity is used + # This is zero unless a shifting technique is used delta_v_ = delta_v(system, particle) for i in 1:ndims(system) @@ -972,17 +972,17 @@ end function check_configuration(system::OpenBoundarySPHSystem, systems, neighborhood_search::PointNeighbors.AbstractNeighborhoodSearch) - (; boundary_model, boundary_zone) = system + (; boundary_model, boundary_zones) = system # Store index of the fluid system. This is necessary for re-linking # in case we use Adapt.jl to create a new semidiscretization. fluid_system_index = findfirst(==(system.fluid_system), systems) system.fluid_system_index[] = fluid_system_index - if boundary_model isa BoundaryModelLastiwka && - boundary_zone isa BoundaryZone{BidirectionalFlow} - throw(ArgumentError("`BoundaryModelLastiwka` needs a specific flow direction. " * - "Please specify inflow and outflow.")) + if boundary_model isa BoundaryModelLastiwkaCharacteristics && + any(zone -> isnothing(zone.flow_direction), boundary_zones) + throw(ArgumentError("`BoundaryModelLastiwkaCharacteristics` needs a specific flow direction. " * + "Please specify `InFlow()` and `OutFlow()`.")) end if first(PointNeighbors.requires_update(neighborhood_search)) @@ -1011,10 +1011,8 @@ function set_system_links(system::OpenBoundarySPHSystem, semi) system.pressure, system.boundary_candidates, system.fluid_candidates, - system.boundary_zone, - system.reference_velocity, - system.reference_pressure, - system.reference_density, + system.boundary_zone_indices, + system.boundary_zones, system.buffer, system.cache) end diff --git a/src/general/system.jl b/src/general/system.jl index 929fb8b42..6c4aed9b0 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -37,17 +37,18 @@ initialize!(system, semi) = system # Number of particles in the system whose positions are to be integrated (corresponds to the size of u and du) @inline n_moving_particles(system) = nparticles(system) -@inline eachparticle(system) = Base.OneTo(nparticles(system)) +@inline eachparticle(system::System) = active_particles(system) +@inline eachparticle(initial_condition) = Base.OneTo(nparticles(initial_condition)) # Wrapper for systems with `SystemBuffer` -@inline each_moving_particle(system) = each_moving_particle(system, system.buffer) +@inline each_moving_particle(system) = each_moving_particle(system, buffer(system)) @inline each_moving_particle(system, ::Nothing) = Base.OneTo(n_moving_particles(system)) -@inline active_coordinates(u, system) = active_coordinates(u, system, system.buffer) +@inline active_coordinates(u, system) = active_coordinates(u, system, buffer(system)) @inline active_coordinates(u, system, ::Nothing) = current_coordinates(u, system) -@inline active_particles(system) = active_particles(system, system.buffer) -@inline active_particles(system, ::Nothing) = eachparticle(system) +@inline active_particles(system) = active_particles(system, buffer(system)) +@inline active_particles(system, ::Nothing) = Base.OneTo(nparticles(system)) # This should not be dispatched by system type. We always expect to get a column of `A`. @propagate_inbounds function extract_svector(A, system, i) @@ -147,9 +148,6 @@ function update_final!(system, v, u, v_ode, u_ode, semi, t) return system end -# Only for systems requiring the use of the `UpdateCallback` -@inline requires_update_callback(system) = false - @inline initial_smoothing_length(system) = smoothing_length(system, nothing) @inline function smoothing_length(system, particle) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index b93e25f32..85311bb2d 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -123,12 +123,12 @@ function trixi2vtk(system_, v_ode_, u_ode_, semi_, t, periodic_box; output_direc write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) # Store particle index - vtk["index"] = active_particles(system) + vtk["index"] = eachparticle(system) vtk["time"] = t vtk["ndims"] = ndims(system) vtk["particle_spacing"] = [particle_spacing(system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] if write_meta_data vtk["solver_version"] = git_hash @@ -262,20 +262,20 @@ end function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true) vtk["velocity"] = view(v, 1:ndims(system), :) vtk["mass"] = [hydrodynamic_mass(system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["radius"] = [particle_radius(system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] return vtk end function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["density"] = [current_density(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] # Indexing the pressure is a workaround for slicing issue (see https://github.com/JuliaSIMD/StrideArrays.jl/issues/88) vtk["pressure"] = [current_pressure(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] if system.surface_normal_method !== nothing vtk["surf_normal"] = [surface_normal(system, particle) @@ -349,9 +349,6 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) else vtk["solver"] = "EDAC" vtk["sound_speed"] = system.sound_speed - vtk["background_pressure_TVF"] = system.transport_velocity isa Nothing ? - "-" : - system.transport_velocity.background_pressure end end @@ -377,7 +374,7 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d n_fixed_particles = nparticles(system) - n_moving_particles(system) vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["jacobian"] = [det(deformation_gradient(system, particle)) for particle in eachparticle(system)] @@ -409,19 +406,11 @@ end function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true) vtk["velocity"] = [current_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["density"] = [current_density(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] vtk["pressure"] = [current_pressure(v, system, particle) - for particle in active_particles(system)] - - if write_meta_data - vtk["boundary_zone"] = type2string(first(typeof(system.boundary_zone).parameters)) - vtk["width"] = round(system.boundary_zone.zone_width, digits=3) - vtk["velocity_function"] = type2string(system.reference_velocity) - vtk["pressure_function"] = type2string(system.reference_pressure) - vtk["density_function"] = type2string(system.reference_density) - end + for particle in eachparticle(system)] return vtk end diff --git a/src/preprocessing/particle_packing/signed_distance.jl b/src/preprocessing/particle_packing/signed_distance.jl index 7e957e124..082deced5 100644 --- a/src/preprocessing/particle_packing/signed_distance.jl +++ b/src/preprocessing/particle_packing/signed_distance.jl @@ -59,7 +59,7 @@ function SignedDistanceField(geometry, particle_spacing; particle_spacing)) grid = rectangular_shape_coords(particle_spacing, n_particles_per_dimension, - min_corner; tlsph=true) + min_corner; place_on_shell=true) points = reinterpret(reshape, SVector{NDIMS, ELTYPE}, grid) end diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index f0772589f..100d0f496 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -6,7 +6,7 @@ smoothing_length_interpolation=smoothing_length, is_boundary=false, boundary_compress_factor=1, neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(), - background_pressure, tlsph=false, fixed_system=false) + background_pressure, place_on_shell=false, fixed_system=false) System to generate body-fitted particles for complex shapes. For more information on the methods, see [particle packing](@ref particle_packing). @@ -18,10 +18,11 @@ For more information on the methods, see [particle packing](@ref particle_packin - `background_pressure`: Constant background pressure to physically pack the particles. A large `background_pressure` can cause high accelerations which requires a properly adjusted time step. -- `tlsph`: With the [`TotalLagrangianSPHSystem`](@ref), particles need to be placed - on the boundary of the shape and not half a particle spacing away, - as for fluids. When `tlsph=true`, particles will be placed - on the boundary of the shape. +- `place_on_shell`: If `place_on_shell=true`, particles will be placed + on the shell of the geometry. For example, + the [`TotalLagrangianSPHSystem`](@ref) requires particles to be placed + on the shell of the geometry and not half a particle spacing away, + as for fluids. - `is_boundary`: When `shape` is inside the geometry that was used to create `signed_distance_field`, set `is_boundary=false`. Otherwise (`shape` is the sampled boundary), set `is_boundary=true`. @@ -64,7 +65,7 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV, smoothing_kernel :: K smoothing_length_interpolation :: ELTYPE background_pressure :: ELTYPE - tlsph :: Bool + place_on_shell :: Bool signed_distance_field :: S is_boundary :: Bool shift_length :: ELTYPE @@ -78,7 +79,8 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV, # See the comments in general/gpu.jl for more details. function ParticlePackingSystem(initial_condition, mass, density, particle_spacing, smoothing_kernel, smoothing_length_interpolation, - background_pressure, tlsph, signed_distance_field, + background_pressure, place_on_shell, + signed_distance_field, is_boundary, shift_length, neighborhood_search, signed_distances, particle_refinement, buffer, fixed_system, cache, advection_velocity) @@ -90,7 +92,7 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV, mass, density, particle_spacing, smoothing_kernel, smoothing_length_interpolation, - background_pressure, tlsph, + background_pressure, place_on_shell, signed_distance_field, is_boundary, shift_length, neighborhood_search, signed_distances, particle_refinement, @@ -105,7 +107,8 @@ function ParticlePackingSystem(shape::InitialCondition; smoothing_length_interpolation=smoothing_length, is_boundary=false, boundary_compress_factor=1, neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(), - background_pressure, tlsph=false, fixed_system=false) + background_pressure, place_on_shell=false, + fixed_system=false) NDIMS = ndims(shape) ELTYPE = eltype(shape) mass = copy(shape.mass) @@ -144,12 +147,12 @@ function ParticlePackingSystem(shape::InitialCondition; # Its value is negative if the particle is inside the geometry. # Otherwise (if outside), the value is positive. if is_boundary - offset = tlsph ? shape.particle_spacing : shape.particle_spacing / 2 + offset = place_on_shell ? shape.particle_spacing : shape.particle_spacing / 2 shift_length = -boundary_compress_factor * signed_distance_field.max_signed_distance - offset else - shift_length = tlsph ? zero(ELTYPE) : shape.particle_spacing / 2 + shift_length = place_on_shell ? zero(ELTYPE) : shape.particle_spacing / 2 end cache = (; create_cache_refinement(shape, particle_refinement, smoothing_length)...) @@ -158,7 +161,7 @@ function ParticlePackingSystem(shape::InitialCondition; return ParticlePackingSystem(shape, mass, density, shape.particle_spacing, smoothing_kernel, smoothing_length_interpolation, - background_pressure, tlsph, signed_distance_field, + background_pressure, place_on_shell, signed_distance_field, is_boundary, shift_length, nhs, fill(zero(ELTYPE), nparticles(shape)), particle_refinement, nothing, fixed_system, cache, advection_velocity) @@ -183,7 +186,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::ParticlePackingSystem) system.neighborhood_search |> typeof |> nameof) summary_line(io, "#particles", nparticles(system)) summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "tlsph", system.tlsph ? "yes" : "no") + summary_line(io, "place_on_shell", system.place_on_shell ? "yes" : "no") summary_line(io, "boundary", system.is_boundary ? "yes" : "no") summary_footer(io) end @@ -215,7 +218,7 @@ end function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true) vtk["velocity"] = [advection_velocity(v, system, particle) - for particle in active_particles(system)] + for particle in eachparticle(system)] if write_meta_data vtk["signed_distances"] = system.signed_distances end @@ -330,8 +333,8 @@ function constrain_particle!(u, system, particle, distance_signed, normal_vector (; shift_length) = system # For fluid particles: - # - `tlsph = true`: `shift_length = 0` - # - `tlsph = false`: `shift_length = particle_spacing / 2` + # - `place_on_shell = true`: `shift_length = 0` + # - `place_on_shell = false`: `shift_length = particle_spacing / 2` # For boundary particles: # `shift_length` is the thickness of the boundary. if distance_signed >= -shift_length @@ -346,7 +349,7 @@ function constrain_particle!(u, system, particle, distance_signed, normal_vector system.is_boundary || return u particle_spacing = system.initial_condition.particle_spacing - shift_length_inner = system.tlsph ? particle_spacing : particle_spacing / 2 + shift_length_inner = system.place_on_shell ? particle_spacing : particle_spacing / 2 if distance_signed < shift_length_inner shift = (distance_signed - shift_length_inner) * normal_vector diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 30e589e1c..6a314cf34 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -90,21 +90,26 @@ bidirectional_flow = BoundaryZone(; plane=plane_points, plane_normal, particle_s !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct BoundaryZone{BT, IC, S, ZO, ZW, FD, PN} +struct BoundaryZone{IC, S, ZO, ZW, FD, PN, R} initial_condition :: IC spanning_set :: S zone_origin :: ZO zone_width :: ZW flow_direction :: FD plane_normal :: PN - boundary_type :: BT + reference_values :: R average_inflow_velocity :: Bool + prescribed_density :: Bool + prescribed_pressure :: Bool + prescribed_velocity :: Bool end function BoundaryZone(; plane, plane_normal, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, - open_boundary_layers::Integer, boundary_type=BidirectionalFlow(), - average_inflow_velocity=true) + open_boundary_layers::Integer, average_inflow_velocity=true, + boundary_type=BidirectionalFlow(), + reference_density=nothing, reference_pressure=nothing, + reference_velocity=nothing) if open_boundary_layers <= 0 throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) end @@ -130,9 +135,118 @@ function BoundaryZone(; plane, plane_normal, density, particle_spacing, extrude_geometry, open_boundary_layers; boundary_type=boundary_type) + NDIMS = ndims(ic) + ELTYPE = eltype(ic) + if !(reference_velocity isa Function || isnothing(reference_velocity) || + (reference_velocity isa Vector && length(reference_velocity) == NDIMS)) + throw(ArgumentError("`reference_velocity` must be either a function mapping " * + "each particle's coordinates and time to its velocity, " * + "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("`velocity` function must be of dimension $NDIMS")) + end + end + # We need this dummy for type stability reasons + velocity_dummy = SVector(ntuple(dim -> ELTYPE(Inf), NDIMS)) + velocity_ref = wrap_reference_function(reference_velocity, velocity_dummy) + end + + if !(reference_pressure isa Function || reference_pressure isa Real || + isnothing(reference_pressure)) + throw(ArgumentError("`reference_pressure` must be either a function mapping " * + "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 + # We need this dummy for type stability reasons + pressure_dummy = ELTYPE(Inf) + pressure_ref = wrap_reference_function(reference_pressure, pressure_dummy) + end + + if !(reference_density isa Function || reference_density isa Real || + isnothing(reference_density)) + throw(ArgumentError("`reference_density` must be either a function mapping " * + "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 + # We need this dummy for type stability reasons + density_dummy = ELTYPE(Inf) + density_ref = wrap_reference_function(reference_density, density_dummy) + end + + prescribed_pressure = isnothing(reference_pressure) ? false : true + prescribed_density = isnothing(reference_density) ? false : true + prescribed_velocity = isnothing(reference_velocity) ? false : true + + reference_values = (reference_velocity=velocity_ref, reference_pressure=pressure_ref, + reference_density=density_ref) + + coordinates_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, ic.coordinates) + + if prescribed_pressure + ic.pressure .= pressure_ref.(coordinates_svector, 0) + end + if prescribed_density + ic.density .= density_ref.(coordinates_svector, 0) + ic.mass .= ic.density * ic.particle_spacing^NDIMS + end + if prescribed_velocity + ic.velocity .= stack(velocity_ref.(coordinates_svector, 0)) + end + return BoundaryZone(ic, spanning_set_, zone_origin, zone_width, - flow_direction, plane_normal_, boundary_type, - average_inflow_velocity) + flow_direction, plane_normal_, reference_values, + average_inflow_velocity, prescribed_density, prescribed_pressure, + prescribed_velocity) +end + +function boundary_type_name(boundary_zone::BoundaryZone) + (; flow_direction, plane_normal) = boundary_zone + + # Bidirectional flow + isnothing(flow_direction) && return "bidirectional_flow" + + # Outflow + signbit(dot(flow_direction, plane_normal)) && return "outflow" + + # Inflow + return "inflow" +end + +function Base.show(io::IO, boundary_zone::BoundaryZone) + @nospecialize boundary_zone # reduce precompilation time + + print(io, "BoundaryZone(") + print(io, ") with ", nparticles(boundary_zone.initial_condition), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone) + @nospecialize boundary_zone # reduce precompilation time + + if get(io, :compact, false) + show(io, boundary_zone) + else + summary_header(io, "BoundaryZone") + summary_line(io, "boundary type", boundary_type_name(boundary_zone)) + summary_line(io, "#particles", nparticles(boundary_zone.initial_condition)) + summary_line(io, "width", round(boundary_zone.zone_width, digits=3)) + summary_footer(io) + end end function set_up_boundary_zone(plane, plane_normal, flow_direction, density, @@ -178,7 +292,7 @@ function set_up_boundary_zone(plane, plane_normal, flow_direction, density, throw(ArgumentError("`plane_normal` is not normal to the boundary plane")) end - if boundary_type isa InFlow + if boundary_type == InFlow() # First vector of `spanning_vectors` is normal to the boundary plane dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction) @@ -186,7 +300,7 @@ function set_up_boundary_zone(plane, plane_normal, flow_direction, density, # Flip the normal vector to point in the opposite direction of `flow_direction`. spanning_set[:, 1] .*= -sign(dot_flow) - elseif boundary_type isa OutFlow + elseif boundary_type == OutFlow() # First vector of `spanning_vectors` is normal to the boundary plane dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction) @@ -194,7 +308,7 @@ function set_up_boundary_zone(plane, plane_normal, flow_direction, density, # Flip the normal vector to point in `flow_direction`. spanning_set[:, 1] .*= sign(dot_flow) - elseif boundary_type isa BidirectionalFlow + elseif boundary_type == BidirectionalFlow() # Flip the normal vector to point opposite to fluid domain spanning_set[:, 1] .*= -sign(dot_plane_normal) end @@ -237,7 +351,7 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) return hcat(c, edge1, edge2) end -@inline function is_in_boundary_zone(boundary_zone::BoundaryZone, particle_coords) +@inline function is_in_boundary_zone(boundary_zone, particle_coords) (; zone_origin, spanning_set) = boundary_zone particle_position = particle_coords - zone_origin @@ -261,6 +375,27 @@ end return true end +function update_boundary_zone_indices!(system, u, boundary_zones, semi) + set_zero!(system.boundary_zone_indices) + + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + for (zone_id, boundary_zone) in enumerate(boundary_zones) + # Check if boundary particle is in the boundary zone + if is_in_boundary_zone(boundary_zone, particle_coords) + system.boundary_zone_indices[particle] = zone_id + end + end + end + + return system +end + +function current_boundary_zone(system, particle) + return system.boundary_zones[system.boundary_zone_indices[particle]] +end + function remove_outside_particles(initial_condition, spanning_set, zone_origin) (; coordinates, density, particle_spacing) = initial_condition @@ -276,3 +411,46 @@ function remove_outside_particles(initial_condition, spanning_set, zone_origin) return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), particle_spacing) end + +function wrap_reference_function(function_::Nothing, ref_dummy) + # Return a dummy value for type stability + return @inline((coords, t)->ref_dummy) +end + +function wrap_reference_function(function_::Function, ref_dummy) + # Already a function + return function_ +end + +function wrap_reference_function(constant_scalar_::Number, ref_dummy) + return @inline((coords, t)->constant_scalar_) +end + +function wrap_reference_function(constant_vector_::AbstractVector, + ref_dummy::SVector{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} + return @inline((coords, t)->SVector{NDIMS, ELTYPE}(constant_vector_)) +end + +function apply_reference_pressure(system, particle, pos, t) + (; pressure_references) = system.cache + + zone_id = system.boundary_zone_indices[particle] + + return apply_ith_function(pressure_references, zone_id, pos, t) +end + +function apply_reference_density(system, particle, pos, t) + (; density_references) = system.cache + + zone_id = system.boundary_zone_indices[particle] + + return apply_ith_function(density_references, zone_id, pos, t) +end + +function apply_reference_velocity(system, particle, pos, t) + (; velocity_references) = system.cache + + zone_id = system.boundary_zone_indices[particle] + + return apply_ith_function(velocity_references, zone_id, pos, t) +end diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index a37a9cd60..154550767 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -1,5 +1,5 @@ @doc raw""" - BoundaryModelLastiwka(; extrapolate_reference_values=nothing) + BoundaryModelLastiwkaCharacteristics(; extrapolate_reference_values=nothing) Boundary model for [`OpenBoundarySPHSystem`](@ref). This model uses the characteristic variables to propagate the appropriate values @@ -19,58 +19,67 @@ For more information about the method see [description below](@ref method_of_cha Note that even without this extrapolation feature, the reference values don't need to be prescribed - they're computed from the characteristics. """ -struct BoundaryModelLastiwka{T} +struct BoundaryModelLastiwkaCharacteristics{T} extrapolate_reference_values::T - function BoundaryModelLastiwka(; extrapolate_reference_values=nothing) + function BoundaryModelLastiwkaCharacteristics(; extrapolate_reference_values=nothing) return new{typeof(extrapolate_reference_values)}(extrapolate_reference_values) end end # Called from update callback via `update_open_boundary_eachstep!` -@inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka, +@inline function update_boundary_quantities!(system, + boundary_model::BoundaryModelLastiwkaCharacteristics, v, u, v_ode, u_ode, semi, t) - (; density, pressure, cache, boundary_zone, - reference_velocity, reference_pressure, reference_density) = system - (; flow_direction) = boundary_zone - - fluid_system = corresponding_fluid_system(system, semi) + (; density, pressure, cache, boundary_zones, fluid_system) = system sound_speed = system_sound_speed(fluid_system) if !isnothing(boundary_model.extrapolate_reference_values) - (; prescribed_pressure, prescribed_velocity, prescribed_density) = cache v_fluid = wrap_v(v_ode, fluid_system, semi) u_fluid = wrap_u(u_ode, fluid_system, semi) @trixi_timeit timer() "extrapolate and correct values" begin extrapolate_values!(system, boundary_model.extrapolate_reference_values, - v, v_fluid, u, u_fluid, semi, t; - prescribed_pressure, prescribed_velocity, - prescribed_density) + v, v_fluid, u, u_fluid, semi) end end # Update quantities based on the characteristic variables @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_velocity, prescribed_density, prescribed_pressure, + flow_direction) = boundary_zone + particle_position = current_coords(u, system, particle) J1 = cache.characteristics[1, particle] J2 = cache.characteristics[2, particle] J3 = cache.characteristics[3, particle] - rho_ref = reference_value(reference_density, density[particle], - particle_position, t) + if prescribed_density + rho_ref = apply_reference_density(system, particle, particle_position, t) + else + rho_ref = current_density(v, system, particle) + end + density[particle] = rho_ref + ((-J1 + (J2 + J3) / 2) / sound_speed^2) - p_ref = reference_value(reference_pressure, pressure[particle], - particle_position, t) + if prescribed_pressure + p_ref = apply_reference_pressure(system, particle, particle_position, t) + else + p_ref = current_pressure(v, system, particle) + end + pressure[particle] = p_ref + (J2 + J3) / 2 - v_current = current_velocity(v, system, particle) - v_ref = reference_value(reference_velocity, v_current, - particle_position, t) - rho = density[particle] + if prescribed_velocity + v_ref = apply_reference_velocity(system, particle, particle_position, t) + else + v_ref = current_velocity(v, system, particle) + end + + rho = current_density(v, system, particle) v_ = v_ref + ((J2 - J3) / (2 * sound_speed * rho)) * flow_direction for dim in 1:ndims(system) @@ -78,19 +87,21 @@ end end end - if boundary_zone.average_inflow_velocity - # Even if the velocity is prescribed, this boundary model computes the velocity for each particle individually. - # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, - # resulting in a potential numerical instability. Averaging mitigates these effects. - average_velocity!(v, u, system, boundary_model, boundary_zone, semi) + for boundary_zone in boundary_zones + if boundary_zone.average_inflow_velocity + # Even if the velocity is prescribed, this boundary model computes the velocity for each particle individually. + # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, + # resulting in a potential numerical instability. Averaging mitigates these effects. + average_velocity!(v, u, system, boundary_model, boundary_zone, semi) + end end return system end # Called from semidiscretization -function update_boundary_model!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, - semi, t) +function update_boundary_model!(system, ::BoundaryModelLastiwkaCharacteristics, + 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 @@ -103,9 +114,8 @@ end # J2: Propagates downstream to the local flow # J3: Propagates upstream to the local flow function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) - (; volume, cache, boundary_zone) = system + (; volume, cache, fluid_system, density, pressure) = system (; characteristics, previous_characteristics) = cache - fluid_system = corresponding_fluid_system(system, semi) @threaded semi for particle in eachparticle(system) previous_characteristics[1, particle] = characteristics[1, particle] @@ -117,14 +127,68 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) set_zero!(volume) # Evaluate the characteristic variables with the fluid system - evaluate_characteristics!(system, fluid_system, v, u, v_ode, u_ode, semi, t) + v_fluid = wrap_v(v_ode, fluid_system, semi) + u_fluid = wrap_u(u_ode, fluid_system, semi) + + system_coords = current_coordinates(u, system) + fluid_coords = current_coordinates(u_fluid, fluid_system) + sound_speed = system_sound_speed(fluid_system) + + # Loop over all fluid neighbors within the kernel cutoff + foreach_point_neighbor(system, fluid_system, system_coords, fluid_coords, semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_velocity, prescribed_density, prescribed_pressure, + flow_direction) = boundary_zone + + neighbor_position = current_coords(u_fluid, fluid_system, neighbor) + + # Determine current and prescribed quantities + rho_b = current_density(v_fluid, fluid_system, neighbor) + + if prescribed_density + rho_ref = apply_reference_density(system, particle, neighbor_position, t) + else + rho_ref = current_density(v, system, particle) + end + + p_b = current_pressure(v_fluid, fluid_system, neighbor) + + if prescribed_pressure + p_ref = apply_reference_pressure(system, particle, neighbor_position, t) + else + p_ref = current_pressure(v, system, particle) + end + + v_b = current_velocity(v_fluid, fluid_system, neighbor) + + if prescribed_velocity + v_neighbor_ref = apply_reference_velocity(system, particle, neighbor_position, + t) + else + v_neighbor_ref = current_velocity(v, system, particle) + end + + # Determine characteristic variables + density_term = -sound_speed^2 * (rho_b - rho_ref) + pressure_term = p_b - p_ref + velocity_term = rho_b * sound_speed * (dot(v_b - v_neighbor_ref, flow_direction)) + + kernel_ = smoothing_kernel(fluid_system, distance, particle) + + characteristics[1, particle] += (density_term + pressure_term) * kernel_ + characteristics[2, particle] += (velocity_term + pressure_term) * kernel_ + characteristics[3, particle] += (-velocity_term + pressure_term) * kernel_ + + volume[particle] += kernel_ + end # Only some of the in-/outlet particles are in the influence of the fluid particles. # Thus, we compute the characteristics for the particles that are outside the influence # of fluid particles by using the average of the values of the previous time step. # See eq. 27 in Negi (2020) https://doi.org/10.1016/j.cma.2020.113119 @threaded semi for particle in each_moving_particle(system) - # Particle is outside of the influence of fluid particles if isapprox(volume[particle], 0) @@ -161,93 +225,44 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) characteristics[2, particle] /= volume[particle] characteristics[3, particle] /= volume[particle] end - prescribe_conditions!(characteristics, particle, boundary_zone) - end - - return system -end - -function evaluate_characteristics!(system, neighbor_system::FluidSystem, - v, u, v_ode, u_ode, semi, t) - (; 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) - u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) - sound_speed = system_sound_speed(neighbor_system) - - # Loop over all fluid neighbors within the kernel cutoff - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - neighbor_position = current_coords(u_neighbor_system, neighbor_system, neighbor) - - # Determine current and prescribed quantities - rho_b = current_density(v_neighbor_system, neighbor_system, neighbor) - rho_ref = reference_value(reference_density, density[particle], - neighbor_position, t) - p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor) - p_ref = reference_value(reference_pressure, pressure[particle], - neighbor_position, t) + boundary_zone = current_boundary_zone(system, particle) + (; flow_direction, plane_normal) = boundary_zone - v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - v_particle = current_velocity(v, system, particle) - v_neighbor_ref = reference_value(reference_velocity, v_particle, - neighbor_position, t) + # Outflow + if signbit(dot(flow_direction, plane_normal)) + # J3 is prescribed (i.e. determined from the exterior of the domain). + # J1 and J2 is transmitted from the domain interior. + characteristics[3, particle] = zero(eltype(characteristics)) - # Determine characteristic variables - density_term = -sound_speed^2 * (rho_b - rho_ref) - pressure_term = p_b - p_ref - velocity_term = rho_b * sound_speed * (dot(v_b - v_neighbor_ref, flow_direction)) - - kernel_ = smoothing_kernel(neighbor_system, distance, particle) - - characteristics[1, particle] += (density_term + pressure_term) * kernel_ - characteristics[2, particle] += (velocity_term + pressure_term) * kernel_ - characteristics[3, particle] += (-velocity_term + pressure_term) * kernel_ - - volume[particle] += kernel_ + else # Inflow + # Allow only J3 to propagate upstream to the boundary + characteristics[1, particle] = zero(eltype(characteristics)) + characteristics[2, particle] = zero(eltype(characteristics)) + end end return system end -@inline function prescribe_conditions!(characteristics, particle, ::BoundaryZone{OutFlow}) - # J3 is prescribed (i.e. determined from the exterior of the domain). - # J1 and J2 is transmitted from the domain interior. - characteristics[3, particle] = zero(eltype(characteristics)) - - return characteristics -end - -@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)) +# Only apply averaging at the inflow +function average_velocity!(v, u, system, ::BoundaryModelLastiwkaCharacteristics, + boundary_zone, semi) + (; flow_direction, plane_normal) = boundary_zone - return characteristics -end + # Outflow + signbit(dot(flow_direction, plane_normal)) && return v -function average_velocity!(v, u, system, ::BoundaryModelLastiwka, boundary_zone, semi) - # Only apply averaging at the inflow - return v -end - -function average_velocity!(v, u, system, ::BoundaryModelLastiwka, ::BoundaryZone{InFlow}, - semi) + particles_in_zone = findall(particle -> boundary_zone == + current_boundary_zone(system, particle), + each_moving_particle(system)) # Division inside the `sum` closure to maintain GPU compatibility - avg_velocity = sum(each_moving_particle(system)) do particle - return current_velocity(v, system, particle) / system.buffer.active_particle_count[] + avg_velocity = sum(particles_in_zone) do particle + return current_velocity(v, system, particle) / length(particles_in_zone) end - @threaded semi for particle in each_moving_particle(system) + @threaded semi for particle in particles_in_zone # Set the velocity of the ghost node to the average velocity of the fluid domain for dim in eachindex(avg_velocity) @inbounds v[dim, particle] = avg_velocity[dim] diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 45a98ad7c..7ff133aa6 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -47,7 +47,7 @@ The interpolated values at the ghost nodes are then assigned to the correspondin struct ZerothOrderMirroring end @doc raw""" - BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) + BoundaryModelTafuniMirroring(; mirror_method=FirstOrderMirroring()) Boundary model for the `OpenBoundarySPHSystem`. This model implements the method of [Tafuni et al. (2018)](@cite Tafuni2018) to extrapolate the properties from the fluid domain @@ -61,62 +61,55 @@ We provide three different mirroring methods: - [`FirstOrderMirroring`](@ref): Uses a first order correction based on the gradient of the interpolated values . - [`SimpleMirroring`](@ref): Similar to the first order mirroring, but does not use the gradient of the interpolated values. """ -struct BoundaryModelTafuni{MM} +struct BoundaryModelTafuniMirroring{MM} mirror_method::MM end -function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) - return BoundaryModelTafuni(mirror_method) +function BoundaryModelTafuniMirroring(; mirror_method=FirstOrderMirroring()) + return BoundaryModelTafuniMirroring(mirror_method) end -function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni, +function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuniMirroring, v, u, v_ode, u_ode, semi, t) - (; reference_pressure, reference_density, reference_velocity, boundary_zone, - pressure, density, cache) = system - (; prescribed_pressure, prescribed_density, prescribed_velocity) = cache + (; boundary_zones, pressure, density, fluid_system, cache) = system @trixi_timeit timer() "extrapolate and correct values" begin - fluid_system = corresponding_fluid_system(system, semi) - v_fluid = wrap_v(v_ode, fluid_system, semi) u_fluid = wrap_u(u_ode, fluid_system, semi) - extrapolate_values!(system, boundary_model.mirror_method, v, v_fluid, - u, u_fluid, semi, t; prescribed_pressure, - prescribed_density, prescribed_velocity) - end - - if !prescribed_velocity && boundary_zone.average_inflow_velocity - # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain. - # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, - # resulting in a potential numerical instability. Averaging mitigates these effects. - average_velocity!(v, u, system, boundary_zone, semi) + extrapolate_values!(system, boundary_model.mirror_method, + v, v_fluid, u, u_fluid, semi) end - if prescribed_pressure - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) + for boundary_zone in boundary_zones + (; average_inflow_velocity, prescribed_velocity) = boundary_zone - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) + if !prescribed_velocity && average_inflow_velocity + # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain. + # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, + # resulting in a potential numerical instability. Averaging mitigates these effects. + average_velocity!(v, u, system, boundary_zone, semi) end end - if prescribed_density - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) + @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_density, prescribed_pressure, prescribed_velocity) = boundary_zone + + particle_coords = current_coords(u, system, particle) - density[particle] = reference_value(reference_density, density[particle], - particle_coords, t) + if prescribed_pressure + pressure[particle] = apply_reference_pressure(system, particle, + particle_coords, t) end - end - if prescribed_velocity - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) - v_particle = current_velocity(v, system, particle) + if prescribed_density + density[particle] = apply_reference_density(system, particle, + particle_coords, t) + end - v_ref = reference_value(reference_velocity, v_particle, particle_coords, t) + if prescribed_velocity + v_ref = apply_reference_velocity(system, particle, particle_coords, t) for dim in eachindex(v_ref) @inbounds v[dim, particle] = v_ref[dim] @@ -125,16 +118,15 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni end end -update_boundary_model!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system +function update_boundary_model!(system, ::BoundaryModelTafuniMirroring, v, u, v_ode, u_ode, + semi, t) + return system +end function extrapolate_values!(system, mirror_method::Union{FirstOrderMirroring, SimpleMirroring}, - v_open_boundary, v_fluid, u_open_boundary, u_fluid, - semi, t; prescribed_density=false, - prescribed_pressure=false, prescribed_velocity=false) - (; pressure, density, boundary_zone) = system - - fluid_system = corresponding_fluid_system(system, semi) + v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi) + (; pressure, density, fluid_system) = system # Static indices to avoid allocations two_to_end = SVector{ndims(system)}(2:(ndims(system) + 1)) @@ -149,6 +141,9 @@ function extrapolate_values!(system, # We can do this because we require the neighborhood search to support querying neighbors # of arbitrary positions (see `PointNeighbors.requires_update`). @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_density, prescribed_pressure, prescribed_velocity) = boundary_zone + particle_coords = current_coords(u_open_boundary, system, particle) ghost_node_position = mirror_position(particle_coords, boundary_zone) @@ -280,12 +275,8 @@ function extrapolate_values!(system, end function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, - v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi, t; - prescribed_density=false, prescribed_pressure=false, - prescribed_velocity=false) - (; pressure, density, boundary_zone) = system - - fluid_system = corresponding_fluid_system(system, semi) + v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi) + (; pressure, density, fluid_system) = system # Use the fluid-fluid nhs, since the boundary particles are mirrored into the fluid domain nhs = get_neighborhood_search(fluid_system, fluid_system, semi) @@ -297,6 +288,9 @@ function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, # We can do this because we require the neighborhood search to support querying neighbors # of arbitrary positions (see `PointNeighbors.requires_update`). @threaded semi for particle in each_moving_particle(system) + boundary_zone = current_boundary_zone(system, particle) + (; prescribed_pressure, prescribed_density, prescribed_velocity) = boundary_zone + particle_coords = current_coords(u_open_boundary, system, particle) ghost_node_position = mirror_position(particle_coords, boundary_zone) @@ -486,10 +480,15 @@ function mirror_position(particle_coords, boundary_zone) return particle_coords - 2 * dist * boundary_zone.plane_normal end -average_velocity!(v, u, system, boundary_zone, semi) = v +# Only for inflow boundary zones +function average_velocity!(v, u, system, boundary_zone, semi) + (; plane_normal, zone_origin, initial_condition, flow_direction) = boundary_zone -function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, semi) - (; plane_normal, zone_origin, initial_condition) = boundary_zone + # Bidirectional flow + isnothing(flow_direction) && return v + + # Outflow + signbit(dot(flow_direction, plane_normal)) && return v # We only use the extrapolated velocity in the vicinity of the transition region. # Otherwise, if the boundary zone is too large, averaging would be excessively influenced @@ -497,15 +496,20 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se max_dist = initial_condition.particle_spacing * 110 / 100 candidates = findall(x -> dot(x - zone_origin, -plane_normal) <= max_dist, - reinterpret(reshape, SVector{ndims(system), eltype(u)}, - active_coordinates(u, system))) + reinterpret(reshape, SVector{ndims(system), eltype(u)}, u)) + + particles_in_zone = findall(particle -> boundary_zone == + current_boundary_zone(system, particle), + each_moving_particle(system)) + + intersect!(candidates, particles_in_zone) # Division inside the `sum` closure to maintain GPU compatibility avg_velocity = sum(candidates) do particle return current_velocity(v, system, particle) / length(candidates) end - @threaded semi for particle in each_moving_particle(system) + @threaded semi for particle in particles_in_zone # Set the velocity of the ghost node to the average velocity of the fluid domain for dim in eachindex(avg_velocity) @inbounds v[dim, particle] = avg_velocity[dim] @@ -515,10 +519,16 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se return v end -project_velocity_on_plane_normal!(v, system, particle, boundary_zone) = v +# Only for inflow boundary zones +function project_velocity_on_plane_normal!(v, system, particle, boundary_zone) + (; plane_normal, flow_direction) = boundary_zone + + # Bidirectional flow + isnothing(flow_direction) && return v + + # Outflow + signbit(dot(flow_direction, plane_normal)) && return v -function project_velocity_on_plane_normal!(v, system, particle, - boundary_zone::BoundaryZone{InFlow}) # Project `vel` on 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, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index bc21bc70d..db52ef98c 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -29,170 +29,135 @@ Open boundary system for in- and outflow particles. !!! note "Note" The reference values (`reference_velocity`, `reference_pressure`, `reference_density`) can also be set to `nothing`. - In this case, they will either be extrapolated from the fluid domain ([BoundaryModelTafuni](@ref BoundaryModelTafuni)) - or evolved using the characteristic flow variables ([BoundaryModelLastiwka](@ref BoundaryModelLastiwka)). + In this case, they will either be extrapolated from the fluid domain ([BoundaryModelTafuniMirroring](@ref BoundaryModelTafuniMirroring)) + or evolved using the characteristic flow variables ([BoundaryModelLastiwkaCharacteristics](@ref BoundaryModelLastiwkaCharacteristics)). !!! warning "Experimental Implementation" This is an experimental feature and may change in future releases. It is GPU-compatible (e.g., with CUDA.jl and AMDGPU.jl), but currently **not** supported with Metal.jl. """ -struct OpenBoundarySPHSystem{BM, ELTYPE, NDIMS, IC, FS, FSI, ARRAY1D, BC, FC, BZ, RV, - RP, RD, B, C} <: System{NDIMS} - boundary_model :: BM - initial_condition :: IC - fluid_system :: FS - fluid_system_index :: FSI - smoothing_length :: ELTYPE - mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] - density :: ARRAY1D # Array{ELTYPE, 1}: [particle] - volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] - pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] - boundary_candidates :: BC # Array{UInt32, 1}: [particle] - fluid_candidates :: FC # Array{UInt32, 1}: [particle] - boundary_zone :: BZ - reference_velocity :: RV - reference_pressure :: RP - reference_density :: RD - buffer :: B - cache :: C +struct OpenBoundarySPHSystem{BM, ELTYPE, NDIMS, IC, FS, FSI, ARRAY1D, BC, FC, BZI, BZ, + B, C} <: System{NDIMS} + boundary_model :: BM + initial_condition :: IC + fluid_system :: FS + fluid_system_index :: FSI + smoothing_length :: ELTYPE + mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] + density :: ARRAY1D # Array{ELTYPE, 1}: [particle] + volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] + pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] + boundary_candidates :: BC # Array{Bool, 1}: [particle] + fluid_candidates :: FC # Array{Bool, 1}: [particle] + boundary_zone_indices :: BZI # Array{UInt8, 1}: [particle] + boundary_zones :: BZ + buffer :: B + cache :: C end function OpenBoundarySPHSystem(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_length, mass, density, volume, pressure, boundary_candidates, fluid_candidates, - boundary_zone, reference_velocity, - reference_pressure, reference_density, buffer, cache) + boundary_zone_indices, boundary_zone, buffer, cache) OpenBoundarySPHSystem{typeof(boundary_model), eltype(mass), ndims(initial_condition), typeof(initial_condition), typeof(fluid_system), typeof(fluid_system_index), typeof(mass), typeof(boundary_candidates), typeof(fluid_candidates), - typeof(boundary_zone), typeof(reference_velocity), - typeof(reference_pressure), typeof(reference_density), + typeof(boundary_zone_indices), typeof(boundary_zone), typeof(buffer), typeof(cache)}(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_length, mass, density, volume, pressure, boundary_candidates, - fluid_candidates, boundary_zone, - reference_velocity, reference_pressure, - reference_density, buffer, cache) + fluid_candidates, boundary_zone_indices, + boundary_zone, buffer, cache) end -function OpenBoundarySPHSystem(boundary_zone::BoundaryZone; - fluid_system::FluidSystem, - buffer_size::Integer, boundary_model, - reference_velocity=nothing, - reference_pressure=nothing, - reference_density=nothing) - (; initial_condition) = boundary_zone +function OpenBoundarySPHSystem(boundary_zones::Union{BoundaryZone, Nothing}...; + fluid_system::FluidSystem, buffer_size::Integer, + boundary_model) + boundary_zones_ = filter(bz -> !isnothing(bz), boundary_zones) + reference_values_ = map(bz -> bz.reference_values, boundary_zones_) - buffer = SystemBuffer(nparticles(initial_condition), buffer_size) + initial_conditions = union((bz.initial_condition for bz in boundary_zones)...) - initial_condition = allocate_buffer(initial_condition, buffer) + buffer = SystemBuffer(nparticles(initial_conditions), buffer_size) - NDIMS = ndims(initial_condition) + initial_conditions = allocate_buffer(initial_conditions, buffer) - pressure = copy(initial_condition.pressure) - mass = copy(initial_condition.mass) - density = copy(initial_condition.density) - volume = similar(initial_condition.density) + pressure = copy(initial_conditions.pressure) + mass = copy(initial_conditions.mass) + density = copy(initial_conditions.density) + volume = similar(initial_conditions.density) - if !(reference_velocity isa Function || isnothing(reference_velocity) || - (reference_velocity isa Vector && length(reference_velocity) == NDIMS)) - throw(ArgumentError("`reference_velocity` must be either a function mapping " * - "each particle's coordinates and time to its velocity, " * - "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 - - if !(reference_pressure isa Function || reference_pressure isa Real || - isnothing(reference_pressure)) - throw(ArgumentError("`reference_pressure` must be either a function mapping " * - "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 - - if !(reference_density isa Function || reference_density isa Real || - isnothing(reference_density)) - throw(ArgumentError("`reference_density` must be either a function mapping " * - "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 - - cache = create_cache_open_boundary(boundary_model, initial_condition, - reference_density, reference_velocity, - reference_pressure) + cache = create_cache_open_boundary(boundary_model, initial_conditions, + reference_values_) fluid_system_index = Ref(0) smoothing_length = initial_smoothing_length(fluid_system) - boundary_candidates = fill(false, nparticles(initial_condition)) + boundary_candidates = fill(false, nparticles(initial_conditions)) fluid_candidates = fill(false, nparticles(fluid_system)) - return OpenBoundarySPHSystem(boundary_model, initial_condition, fluid_system, + boundary_zone_indices = zeros(UInt8, nparticles(initial_conditions)) + + boundary_zones_new = map(zone -> BoundaryZone(zone.initial_condition, + zone.spanning_set, + zone.zone_origin, + zone.zone_width, + zone.flow_direction, + zone.plane_normal, + nothing, + zone.average_inflow_velocity, + zone.prescribed_density, + zone.prescribed_pressure, + zone.prescribed_velocity), + boundary_zones) + + return OpenBoundarySPHSystem(boundary_model, initial_conditions, fluid_system, fluid_system_index, smoothing_length, mass, density, volume, pressure, boundary_candidates, fluid_candidates, - boundary_zone, reference_velocity_, - reference_pressure_, reference_density_, buffer, cache) + boundary_zone_indices, boundary_zones_new, buffer, cache) end -function create_cache_open_boundary(boundary_model, initial_condition, - reference_density, reference_velocity, - reference_pressure) - ELTYPE = eltype(initial_condition) +function initialize!(system::OpenBoundarySPHSystem, semi) + (; boundary_zones) = system - prescribed_pressure = isnothing(reference_pressure) ? false : true - prescribed_velocity = isnothing(reference_velocity) ? false : true - prescribed_density = isnothing(reference_density) ? false : true + update_boundary_zone_indices!(system, initial_coordinates(system), boundary_zones, semi) - if boundary_model isa BoundaryModelTafuni - return (; prescribed_pressure=prescribed_pressure, - prescribed_density=prescribed_density, - prescribed_velocity=prescribed_velocity) - end + return system +end + +function create_cache_open_boundary(boundary_model, initial_condition, reference_values) + ELTYPE = eltype(initial_condition) + + pressure_references = map(ref -> ref.reference_pressure, reference_values) + density_references = map(ref -> ref.reference_density, reference_values) + velocity_references = map(ref -> ref.reference_velocity, reference_values) - characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) - previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) + if boundary_model isa BoundaryModelLastiwkaCharacteristics + characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) + previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition)) - return (; characteristics=characteristics, - previous_characteristics=previous_characteristics, - prescribed_pressure=prescribed_pressure, - prescribed_density=prescribed_density, prescribed_velocity=prescribed_velocity) + return (; characteristics=characteristics, + previous_characteristics=previous_characteristics, + pressure_references=pressure_references, + density_references=density_references, + velocity_references=velocity_references) + else + return (; pressure_references=pressure_references, + density_references=density_references, + velocity_references=velocity_references) + end 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, boundary_type_name(system.boundary_zone)) print(io, ") with ", nparticles(system), " particles") end @@ -205,13 +170,9 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem) summary_header(io, "OpenBoundarySPHSystem{$(ndims(system))}") summary_line(io, "#particles", nparticles(system)) summary_line(io, "#buffer_particles", system.buffer.buffer_size) + summary_line(io, "#boundary_zones", length(system.boundary_zones)) summary_line(io, "fluid system", type2string(system.fluid_system)) summary_line(io, "boundary model", type2string(system.boundary_model)) - 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)) - summary_line(io, "width", round(system.boundary_zone.zone_width, digits=3)) summary_footer(io) end end @@ -220,13 +181,11 @@ end return ELTYPE end +@inline buffer(system::OpenBoundarySPHSystem) = system.buffer + # The `UpdateCallback` is required to update particle positions between time steps @inline requires_update_callback(system::OpenBoundarySPHSystem) = true -function corresponding_fluid_system(system::OpenBoundarySPHSystem, semi) - return system.fluid_system -end - function smoothing_length(system::OpenBoundarySPHSystem, particle) return system.smoothing_length end @@ -249,19 +208,17 @@ end # This function is called by the `UpdateCallback`, as the integrator array might be modified function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ode, semi, t) + (; boundary_model) = system + u = wrap_u(u_ode, system, semi) v = wrap_v(v_ode, system, semi) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) - # 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 boundary quantities" update_boundary_quantities!(system, - system.boundary_model, - v, u, - v_ode, - u_ode, - semi, t) + # Update density, pressure and velocity based on the specific boundary model + @trixi_timeit timer() "update boundary quantities" begin + update_boundary_quantities!(system, boundary_model, v, u, v_ode, u_ode, semi, t) + end return system end @@ -269,8 +226,7 @@ end update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) = system function check_domain!(system, v, u, v_ode, u_ode, semi) - (; boundary_zone, boundary_candidates, fluid_candidates) = system - fluid_system = corresponding_fluid_system(system, semi) + (; boundary_zones, boundary_candidates, fluid_candidates, fluid_system) = system u_fluid = wrap_u(u_ode, fluid_system, semi) v_fluid = wrap_v(v_ode, fluid_system, semi) @@ -282,6 +238,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle_coords = current_coords(u, system, particle) # Check if boundary particle is outside the boundary zone + boundary_zone = current_boundary_zone(system, particle) if !is_in_boundary_zone(boundary_zone, particle_coords) boundary_candidates[particle] = true end @@ -297,6 +254,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle = crossed_boundary_particles[i] particle_new = available_fluid_particles[i] + boundary_zone = current_boundary_zone(system, particle) convert_particle!(system, fluid_system, boundary_zone, particle, particle_new, v, u, v_fluid, u_fluid) end @@ -310,9 +268,11 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) @threaded semi for fluid_particle in each_moving_particle(fluid_system) fluid_coords = current_coords(u_fluid, fluid_system, fluid_particle) - # Check if fluid particle is in boundary zone - if is_in_boundary_zone(boundary_zone, fluid_coords) - fluid_candidates[fluid_particle] = true + # Check if fluid particle is in any boundary zone + for boundary_zone in boundary_zones + if is_in_boundary_zone(boundary_zone, fluid_coords) + fluid_candidates[fluid_particle] = true + end end end @@ -326,7 +286,7 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) particle = crossed_fluid_particles[i] particle_new = available_boundary_particles[i] - convert_particle!(fluid_system, system, boundary_zone, particle, particle_new, + convert_particle!(fluid_system, system, particle, particle_new, v, u, v_fluid, u_fluid) end @@ -336,39 +296,15 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) # Since particles have been transferred, the neighborhood searches must be updated update_nhs!(semi, u_ode) - return system -end - -# Outflow particle is outside the boundary zone -@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::BoundaryZone{OutFlow}, particle, - particle_new, v, u, v_fluid, u_fluid) - deactivate_particle!(system, particle, u) - - return system -end - -# Inflow particle is outside the boundary zone -@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::BoundaryZone{InFlow}, particle, - particle_new, v, u, v_fluid, u_fluid) - (; spanning_set) = boundary_zone - - # Activate a new particle in simulation domain - transfer_particle!(fluid_system, system, particle, particle_new, v_fluid, u_fluid, v, u) - - # Reset position of boundary particle - for dim in 1:ndims(system) - u[dim, particle] += spanning_set[1][dim] - end + update_boundary_zone_indices!(system, u, boundary_zones, semi) return system end # Buffer particle is outside the boundary zone @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, - boundary_zone::BoundaryZone{BidirectionalFlow}, - particle, particle_new, v, u, v_fluid, u_fluid) + boundary_zone, particle, particle_new, + v, u, v_fluid, u_fluid) relative_position = current_coords(u, system, particle) - boundary_zone.zone_origin # Check if particle is in- or outside the fluid domain. @@ -392,8 +328,7 @@ end # Fluid particle is in boundary zone @inline function convert_particle!(fluid_system::FluidSystem, system, - boundary_zone, particle, particle_new, - v, u, v_fluid, u_fluid) + particle, particle_new, v, u, v_fluid, u_fluid) # Activate particle in boundary zone transfer_particle!(system, fluid_system, particle, particle_new, v, u, v_fluid, u_fluid) @@ -443,32 +378,6 @@ function write_u0!(u0, system::OpenBoundarySPHSystem) return u0 end -wrap_reference_function(::Nothing, ::Val) = nothing - -function wrap_reference_function(function_::Function, ::Val) - # Already a function - return function_ -end - -# Name the function so that the summary box does know which kind of function this is -function wrap_reference_function(constant_scalar_::Number, ::Val) - return constant_scalar(coords, t) = constant_scalar_ -end - -# For vectors and tuples -# Name the function so that the summary box does know which kind of function this is -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, 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, position, t) = quantity - # To account for boundary effects in the viscosity term of the RHS, use the viscosity model # of the neighboring particle systems. @inline function viscosity_model(system::OpenBoundarySPHSystem, diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 08139e2ca..a4d7d280a 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -51,13 +51,12 @@ 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) - # Add convection term (only when using `TransportVelocityAdami`) - dv_tvf = dv_transport_velocity(transport_velocity(particle_system), - particle_system, neighbor_system, - particle, neighbor, - v_particle_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) + # Extra terms in the momentum equation when using a shifting technique + dv_tvf = dv_shifting(shifting_technique(particle_system), + particle_system, neighbor_system, particle, neighbor, + v_particle_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 19e50c441..5bb8feddc 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -3,7 +3,7 @@ smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), - transport_velocity=nothing, + shifting_technique=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, surface_normal_method=nothing, buffer_size=nothing, @@ -30,10 +30,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). - Default is no TVF. +- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity + formulation](@ref transport_velocity_formulation) to use + with this system. Default is no shifting. - `average_pressure_reduction`: Whether to subtract the average pressure of neighboring particles - from the local pressure (default: `true` when using TVF, `false` otherwise). + from the local pressure (default: `true` when using shifting, `false` otherwise). - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -67,7 +68,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, acceleration :: SVector{NDIMS, ELTYPE} correction :: COR pressure_acceleration_formulation :: PF - transport_velocity :: TV + shifting_technique :: TV average_pressure_reduction :: AVGP source_terms :: ST surface_tension :: SRFT @@ -83,8 +84,8 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), - transport_velocity=nothing, - average_pressure_reduction=(!isnothing(transport_velocity)), + shifting_technique=nothing, + average_pressure_reduction=(!isnothing(shifting_technique)), alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), @@ -137,7 +138,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = (; create_cache_density(initial_condition, density_calculator)..., - create_cache_tvf(initial_condition, transport_velocity)..., + create_cache_shifting(initial_condition, shifting_technique)..., create_cache_avg_pressure_reduction(initial_condition, avg_pressure_reduction)..., create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, @@ -161,14 +162,14 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), - typeof(pressure_acceleration), typeof(transport_velocity), + typeof(pressure_acceleration), typeof(shifting_technique), typeof(avg_pressure_reduction), typeof(source_terms), typeof(surface_tension), typeof(surface_normal_method), typeof(buffer), Nothing, typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, acceleration_, correction, - pressure_acceleration, transport_velocity, + pressure_acceleration, shifting_technique, avg_pressure_reduction, source_terms, surface_tension, surface_normal_method, buffer, @@ -218,8 +219,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "tansport velocity formulation", - system.transport_velocity |> typeof |> nameof) + summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "average pressure reduction", typeof(system.average_pressure_reduction).parameters[1] ? "yes" : "no") summary_line(io, "acceleration", system.acceleration) @@ -245,6 +245,8 @@ end return ndims(system) + 2 end +@inline buffer(system::EntropicallyDampedSPHSystem) = system.buffer + system_correction(system::EntropicallyDampedSPHSystem) = system.correction @inline function current_pressure(v, system::EntropicallyDampedSPHSystem, particle) @@ -271,7 +273,7 @@ end @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed -@inline transport_velocity(system::EntropicallyDampedSPHSystem) = system.transport_velocity +@inline shifting_technique(system::EntropicallyDampedSPHSystem) = system.shifting_technique @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) average_pressure(system, system.average_pressure_reduction, particle) @@ -321,7 +323,7 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, 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) update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi) - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) end # No average pressure reduction is used diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index d3e6adf97..59703f9b6 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -193,7 +193,7 @@ end include("pressure_acceleration.jl") include("viscosity.jl") -include("transport_velocity.jl") +include("shifting_techniques.jl") include("surface_tension.jl") include("surface_normal_sph.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl new file mode 100644 index 000000000..98775b3be --- /dev/null +++ b/src/schemes/fluid/shifting_techniques.jl @@ -0,0 +1,291 @@ +abstract type AbstractShiftingTechnique end + +# No shifting for a system by default +@inline shifting_technique(system) = nothing + +# WARNING: Be careful if defining this function for a specific system type. +# The version for a specific system type will override this generic version. +requires_update_callback(system) = requires_update_callback(shifting_technique(system)) +requires_update_callback(::Nothing) = false +requires_update_callback(::AbstractShiftingTechnique) = true + +# This is called from the `UpdateCallback` +particle_shifting_from_callback!(u_ode, shifting, system, v_ode, semi, dt) = u_ode + +create_cache_shifting(initial_condition, ::Nothing) = (;) + +function create_cache_shifting(initial_condition, ::AbstractShiftingTechnique) + delta_v = zeros(eltype(initial_condition), ndims(initial_condition), + nparticles(initial_condition)) + + return (; delta_v) +end + +# `δv` is the correction to the particle velocity due to the shifting. +# Particles are advected with the velocity `v + δv`. +@propagate_inbounds function delta_v(system, particle) + return delta_v(system, shifting_technique(system), particle) +end + +# Zero when no shifting is used +@inline function delta_v(system, shifting, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +@propagate_inbounds function delta_v(system, ::AbstractShiftingTechnique, particle) + return extract_svector(system.cache.delta_v, system, particle) +end + +function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi) + return system +end + +# Additional term in the momentum equation due to the shifting technique +@inline function dv_shifting(shifting, system, neighbor_system, + particle, neighbor, v_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + return zero(grad_kernel) +end + +@doc raw""" + ParticleShiftingTechnique() + +Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017). +Following the original paper, the callback is applied in every time step and not +in every stage of a multi-stage time integration method to reduce the computational +cost and improve the stability of the scheme. + +See [Particle Shifting Technique](@ref shifting) for more information on the method. + +!!! warning + The Particle Shifting Technique needs to be disabled close to the free surface + and therefore requires a free surface detection method. This is not yet implemented. + **This technique cannot be used in a free surface simulation.** +""" +struct ParticleShiftingTechnique <: AbstractShiftingTechnique end + +# Zero because PST is applied in a callback +@inline function delta_v(system, ::ParticleShiftingTechnique, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +function particle_shifting_from_callback!(u_ode, shifting::ParticleShiftingTechnique, + system, v_ode, semi, dt) + @trixi_timeit timer() "particle shifting" begin + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + # Update the shifting velocity + update_shifting_from_callback!(system, shifting, v, u, v_ode, u_ode, semi) + + # Update the particle positions with the shifting velocity + particle_shifting!(u_ode, shifting, system, semi, dt) + end +end + +function update_shifting_from_callback!(system, ::ParticleShiftingTechnique, + v, u, v_ode, u_ode, semi) + (; cache) = system + (; delta_v) = cache + + set_zero!(delta_v) + + # This has similar performance to `maximum(..., eachparticle(system))`, + # but is GPU-compatible. + v_max = maximum(x -> sqrt(dot(x, x)), + reinterpret(reshape, SVector{ndims(system), eltype(v)}, + current_velocity(v, system))) + + # TODO this needs to be adapted to multi-resolution. + # Section 3.2 explains what else needs to be changed. + dx = particle_spacing(system, 1) + Wdx = smoothing_kernel(system, dx, 1) + h = smoothing_length(system, 1) + + foreach_system(semi) do neighbor_system + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_a = current_density(v, system, particle) + rho_b = current_density(v_neighbor, neighbor_system, neighbor) + + kernel = smoothing_kernel(system, distance, particle) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # According to p. 29 below Eq. 9 + R = 2 // 10 + n = 4 + + # Eq. 7 in Sun et al. (2017). + # According to the paper, CFL * Ma can be rewritten as Δt * v_max / h + # (see p. 29, right above Eq. 9), but this does not work when scaling h. + # When setting CFL * Ma = Δt * v_max / (2 * Δx), PST works as expected + # for both small and large smoothing length factors. + # We need to scale + # - quadratically with the smoothing length, + # - linearly with the particle spacing, + # - linearly with the time step. + # See https://github.com/trixi-framework/TrixiParticles.jl/pull/834. + delta_v_ = -v_max * (2 * h)^2 / (2 * dx) * (1 + R * (kernel / Wdx)^n) * + m_b / (rho_a + rho_b) * grad_kernel + + # Write into the buffer + for i in eachindex(delta_v_) + @inbounds delta_v[i, particle] += delta_v_[i] + end + end + end + + return system +end + +function particle_shifting!(u_ode, ::ParticleShiftingTechnique, system, semi, dt) + (; cache) = system + (; delta_v) = cache + + u = wrap_u(u_ode, system, semi) + + # Add δr from the cache to the current coordinates + @threaded semi for particle in eachparticle(system) + for i in axes(delta_v, 1) + @inbounds u[i, particle] += dt * delta_v[i, particle] + end + end + + return u +end + +""" + TransportVelocityAdami(background_pressure::Real) + +Transport Velocity Formulation (TVF) by [Adami et al. (2013)](@cite Adami2013) +to suppress pairing and tensile instability. +See [TVF](@ref transport_velocity_formulation) for more details of the method. + +# Arguments +- `background_pressure`: Background pressure. Suggested is a background pressure which is + on the order of the reference pressure. + +!!! warning + The Transport Velocity Formulation needs to be disabled close to the free surface + and therefore requires a free surface detection method. This is not yet implemented. + **This technique cannot be used in a free surface simulation.** +""" +struct TransportVelocityAdami{T <: Real} <: AbstractShiftingTechnique + background_pressure::T +end + +@inline function dv_shifting(::TransportVelocityAdami, system, neighbor_system, + particle, neighbor, v_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + v_a = current_velocity(v_system, system, particle) + delta_v_a = delta_v(system, particle) + + v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + delta_v_b = delta_v(neighbor_system, neighbor) + + A_a = rho_a * v_a * delta_v_a' + A_b = rho_b * v_b * delta_v_b' + + # The following term depends on the pressure acceleration formulation. + # See the large comment below. In the original paper (Adami et al., 2013), this is + # (V_a^2 + V_b^2) / m_a * ((A_a + A_b) / 2) * ∇W_ab. + # With the most common pressure acceleration formulation, this is + # m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab. + # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the + # `pressure_acceleration` function. + return pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) +end + +function update_shifting!(system, shifting::TransportVelocityAdami, v, u, v_ode, + u_ode, semi) + (; cache, correction) = system + (; delta_v) = cache + (; background_pressure) = shifting + + sound_speed = system_sound_speed(system) + + set_zero!(delta_v) + + foreach_system(semi) do neighbor_system + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + m_a = @inbounds hydrodynamic_mass(system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + + rho_a = @inbounds current_density(v, system, particle) + rho_b = @inbounds current_density(v_neighbor, neighbor_system, neighbor) + + h = smoothing_length(system, particle) + + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # In the original paper (Adami et al., 2013), the transport velocity is applied + # as follows: + # v_{1/2} = v_0 + Δt/2 * a, + # where a is the regular SPH acceleration term (pressure, viscosity, etc.). + # r_1 = r_0 + Δt * (v_{1/2}, + # where ̃v_{1/2} = v_{1/2} + Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ] + # is the transport velocity. + # We call δv_{1/2} = ̃v_{1/2} - v_{1/2} the shifting velocity. + # We will call δv_{1/2} the shifting velocity, which is given by + # δv = -Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ], + # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b. + # This term depends on the pressure acceleration formulation. + # In Zhang et al. (2017), the pressure acceleration term + # m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab + # is used. They consequently changed the shifting velocity to + # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab ]. + # We therefore use the function `pressure_acceleration` to compute the + # shifting velocity according to the used pressure acceleration formulation. + # In most cases, this will be + # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. + # + # In these papers, the shifting velocity is scaled by the time step Δt. + # We generally want the spatial discretization to be independent of the time step. + # Scaling the shifting velocity by the time step would lead to less shifting + # when very small time steps are used for testing/debugging purposes. + # This is especially problematic in TrixiParticles.jl, as the time step can vary + # significantly between different time integration methods (low vs high order). + # In order to eliminate the time step from the shifting velocity, we apply the + # CFL condition used in Adami et al. (2013): + # Δt <= 0.25 * h / c, + # where h is the smoothing length and c is the sound speed. + # Applying this equation as equality yields the shifting velocity + # δv = -p_0 / 8 * h / c * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. + # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the + # `pressure_acceleration` function. + delta_v_ = background_pressure / 8 * h / sound_speed * + pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) + + # Write into the buffer + for i in eachindex(delta_v_) + @inbounds delta_v[i, particle] += delta_v_[i] + end + end + end + + return system +end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl deleted file mode 100644 index fce7375dd..000000000 --- a/src/schemes/fluid/transport_velocity.jl +++ /dev/null @@ -1,158 +0,0 @@ -""" - TransportVelocityAdami(background_pressure::Real) - -Transport Velocity Formulation (TVF) by [Adami et al. (2013)](@cite Adami2013) -to suppress pairing and tensile instability. -See [TVF](@ref transport_velocity_formulation) for more details of the method. - -# Arguments -- `background_pressure`: Background pressure. Suggested is a background pressure which is - on the order of the reference pressure. -""" -struct TransportVelocityAdami{T <: Real} - background_pressure::T -end - -# No TVF for a system by default -@inline transport_velocity(system) = nothing - -create_cache_tvf(initial_condition, ::Nothing) = (;) - -function create_cache_tvf(initial_condition, ::TransportVelocityAdami) - delta_v = zeros(eltype(initial_condition), ndims(initial_condition), - nparticles(initial_condition)) - - return (; delta_v) -end - -# `δv` is the correction to the particle velocity due to the TVF. -# Particles are advected with the velocity `v + δv`. -@propagate_inbounds function delta_v(system, particle) - return delta_v(system, transport_velocity(system), particle) -end - -@propagate_inbounds function delta_v(system, ::TransportVelocityAdami, particle) - return extract_svector(system.cache.delta_v, system, particle) -end - -# Zero when no TVF is used -@inline function delta_v(system, transport_velocity, particle) - return zero(SVector{ndims(system), eltype(system)}) -end - -@inline function dv_transport_velocity(::Nothing, system, neighbor_system, - particle, neighbor, v_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) - return zero(grad_kernel) -end - -@inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, - particle, neighbor, v_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) - v_a = current_velocity(v_system, system, particle) - delta_v_a = delta_v(system, particle) - - v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - delta_v_b = delta_v(neighbor_system, neighbor) - - A_a = rho_a * v_a * delta_v_a' - A_b = rho_b * v_b * delta_v_b' - - # The following term depends on the pressure acceleration formulation. - # See the large comment below. In the original paper (Adami et al., 2013), this is - # (V_a^2 + V_b^2) / m_a * ((A_a + A_b) / 2) * ∇W_ab. - # With the most common pressure acceleration formulation, this is - # m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab. - # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the - # `pressure_acceleration` function. - return pressure_acceleration(system, neighbor_system, particle, neighbor, - m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) -end - -function update_tvf!(system, transport_velocity, v, u, v_ode, u_ode, semi, t) - return system -end - -function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v_ode, - u_ode, semi, t) - (; cache, correction) = system - (; delta_v) = cache - (; background_pressure) = transport_velocity - - sound_speed = system_sound_speed(system) - - set_zero!(delta_v) - - foreach_system(semi) do neighbor_system - v_neighbor = wrap_v(v_ode, neighbor_system, semi) - u_neighbor = wrap_u(u_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - m_a = @inbounds hydrodynamic_mass(system, particle) - m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) - - rho_a = @inbounds current_density(v, system, particle) - rho_b = @inbounds current_density(v_neighbor, neighbor_system, neighbor) - - h = smoothing_length(system, particle) - - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - # In the original paper (Adami et al., 2013), the transport velocity is applied - # as follows: - # v_{1/2} = v_0 + Δt/2 * a, - # where a is the regular SPH acceleration term (pressure, viscosity, etc.). - # r_1 = r_0 + Δt * (v_{1/2}, - # where ̃v_{1/2} = v_{1/2} + Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ] - # is the transport velocity. - # We call δv_{1/2} = ̃v_{1/2} - v_{1/2} the shifting velocity. - # We will call δv_{1/2} the shifting velocity, which is given by - # δv = -Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ], - # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b. - # This term depends on the pressure acceleration formulation. - # In Zhang et al. (2017), the pressure acceleration term - # m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab - # is used. They consequently changed the shifting velocity to - # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab ]. - # We therefore use the function `pressure_acceleration` to compute the - # shifting velocity according to the used pressure acceleration formulation. - # In most cases, this will be - # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. - # - # In these papers, the shifting velocity is scaled by the time step Δt. - # We generally want the spatial discretization to be independent of the time step. - # Scaling the shifting velocity by the time step would lead to less shifting - # when very small time steps are used for testing/debugging purposes. - # This is especially problematic in TrixiParticles.jl, as the time step can vary - # significantly between different time integration methods (low vs high order). - # In order to eliminate the time step from the shifting velocity, we apply the - # CFL condition used in Adami et al. (2013): - # Δt <= 0.25 * h / c, - # where h is the smoothing length and c is the sound speed. - # Applying this equation as equality yields the shifting velocity - # δv = -p_0 / 8 * h / c * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. - # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the - # `pressure_acceleration` function. - delta_v_ = background_pressure / 8 * h / sound_speed * - pressure_acceleration(system, neighbor_system, particle, neighbor, - m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) - - # Write into the buffer - for i in eachindex(delta_v_) - @inbounds delta_v[i, particle] += delta_v_[i] - end - end - end - - return system -end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index af088e0a1..b9cd27f07 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -70,13 +70,12 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Add convection term (only when using `TransportVelocityAdami`) - dv_tvf = dv_transport_velocity(transport_velocity(particle_system), - particle_system, neighbor_system, - particle, neighbor, - v_particle_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) + # Extra terms in the momentum equation when using a shifting technique + dv_tvf = dv_shifting(shifting_technique(particle_system), + particle_system, neighbor_system, particle, neighbor, + v_particle_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index a5f56499b..46a51cdb8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -5,7 +5,7 @@ acceleration=ntuple(_ -> 0.0, NDIMS), viscosity=nothing, density_diffusion=nothing, pressure_acceleration=nothing, - transport_velocity=nothing, + shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, @@ -36,8 +36,9 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. density calculator and the correction method. To use [Tensile Instability Control](@ref tic), pass [`tensile_instability_control`](@ref) here. -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). - Default is no TVF. +- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity + formulation](@ref transport_velocity_formulation) to use + with this system. Default is no shifting. - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -59,7 +60,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. - `color_value`: The value used to initialize the color of particles in the system. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, - PF, TV, ST, B, SRFT, SRFN, PR, C} <: FluidSystem{NDIMS} + PF, SC, ST, B, SRFT, SRFN, PR, C} <: FluidSystem{NDIMS} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -71,7 +72,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF - transport_velocity :: TV + shifting_technique :: SC source_terms :: ST surface_tension :: SRFT surface_normal_method :: SRFN @@ -89,7 +90,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, ndims(smoothing_kernel)), viscosity=nothing, density_diffusion=nothing, pressure_acceleration=nothing, - transport_velocity=nothing, + shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, @@ -147,7 +148,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, smoothing_length)..., - create_cache_tvf(initial_condition, transport_velocity)..., + create_cache_shifting(initial_condition, shifting_technique)..., color=Int(color_value)) # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` @@ -162,7 +163,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, state_equation, smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, - transport_velocity, source_terms, surface_tension, + shifting_technique, source_terms, surface_tension, surface_normal_method, buffer, particle_refinement, cache) end @@ -177,6 +178,7 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.viscosity) print(io, ", ", system.density_diffusion) + print(io, ", ", system.shifting_technique) print(io, ", ", system.surface_tension) print(io, ", ", system.surface_normal_method) if system.surface_normal_method isa ColorfieldSurfaceNormal @@ -207,9 +209,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "state equation", system.state_equation |> typeof |> nameof) summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity) - summary_line(io, "tansport velocity formulation", - system.transport_velocity |> typeof |> nameof) summary_line(io, "density diffusion", system.density_diffusion) + summary_line(io, "shifting technique", system.shifting_technique) 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 @@ -237,6 +238,8 @@ end return ndims(system) + 1 end +@inline buffer(system::WeaklyCompressibleSPHSystem) = system.buffer + system_correction(system::WeaklyCompressibleSPHSystem) = system.correction @inline function current_velocity(v, system::WeaklyCompressibleSPHSystem) @@ -278,7 +281,7 @@ end @inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed -@inline transport_velocity(system::WeaklyCompressibleSPHSystem) = system.transport_velocity +@inline shifting_technique(system::WeaklyCompressibleSPHSystem) = system.shifting_technique function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) @@ -316,7 +319,7 @@ function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, # 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) - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) end function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl index 1b65c7e26..8398d03b2 100644 --- a/src/schemes/solid/discrete_element_method/system.jl +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -37,7 +37,6 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST, CM} <: SolidSystem{NDIM acceleration :: SVector{NDIMS, ELTYPE} source_terms :: ST contact_model :: CM - buffer :: Nothing function DEMSystem(initial_condition, contact_model; damping_coefficient=0.0001, acceleration=ntuple(_ -> 0.0, @@ -65,7 +64,7 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST, CM} <: SolidSystem{NDIM typeof(mass), typeof(source_terms), typeof(contact_model)}(initial_condition, mass, radius, damping_coefficient, acceleration_, source_terms, - contact_model, nothing) + contact_model) end end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index b465ca2d7..69c73c4f6 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -75,7 +75,6 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, boundary_model :: BM penalty_force :: PF source_terms :: ST - buffer :: Nothing end function TotalLagrangianSPHSystem(initial_condition, @@ -119,7 +118,7 @@ function TotalLagrangianSPHSystem(initial_condition, n_moving_particles, young_modulus, poisson_ratio, lame_lambda, lame_mu, smoothing_kernel, smoothing_length, acceleration_, boundary_model, - penalty_force, source_terms, nothing) + penalty_force, source_terms) end function Base.show(io::IO, system::TotalLagrangianSPHSystem) diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index cdb8e146b..032250ce5 100644 --- a/src/setups/complex_shape.jl +++ b/src/setups/complex_shape.jl @@ -79,7 +79,7 @@ end """ sample_boundary(signed_distance_field; - boundary_density, boundary_thickness, tlsph=true) + boundary_density, boundary_thickness, place_on_shell=true) Sample boundary particles of a complex geometry by using the [`SignedDistanceField`](@ref) of the geometry. @@ -90,9 +90,9 @@ of the geometry. # Keywords - `boundary_thickness`: Thickness of the boundary - `boundary_density`: Density of each boundary particle. -- `tlsph` : When `tlsph=true`, boundary particles will be placed +- `place_on_shell`: When `place_on_shell=true`, boundary particles will be placed one particle spacing from the surface of the geometry. - Otherwise when `tlsph=true` (simulating fluid particles), + Otherwise when `place_on_shell=true` (simulating fluid particles), boundary particles will be placed half particle spacing away from the surface. @@ -118,7 +118,7 @@ boundary_sampled = sample_boundary(signed_distance_field; boundary_density=1.0, ``` """ function sample_boundary(signed_distance_field; - boundary_density, boundary_thickness, tlsph=true) + boundary_density, boundary_thickness, place_on_shell=true) (; max_signed_distance, boundary_packing, positions, distances, particle_spacing) = signed_distance_field @@ -158,6 +158,6 @@ function particle_grid(geometry, particle_spacing; end grid = rectangular_shape_coords(particle_spacing, n_particles_per_dimension, - min_corner; tlsph=true) + min_corner; place_on_shell=true) return reinterpret(reshape, SVector{ndims(geometry), eltype(geometry)}, grid) end diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index 498266ecd..6d169762c 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -30,9 +30,11 @@ Returns an [`InitialCondition`](@ref). - `pressure`: Scalar to set the pressure of all particles to this value. This is only used by the [`EntropicallyDampedSPHSystem`](@ref) and will be overwritten when using an initial pressure function in the system. -- `tlsph`: With the [`TotalLagrangianSPHSystem`](@ref), particles need to be placed - on the boundary of the shape and not one particle radius away, as for fluids. - When `tlsph=true`, particles will be placed on the boundary of the shape. +- `place_on_shell`: If `place_on_shell=true`, particles will be placed + on the shell of the geometry. For example, + the [`TotalLagrangianSPHSystem`](@ref) requires particles to be placed + on the shell of the geometry and not half a particle spacing away, + as for fluids. # Examples ```jldoctest; output = false @@ -79,7 +81,7 @@ shape = extrude_geometry(shape; direction, particle_spacing=0.1, n_extrude=4, de This is an experimental feature and may change in any future releases. """ function extrude_geometry(geometry; particle_spacing=-1, direction, n_extrude::Integer, - velocity=zeros(length(direction)), tlsph=false, + velocity=zeros(length(direction)), place_on_shell=false, mass=nothing, density=nothing, pressure=0.0) direction_ = normalize(direction) NDIMS = length(direction_) @@ -95,9 +97,11 @@ function extrude_geometry(geometry; particle_spacing=-1, direction, n_extrude::I throw(ArgumentError("`particle_spacing` must be specified when not extruding an `InitialCondition`")) end - geometry = shift_plane_corners(geometry, direction_, particle_spacing, tlsph) + geometry = shift_plane_corners(geometry, direction_, particle_spacing, place_on_shell) - face_coords, particle_spacing_ = sample_plane(geometry, particle_spacing; tlsph=tlsph) + face_coords, + particle_spacing_ = sample_plane(geometry, particle_spacing; + place_on_shell=place_on_shell) if !isapprox(particle_spacing, particle_spacing_, rtol=5e-2) @info "The desired size is not a multiple of the particle spacing $particle_spacing." * @@ -119,12 +123,13 @@ end # For corners/endpoints of a plane/line, sample the plane/line with particles. # For 2D coordinates or an `InitialCondition`, add a third dimension. -function sample_plane(geometry::AbstractMatrix, particle_spacing; tlsph) +function sample_plane(geometry::AbstractMatrix, particle_spacing; place_on_shell) if size(geometry, 1) == 2 # Extruding a 2D shape results in a 3D shape - # When `tlsph=true`, particles will be placed on the x-y plane - coords = vcat(geometry, fill(tlsph ? 0 : particle_spacing / 2, size(geometry, 2))') + # When `place_on_shell=true`, particles will be placed on the x-y plane + coords = vcat(geometry, + fill(place_on_shell ? 0 : particle_spacing / 2, size(geometry, 2))') # TODO: 2D shapes not only in x-y plane but in any user-defined plane return coords, particle_spacing @@ -133,13 +138,14 @@ function sample_plane(geometry::AbstractMatrix, particle_spacing; tlsph) return geometry, particle_spacing end -function sample_plane(shape::InitialCondition, particle_spacing; tlsph) +function sample_plane(shape::InitialCondition, particle_spacing; place_on_shell) if ndims(shape) == 2 # Extruding a 2D shape results in a 3D shape - # When `tlsph=true`, particles will be placed on the x-y plane + # When `place_on_shell=true`, particles will be placed on the x-y plane coords = vcat(shape.coordinates, - fill(tlsph ? 0 : particle_spacing / 2, size(shape.coordinates, 2))') + fill(place_on_shell ? 0 : particle_spacing / 2, + size(shape.coordinates, 2))') # TODO: 2D shapes not only in x-y plane but in any user-defined plane return coords, particle_spacing @@ -148,13 +154,13 @@ function sample_plane(shape::InitialCondition, particle_spacing; tlsph) return shape.coordinates, particle_spacing end -function sample_plane(plane_points, particle_spacing; tlsph=nothing) +function sample_plane(plane_points, particle_spacing; place_on_shell=nothing) # Convert to tuple - return sample_plane(tuple(plane_points...), particle_spacing; tlsph=nothing) + return sample_plane(tuple(plane_points...), particle_spacing; place_on_shell=nothing) end -function sample_plane(plane_points::NTuple{2}, particle_spacing; tlsph=nothing) +function sample_plane(plane_points::NTuple{2}, particle_spacing; place_on_shell=nothing) # Verify that points are in 2D space if any(length.(plane_points) .!= 2) throw(ArgumentError("all points must be 2D coordinates")) @@ -168,7 +174,7 @@ function sample_plane(plane_points::NTuple{2}, particle_spacing; tlsph=nothing) return coords, particle_spacing_new end -function sample_plane(plane_points::NTuple{3}, particle_spacing; tlsph=nothing) +function sample_plane(plane_points::NTuple{3}, particle_spacing; place_on_shell=nothing) # Verify that points are in 3D space if any(length.(plane_points) .!= 3) throw(ArgumentError("all points must be 3D coordinates")) @@ -209,21 +215,22 @@ function sample_plane(plane_points::NTuple{3}, particle_spacing; tlsph=nothing) return coords, particle_spacing_new end -# Shift corners of the plane/line inwards by half a particle spacing with `tlsph=false` +# Shift corners of the plane/line inwards by half a particle spacing with `place_on_shell=false` # because fluid particles need to be half a particle spacing away from the boundary of the shape. function shift_plane_corners(geometry::Union{AbstractMatrix, InitialCondition}, - direction, particle_spacing, tlsph) + direction, particle_spacing, place_on_shell) return geometry end -function shift_plane_corners(plane_points, direction, particle_spacing, tlsph) - shift_plane_corners(tuple(plane_points...), direction, particle_spacing, tlsph) +function shift_plane_corners(plane_points, direction, particle_spacing, place_on_shell) + shift_plane_corners(tuple(plane_points...), direction, particle_spacing, place_on_shell) end -function shift_plane_corners(plane_points::NTuple{2}, direction, particle_spacing, tlsph) - # With TLSPH, particles need to be AT the min coordinates and not half a particle +function shift_plane_corners(plane_points::NTuple{2}, direction, particle_spacing, + place_on_shell) + # With `place_on_shell`, particles need to be AT the min coordinates and not half a particle # spacing away from it. - (tlsph) && (return plane_points) + (place_on_shell) && (return plane_points) plane_point1 = copy(plane_points[1]) plane_point2 = copy(plane_points[2]) @@ -238,10 +245,11 @@ function shift_plane_corners(plane_points::NTuple{2}, direction, particle_spacin return (plane_point1, plane_point2) end -function shift_plane_corners(plane_points::NTuple{3}, direction, particle_spacing, tlsph) - # With TLSPH, particles need to be AT the min coordinates and not half a particle +function shift_plane_corners(plane_points::NTuple{3}, direction, particle_spacing, + place_on_shell) + # With `place_on_shell`, particles need to be AT the min coordinates and not half a particle # spacing away from it. - (tlsph) && (return plane_points) + (place_on_shell) && (return plane_points) plane_point1 = copy(plane_points[1]) plane_point2 = copy(plane_points[2]) diff --git a/src/setups/rectangular_shape.jl b/src/setups/rectangular_shape.jl index 98160bcac..1c30ef5c1 100644 --- a/src/setups/rectangular_shape.jl +++ b/src/setups/rectangular_shape.jl @@ -3,7 +3,7 @@ velocity=zeros(length(n_particles_per_dimension)), mass=nothing, density=nothing, pressure=0.0, acceleration=nothing, state_equation=nothing, - tlsph=false, loop_order=nothing) + place_on_shell=false, loop_order=nothing) Rectangular shape filled with particles. Returns an [`InitialCondition`](@ref). @@ -40,9 +40,10 @@ Rectangular shape filled with particles. Returns an [`InitialCondition`](@ref). - `state_equation`: When calculating a hydrostatic pressure gradient by setting `acceleration`, the `state_equation` will be used to set the corresponding density. Cannot be used together with `density`. -- `tlsph`: With the [`TotalLagrangianSPHSystem`](@ref), particles need to be placed - on the boundary of the shape and not one particle radius away, as for fluids. - When `tlsph=true`, particles will be placed on the boundary of the shape. +- `place_on_shell`: If `place_on_shell=true`, particles will be placed on the shell of the shape. + For example, the [`TotalLagrangianSPHSystem`](@ref) requires particles + to be placed on the shell of the shape and not half a particle spacing away, + as for fluids. - `coordinates_perturbation`: Add a small random displacement to the particle positions, where the amplitude is `coordinates_perturbation * particle_spacing`. @@ -75,7 +76,7 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord coordinates_perturbation=nothing, mass=nothing, density=nothing, pressure=0.0, acceleration=nothing, state_equation=nothing, - tlsph=false, loop_order=nothing) + place_on_shell=false, loop_order=nothing) if particle_spacing < eps() throw(ArgumentError("`particle_spacing` needs to be positive and larger than $(eps())")) end @@ -95,7 +96,7 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord n_particles = prod(n_particles_per_dimension) coordinates = rectangular_shape_coords(particle_spacing, n_particles_per_dimension, - min_coordinates, tlsph=tlsph, + min_coordinates, place_on_shell=place_on_shell, loop_order=loop_order) if !isnothing(coordinates_perturbation) @@ -190,15 +191,15 @@ function loop_permutation(loop_order, NDIMS::Val{3}) end function rectangular_shape_coords(particle_spacing, n_particles_per_dimension, - min_coordinates; tlsph=false, loop_order=nothing) + min_coordinates; place_on_shell=false, loop_order=nothing) ELTYPE = eltype(particle_spacing) NDIMS = length(n_particles_per_dimension) coordinates = Array{ELTYPE, 2}(undef, NDIMS, prod(n_particles_per_dimension)) - # With TLSPH, particles need to be AT the min coordinates and not half a particle + # With place_on_shell, particles need to be AT the min coordinates and not half a particle # spacing away from it. - if tlsph + if place_on_shell min_coordinates = min_coordinates .- 0.5particle_spacing end diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl index 0233fc899..3c6a7466d 100644 --- a/src/setups/sphere_shape.jl +++ b/src/setups/sphere_shape.jl @@ -1,7 +1,7 @@ """ SphereShape(particle_spacing, radius, center_position, density; sphere_type=VoxelSphere(), n_layers=-1, layer_outwards=false, - cutout_min=(0.0, 0.0), cutout_max=(0.0, 0.0), tlsph=false, + cutout_min=(0.0, 0.0), cutout_max=(0.0, 0.0), place_on_shell=false, velocity=zeros(length(center_position)), mass=nothing, pressure=0.0) Generate a sphere that is either completely filled (by default) @@ -35,18 +35,19 @@ coordinate directions as `cutout_min` and `cutout_max`. cut out of the sphere. - `cutout_max`: Corner in positive coordinate directions of a cuboid that is to be cut out of the sphere. -- `tlsph`: With the [`TotalLagrangianSPHSystem`](@ref), particles need to be placed - on the boundary of the shape and not one particle radius away, as for fluids. - When `tlsph=true`, particles will be placed on the boundary of the shape. -- `velocity`: Either a function mapping each particle's coordinates to its velocity, - or, for a constant fluid velocity, a vector holding this velocity. - Velocity is constant zero by default. -- `mass`: Either `nothing` (default) to automatically compute particle mass from particle - density and spacing, or a function mapping each particle's coordinates to its mass, - or a scalar for a constant mass over all particles. -- `pressure`: Either a function mapping each particle's coordinates to its pressure, - or a scalar for a constant pressure over all particles. This is optional and - only needed when using the [`EntropicallyDampedSPHSystem`](@ref). +- `place_on_shell`: If `place_on_shell=true`, particles will be placed on the shell of the shape. + For example, the [`TotalLagrangianSPHSystem`](@ref) requires particles + to be placed on the shell of the shape and not half a particle spacing away, + as for fluids. +- `velocity`: Either a function mapping each particle's coordinates to its velocity, + or, for a constant fluid velocity, a vector holding this velocity. + Velocity is constant zero by default. +- `mass`: Either `nothing` (default) to automatically compute particle mass from particle + density and spacing, or a function mapping each particle's coordinates to its mass, + or a scalar for a constant mass over all particles. +- `pressure`: Either a function mapping each particle's coordinates to its pressure, + or a scalar for a constant pressure over all particles. This is optional and + only needed when using the [`EntropicallyDampedSPHSystem`](@ref). # Examples ```jldoctest; output = false @@ -89,7 +90,7 @@ SphereShape(0.1, 0.5, (0.2, 0.4, 0.3), 1000.0, sphere_type=RoundSphere()) """ function SphereShape(particle_spacing, radius, center_position, density; sphere_type=VoxelSphere(), n_layers=-1, layer_outwards=false, - cutout_min=(0.0, 0.0), cutout_max=(0.0, 0.0), tlsph=false, + cutout_min=(0.0, 0.0), cutout_max=(0.0, 0.0), place_on_shell=false, velocity=zeros(length(center_position)), mass=nothing, pressure=0) if particle_spacing < eps() throw(ArgumentError("`particle_spacing` needs to be positive and larger than $(eps())")) @@ -99,7 +100,7 @@ function SphereShape(particle_spacing, radius, center_position, density; coordinates = sphere_shape_coords(sphere_type, particle_spacing, radius, SVector{NDIMS}(center_position), - n_layers, layer_outwards, tlsph) + n_layers, layer_outwards, place_on_shell) # Convert tuples to vectors cutout_min_ = collect(cutout_min) @@ -169,13 +170,13 @@ struct RoundSphere{AR} end function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_position, - n_layers, layer_outwards, tlsph) + n_layers, layer_outwards, place_on_shell) if n_layers > 0 if layer_outwards inner_radius = radius outer_radius = radius + n_layers * particle_spacing - if !tlsph + if !place_on_shell # Put first layer of particles half a particle spacing outside of `radius` inner_radius += particle_spacing / 2 outer_radius += particle_spacing / 2 @@ -184,7 +185,7 @@ function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_pos inner_radius = radius - n_layers * particle_spacing outer_radius = radius - if !tlsph + if !place_on_shell # Put first layer of particles half a particle spacing inside of `radius` inner_radius -= particle_spacing / 2 outer_radius -= particle_spacing / 2 @@ -194,7 +195,7 @@ function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_pos outer_radius = radius inner_radius = -1 - if !tlsph + if !place_on_shell # Put first layer of particles half a particle spacing inside of `radius` outer_radius -= particle_spacing / 2 end @@ -225,7 +226,7 @@ function sphere_shape_coords(::VoxelSphere, particle_spacing, radius, center_pos end function sphere_shape_coords(sphere::RoundSphere, particle_spacing, radius, center, - n_layers, layer_outwards, tlsph) + n_layers, layer_outwards, place_on_shell) if n_layers > 0 if layer_outwards inner_radius = radius @@ -233,12 +234,12 @@ function sphere_shape_coords(sphere::RoundSphere, particle_spacing, radius, cent inner_radius = radius - n_layers * particle_spacing end - if !tlsph + if !place_on_shell # Put first layer of particles half a particle spacing outside of inner radius inner_radius += particle_spacing / 2 end else - if tlsph + if place_on_shell # Just create a sphere that is 0.5 particle spacing larger radius += particle_spacing / 2 end diff --git a/src/util.jl b/src/util.jl index 88be2096a..d551761ec 100644 --- a/src/util.jl +++ b/src/util.jl @@ -11,6 +11,16 @@ end @inline foreach_noalloc(func, collection::Tuple{}) = nothing +@inline function apply_ith_function(functions, index, args...) + if index == 1 + # Found the function to apply, apply it and return + return first(functions)(args...) + end + + # Process remaining functions + apply_ith_function(Base.tail(functions), index - 1, args...) +end + # Print informative message at startup function print_startup_message() s = """ diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index ede9e6bbf..97b839c51 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -273,7 +273,19 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - extra_callback=ParticleShiftingCallback()) + shifting_technique=ParticleShiftingTechnique(), + extra_callback=UpdateCallback()) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/periodic_channel_2d.jl with TVF" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "periodic_channel_2d.jl"), + tspan=(0.0, 0.2), + shifting_technique=TransportVelocityAdami(50_000.0), + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @@ -283,68 +295,59 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - extra_callback=ParticleShiftingCallback(), - pressure_acceleration=tensile_instability_control) + shifting_technique=ParticleShiftingTechnique(), + pressure_acceleration=tensile_instability_control, + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwkaCharacteristics (WCSPH)" begin @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + open_boundary_model=BoundaryModelLastiwkaCharacteristics(), joinpath(examples_dir(), "fluid", - "pipe_flow_2d.jl"), - wcsph=true) + "pipe_flow_2d.jl")) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (EDAC)" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwkaCharacteristics (EDAC)" begin + @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), wcsph=false, + open_boundary_model=BoundaryModelLastiwkaCharacteristics(), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl")) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (EDAC)" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuniMirroring (EDAC)" begin + @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), wcsph=false, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - open_boundary_model=BoundaryModelTafuni(), boundary_type_in=BidirectionalFlow(), - boundary_type_out=BidirectionalFlow(), - reference_density_in=nothing, - reference_pressure_in=nothing, - reference_density_out=nothing, - reference_velocity_out=nothing) + boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuniMirroring (WCSPH)" begin @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - wcsph=true, sound_speed=20.0, pressure=0.0, - open_boundary_model=BoundaryModelTafuni(), boundary_type_in=BidirectionalFlow(), - boundary_type_out=BidirectionalFlow(), - reference_density_in=nothing, - reference_pressure_in=nothing, - reference_density_out=nothing, - reference_pressure_out=nothing, - reference_velocity_out=nothing) + boundary_type_out=BidirectionalFlow()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin - steady_state_reached = SteadyStateReachedCallback(; dt=0.002, interval_size=10, + steady_state_reached = SteadyStateReachedCallback(; dt=0.002, interval_size=5, reltol=1e-3) @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelLastiwkaCharacteristics(), extra_callback=steady_state_reached, tspan=(0.0, 1.5), viscosity_boundary=nothing) @@ -354,11 +357,12 @@ end @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`interval`)" begin - steady_state_reached = SteadyStateReachedCallback(; interval=1, interval_size=10, + steady_state_reached = SteadyStateReachedCallback(; interval=1, interval_size=5, reltol=1e-3) @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + open_boundary_model=BoundaryModelLastiwkaCharacteristics(), extra_callback=steady_state_reached, dtmax=2e-3, tspan=(0.0, 1.5), viscosity_boundary=nothing) @@ -368,7 +372,7 @@ end @trixi_testset "fluid/pipe_flow_3d.jl" begin - @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0.0, 0.5), + @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_3d.jl")) @test sol.retcode == ReturnCode.Success diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index 038e0a2b8..4231f86af 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -305,7 +305,7 @@ end end # Test open boundaries and steady-state callback - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwkaCharacteristics (WCSPH)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), @@ -318,7 +318,7 @@ end @test backend == Main.parallelization_backend end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwka (EDAC)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelLastiwkaCharacteristics (EDAC)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), @@ -330,13 +330,13 @@ end @test backend == Main.parallelization_backend end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (EDAC)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuniMirroring (EDAC)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - open_boundary_model=BoundaryModelTafuni(), + open_boundary_model=BoundaryModelTafuniMirroring(), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, @@ -349,7 +349,7 @@ end @test backend == Main.parallelization_backend end - @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuni (WCSPH)" begin + @trixi_testset "fluid/pipe_flow_2d.jl - BoundaryModelTafuniMirroring (WCSPH)" begin @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__, tspan=(0.0f0, 0.5f0), joinpath(examples_dir(), @@ -357,8 +357,8 @@ end "pipe_flow_2d.jl"), wcsph=true, sound_speed=20.0f0, pressure=0.0f0, - open_boundary_model=BoundaryModelTafuni(; - mirror_method=ZerothOrderMirroring()), + open_boundary_model=BoundaryModelTafuniMirroring(; + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, diff --git a/test/general/buffer.jl b/test/general/buffer.jl index 18758590e..50f79ce33 100644 --- a/test/general/buffer.jl +++ b/test/general/buffer.jl @@ -6,15 +6,13 @@ 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()) + reference_density=1.0, reference_pressure=0.0, + reference_velocity=[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) + boundary_model=BoundaryModelLastiwkaCharacteristics(), + buffer_size=0) system_buffer = OpenBoundarySPHSystem(zone; buffer_size=5, - reference_density=0.0, reference_pressure=0.0, - reference_velocity=[0, 0], - boundary_model=BoundaryModelLastiwka(), + boundary_model=BoundaryModelLastiwkaCharacteristics(), fluid_system=FluidSystemMock3()) n_particles = nparticles(system) diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index dea8c356a..862f2ba37 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -140,7 +140,7 @@ v_ode = vcat(vec(v1), v2) # Avoid `SystemBuffer` barrier - TrixiParticles.each_moving_particle(system::Union{System1, System2}) = TrixiParticles.eachparticle(system) + TrixiParticles.each_moving_particle(system::Union{System1, System2}) = Base.OneTo(nparticles(system)) TrixiParticles.add_source_terms!(dv_ode, v_ode, u_ode, semi, 0.0) diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index b14906b63..7efbd9cc1 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -1,4 +1,95 @@ @testset verbose=true "Boundary Zone" begin + @testset "`show`" begin + inflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + plane_normal=(1.0, 0.0), density=1.0, + reference_density=0.0, + reference_pressure=0.0, + reference_velocity=[0.0, 0.0], + open_boundary_layers=4, boundary_type=InFlow()) + + show_compact = "BoundaryZone() with 80 particles" + @test repr(inflow) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ BoundaryZone │ + │ ════════════ │ + │ boundary type: ………………………………………… inflow │ + │ #particles: ………………………………………………… 80 │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", inflow) == show_box + + outflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + reference_density=0.0, + reference_pressure=0.0, + reference_velocity=[0.0, 0.0], + plane_normal=(1.0, 0.0), density=1.0, open_boundary_layers=4, + boundary_type=OutFlow()) + + show_compact = "BoundaryZone() with 80 particles" + @test repr(outflow) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ BoundaryZone │ + │ ════════════ │ + │ boundary type: ………………………………………… outflow │ + │ #particles: ………………………………………………… 80 │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", outflow) == show_box + end + + @testset verbose=true "Illegal Inputs" begin + plane = ([0.0, 0.0], [0.0, 1.0]) + flow_direction = (1.0, 0.0) + + error_str = "`reference_velocity` must be either a function mapping " * + "each particle's coordinates and time to its velocity, " * + "an array where the ``i``-th column holds the velocity of particle ``i`` " * + "or, for a constant fluid velocity, a vector of length 2 for a 2D problem holding this velocity" + + reference_velocity = 1.0 + @test_throws ArgumentError(error_str) BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + reference_density=0, + reference_pressure=0, + reference_velocity, + open_boundary_layers=2, + boundary_type=InFlow()) + + error_str = "`reference_pressure` must be either a function mapping " * + "each particle's coordinates and time to its pressure, " * + "a vector holding the pressure of each particle, or a scalar" + + reference_pressure = [1.0, 1.0] + + @test_throws ArgumentError(error_str) BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + reference_density=0, + reference_velocity=[1.0, + 1.0], reference_pressure, + open_boundary_layers=2, + boundary_type=InFlow()) + + error_str = "`reference_density` must be either a function mapping " * + "each particle's coordinates and time to its density, " * + "a vector holding the density of each particle, or a scalar" + + reference_density = [1.0, 1.0] + @test_throws ArgumentError(error_str) BoundaryZone(; plane, particle_spacing=0.1, + plane_normal=flow_direction, + density=1.0, + reference_density, + reference_velocity=[1.0, + 1.0], + reference_pressure=0, + open_boundary_layers=2, + boundary_type=InFlow()) + end @testset verbose=true "Boundary Zone 2D" begin particle_spacing = 0.2 open_boundary_layers = 4 @@ -35,8 +126,8 @@ zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing - sign_ = (first(typeof(boundary_zone).parameters) === - TrixiParticles.InFlow) ? -1 : 1 + sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? + -1 : 1 @test plane_points_1[i] == boundary_zone.zone_origin @test plane_points_2[i] - boundary_zone.zone_origin == @@ -99,8 +190,8 @@ zone_width = open_boundary_layers * boundary_zone.initial_condition.particle_spacing - sign_ = (first(typeof(boundary_zone).parameters) === - TrixiParticles.InFlow) ? -1 : 1 + sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? + -1 : 1 @test plane_points_1[i] == boundary_zone.zone_origin @test plane_points_2[i] - boundary_zone.zone_origin == @@ -137,7 +228,7 @@ @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones - perturb_ = first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow ? + perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ? sqrt(eps()) : -sqrt(eps()) @@ -183,7 +274,7 @@ @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones - perturb_ = first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow ? + perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ? eps() : -eps() point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 570818f8f..829e06c01 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -15,7 +15,8 @@ # Prescribed quantities reference_velocity = (pos, t) -> SVector(t, 0.0) reference_pressure = (pos, t) -> 50_000.0 * t - reference_density = (pos, t) -> 1000.0 * t + # Add small offset to avoid "ArgumentError: density must be positive and larger than `eps()`" + reference_density = (pos, t) -> 1000.0 * (t + sqrt(eps())) # Plane points of open boundary plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] @@ -36,10 +37,12 @@ flow_direction = flow_directions[j] inflow = BoundaryZone(; plane=plane_points, particle_spacing, density, plane_normal=flow_direction, open_boundary_layers, - boundary_type=InFlow()) + boundary_type=InFlow(), reference_velocity, + reference_pressure, reference_density) outflow = BoundaryZone(; plane=plane_points, particle_spacing, density, plane_normal=(-flow_direction), open_boundary_layers, - boundary_type=OutFlow()) + boundary_type=OutFlow(), reference_velocity, + reference_pressure, reference_density) boundary_zones = [ inflow, @@ -49,7 +52,7 @@ @testset "`$(TrixiParticles.boundary_type_name(boundary_zone))`" for boundary_zone in boundary_zones - sign_ = (first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow) ? + sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ? 1 : -1 fluid = extrude_geometry(plane_points; particle_spacing, n_extrude=4, density, pressure, @@ -62,10 +65,7 @@ boundary_system = OpenBoundarySPHSystem(boundary_zone; fluid_system, buffer_size=0, - boundary_model=BoundaryModelLastiwka(), - reference_velocity, - reference_pressure, - reference_density) + boundary_model=BoundaryModelLastiwkaCharacteristics()) semi = Semidiscretization(fluid_system, boundary_system) @@ -101,13 +101,13 @@ v, u, v0_ode, u0_ode, semi, t1) evaluated_vars1 = boundary_system.cache.characteristics - if first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow + if TrixiParticles.boundary_type_name(boundary_zone) == "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 first(typeof(boundary_zone).parameters) === TrixiParticles.OutFlow + elseif TrixiParticles.boundary_type_name(boundary_zone) == "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)) @@ -121,13 +121,13 @@ v, u, v0_ode, u0_ode, semi, t2) evaluated_vars2 = boundary_system.cache.characteristics - if first(typeof(boundary_zone).parameters) === TrixiParticles.InFlow + if TrixiParticles.boundary_type_name(boundary_zone) == "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 first(typeof(boundary_zone).parameters) === TrixiParticles.OutFlow + elseif TrixiParticles.boundary_type_name(boundary_zone) == "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/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 6e3fb9cd6..d50bbc25f 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -61,11 +61,12 @@ open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelTafuniMirroring(), buffer_size=0) semi = Semidiscretization(fluid_system, open_boundary) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary, semi) v_open_boundary = zero(inflow.initial_condition.velocity) v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') @@ -75,9 +76,7 @@ TrixiParticles.extrapolate_values!(open_boundary, FirstOrderMirroring(), v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false, - prescribed_velocity=false) + domain_fluid.coordinates, semi) # Checked visually in ParaView: # trixi2vtk(fluid_system.initial_condition, filename="fluid", # v=domain_fluid.velocity, p=domain_fluid.pressure) @@ -158,11 +157,12 @@ open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelTafuniMirroring(), buffer_size=0) semi = Semidiscretization(fluid_system, open_boundary) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary, semi) v_open_boundary = zero(inflow.initial_condition.velocity) v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') @@ -172,9 +172,7 @@ TrixiParticles.extrapolate_values!(open_boundary, FirstOrderMirroring(), v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false, - prescribed_velocity=false) + domain_fluid.coordinates, semi) # Checked visually in ParaView: # trixi2vtk(fluid_system.initial_condition, filename="fluid", # v=domain_fluid.velocity, p=domain_fluid.pressure) @@ -231,11 +229,12 @@ open_boundary_layers=open_boundary_layers, density=1000.0, particle_spacing, average_inflow_velocity=true) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelTafuniMirroring(), buffer_size=0) semi = Semidiscretization(fluid_system, open_boundary_in) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary_in, semi) v_open_boundary = zero(inflow.initial_condition.velocity) u_open_boundary = inflow.initial_condition.coordinates @@ -246,10 +245,10 @@ TrixiParticles.extrapolate_values!(open_boundary_in, FirstOrderMirroring(), v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0) + domain_fluid.coordinates, semi) TrixiParticles.average_velocity!(v_open_boundary, u_open_boundary, open_boundary_in, - inflow, semi) + first(open_boundary_in.boundary_zones), semi) # Since the velocity profile increases linearly in positive x-direction, # we can use the first velocity entry as a representative value. @@ -282,12 +281,13 @@ open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelTafuniMirroring(), buffer_size=0) # Temporary semidiscretization just to extrapolate the pressure into the outflow system semi = Semidiscretization(fluid_system, open_boundary_out) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary_out, semi) v_open_boundary = zero(outflow.initial_condition.velocity) v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') @@ -297,8 +297,7 @@ TrixiParticles.extrapolate_values!(open_boundary_out, mirror_method, v_open_boundary, v_fluid, outflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false) + domain_fluid.coordinates, semi) plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) @@ -306,12 +305,13 @@ plane_normal=[1.0, 0.0], open_boundary_layers=10, density=1000.0, particle_spacing) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=BoundaryModelTafuni(), + boundary_model=BoundaryModelTafuniMirroring(), buffer_size=0) # Temporary semidiscretization just to extrapolate the pressure into the outflow system semi = Semidiscretization(fluid_system, open_boundary_in) TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize!(open_boundary_in, semi) v_open_boundary = zero(inflow.initial_condition.velocity) @@ -320,8 +320,7 @@ TrixiParticles.extrapolate_values!(open_boundary_in, mirror_method, v_open_boundary, v_fluid, inflow.initial_condition.coordinates, - domain_fluid.coordinates, semi, 0.0; - prescribed_pressure=false) + domain_fluid.coordinates, semi) return fluid_system, open_boundary_in, open_boundary_out, v_fluid end @@ -335,7 +334,7 @@ v_fluid = mirror(pressure_func, mirror_method) p_fluid = [TrixiParticles.current_pressure(v_fluid, fluid_system, particle) - for particle in TrixiParticles.active_particles(fluid_system)] + for particle in TrixiParticles.eachparticle(fluid_system)] fluid_system.initial_condition.pressure .= p_fluid open_boundary_in.initial_condition.pressure .= open_boundary_in.pressure diff --git a/test/setups/extrude_geometry.jl b/test/setups/extrude_geometry.jl index e695f2e17..5a927dbee 100644 --- a/test/setups/extrude_geometry.jl +++ b/test/setups/extrude_geometry.jl @@ -28,7 +28,8 @@ ] @testset "Direction $i" for i in eachindex(directions) - shape = extrude_geometry((point1, point2); direction=directions[i], tlsph=true, + shape = extrude_geometry((point1, point2); direction=directions[i], + place_on_shell=true, particle_spacing=0.15, n_extrude=5, density=1.0) @test shape.coordinates ≈ expected_coords[i] @@ -68,7 +69,7 @@ end @testset "Direction $i" for i in eachindex(directions) shape = extrude_geometry(geometry; direction=directions[i], particle_spacing, - n_extrude=5, tlsph=true, density=1.0) + n_extrude=5, place_on_shell=true, density=1.0) @test shape.coordinates ≈ expected_coords[i] end @@ -86,7 +87,7 @@ end 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.06666666666666667 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.13333333333333333 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.26666666666666666 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.3333333333333333 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.4666666666666667 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.5333333333333333 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.6666666666666666 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.7333333333333333 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.8666666666666667 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 0.9333333333333333 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.092709445701687 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.15937611236835367 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.22604277903502035 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.292709445701687 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.35937611236835365 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.4260427790350203 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.492709445701687 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.5593761123683537 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.6260427790350204 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.692709445701687 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.7593761123683537 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8260427790350203 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.8927094457016871 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 0.9593761123683537 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.0260427790350204 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 1.092709445701687 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.185418891403374 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.2520855580700407 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.31875222473670733 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.38541889140337404 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.4520855580700407 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.5187522247367073 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.585418891403374 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.6520855580700406 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.7187522247367073 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.785418891403374 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.8520855580700406 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.9187522247367073 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 0.985418891403374 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.0520855580700408 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.1187522247367074 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 1.185418891403374 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.278128337105061 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.34479500377172767 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.4114616704383943 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.47812833710506103 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.5447950037717277 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.6114616704383944 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.678128337105061 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.7447950037717277 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.8114616704383943 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.878128337105061 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 0.9447950037717276 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0114616704383943 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.0781283371050612 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.1447950037717276 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.2114616704383945 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061 1.278128337105061] shape = extrude_geometry((p1, p2, p3); direction, particle_spacing=0.1, n_extrude=4, - density=1000.0, tlsph=true) + density=1000.0, place_on_shell=true) @test shape.coordinates ≈ expected_coords end diff --git a/test/setups/rectangular_shape.jl b/test/setups/rectangular_shape.jl index f3e77043b..7c2fd6875 100644 --- a/test/setups/rectangular_shape.jl +++ b/test/setups/rectangular_shape.jl @@ -41,7 +41,8 @@ ] @testset "$(loop_orders[i])" for i in eachindex(loop_orders) - shape = RectangularShape(1.0, (2, 2), (0.0, 0.0), density=1.0, tlsph=true, + shape = RectangularShape(1.0, (2, 2), (0.0, 0.0), density=1.0, + place_on_shell=true, loop_order=loop_orders[i]) @test shape.coordinates == expected_coords[i] @@ -242,7 +243,7 @@ end @testset "$(loop_orders[i])" for i in eachindex(loop_orders) shape = RectangularShape(1.0, (2, 2, 2), (0.0, 0.0, 0.0), density=1.0, - tlsph=true, loop_order=loop_orders[i]) + place_on_shell=true, loop_order=loop_orders[i]) @test shape.coordinates == expected_coords[i] end diff --git a/test/setups/sphere_shape.jl b/test/setups/sphere_shape.jl index 57d904a37..c8b87d76d 100644 --- a/test/setups/sphere_shape.jl +++ b/test/setups/sphere_shape.jl @@ -106,14 +106,14 @@ SphereShape(1.0, 1.1, (0.2, -1.0, 0.3), 1000.0, sphere_type=RoundSphere()), SphereShape(1.0, 1.2, (-0.3, 0.1, 0.8), 1000.0, sphere_type=RoundSphere()), SphereShape(0.1, 0.5, (0.3, 0.4, 0.5), 1000.0, cutout_min=(0.18, 0.4, 0.5), - cutout_max=(0.42, 10.0, 1.0), tlsph=true), - SphereShape(0.1, 0.5, (0.3, 0.4, 0.5), 1000.0, n_layers=2, tlsph=true), + cutout_max=(0.42, 10.0, 1.0), place_on_shell=true), + SphereShape(0.1, 0.5, (0.3, 0.4, 0.5), 1000.0, n_layers=2, place_on_shell=true), SphereShape(0.1, 0.5, (0.3, 0.4, 0.5), 1000.0, n_layers=2, - layer_outwards=true, tlsph=true), + layer_outwards=true, place_on_shell=true), SphereShape(0.1, 0.5, (0.3, 0.4, 0.5), 1000.0, n_layers=2, sphere_type=RoundSphere()), SphereShape(0.1, 0.55, (0.3, 0.4, 0.5), 1000.0, n_layers=2, layer_outwards=true, - sphere_type=RoundSphere(), tlsph=true) + sphere_type=RoundSphere(), place_on_shell=true) ] expected_coords = [ diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 3d8cf9bab..ed37b9ea5 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -32,7 +32,7 @@ @test system.mass == mass @test system.smoothing_kernel == smoothing_kernel @test TrixiParticles.initial_smoothing_length(system) == smoothing_length - @test system.transport_velocity isa Nothing + @test system.shifting_technique isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -87,7 +87,7 @@ @test system.mass == setup.mass @test system.smoothing_kernel == smoothing_kernel @test TrixiParticles.initial_smoothing_length(system) == smoothing_length - @test system.transport_velocity isa Nothing + @test system.shifting_technique isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -140,7 +140,7 @@ │ viscosity: …………………………………………………… Nothing │ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ - │ tansport velocity formulation: Nothing │ + │ shifting technique: …………………………… nothing │ │ average pressure reduction: ……… no │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ surface tension: …………………………………… nothing │ @@ -225,7 +225,7 @@ names = ["No TVF", "TransportVelocityAdami"] @testset "$(names[i])" for i in eachindex(transport_velocity) system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, - transport_velocity=transport_velocity[i], + shifting_technique=transport_velocity[i], average_pressure_reduction=true, smoothing_length, 0.0) semi = Semidiscretization(system) diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index 0fdf98f90..a31b5b366 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -1,71 +1,19 @@ @testset verbose=true "`OpenBoundarySPHSystem`" begin - @testset verbose=true "Illegal Inputs" begin - plane = ([0.0, 0.0], [0.0, 1.0]) - flow_direction = (1.0, 0.0) + @testset "`show`" begin # Mock fluid system struct FluidSystemMock2 <: TrixiParticles.FluidSystem{2} end TrixiParticles.initial_smoothing_length(system::FluidSystemMock2) = 1.0 TrixiParticles.nparticles(system::FluidSystemMock2) = 1 - 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, " * - "an array where the ``i``-th column holds the velocity of particle ``i`` " * - "or, for a constant fluid velocity, a vector of length 2 for a 2D problem holding this velocity" - - reference_velocity = 1.0 - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; - boundary_model=BoundaryModelLastiwka(), - buffer_size=0, - fluid_system=FluidSystemMock2(), - reference_density=0, - reference_pressure=0, - reference_velocity) - - error_str = "`reference_pressure` must be either a function mapping " * - "each particle's coordinates and time to its pressure, " * - "a vector holding the pressure of each particle, or a scalar" - - reference_pressure = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; - boundary_model=BoundaryModelLastiwka(), - buffer_size=0, - fluid_system=FluidSystemMock2(), - reference_density=0, - reference_velocity=[1.0, - 1.0], - reference_pressure) - - error_str = "`reference_density` must be either a function mapping " * - "each particle's coordinates and time to its density, " * - "a vector holding the density of each particle, or a scalar" - - reference_density = [1.0, 1.0] - @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; - boundary_model=BoundaryModelLastiwka(), - buffer_size=0, - fluid_system=FluidSystemMock2(), - reference_density, - reference_velocity=[1.0, - 1.0], - reference_pressure=0) - end - @testset "`show`" begin 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, - reference_pressure=0.0, - reference_velocity=[0.0, 0.0], + boundary_model=BoundaryModelLastiwkaCharacteristics(), fluid_system=FluidSystemMock2()) - show_compact = "OpenBoundarySPHSystem{2}(InFlow) with 80 particles" + show_compact = "OpenBoundarySPHSystem{2}() with 80 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -73,28 +21,21 @@ │ ════════════════════════ │ │ #particles: ………………………………………………… 80 │ │ #buffer_particles: ……………………………… 0 │ + │ #boundary_zones: …………………………………… 1 │ │ fluid system: …………………………………………… FluidSystemMock2 │ - │ boundary model: ……………………………………… BoundaryModelLastiwka │ - │ boundary type: ………………………………………… InFlow │ - │ prescribed velocity: ………………………… constant_vector │ - │ prescribed pressure: ………………………… constant_scalar │ - │ prescribed density: …………………………… constant_scalar │ - │ width: ……………………………………………………………… 0.2 │ + │ boundary model: ……………………………………… BoundaryModelLastiwkaCharacteristics │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box - outflow = BoundaryZone(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + outflow = BoundaryZone(; plane=([5.0, 0.0], [5.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, - reference_pressure=0.0, - reference_velocity=[0.0, 0.0], + boundary_model=BoundaryModelTafuniMirroring(), fluid_system=FluidSystemMock2()) - show_compact = "OpenBoundarySPHSystem{2}(OutFlow) with 80 particles" + show_compact = "OpenBoundarySPHSystem{2}() with 80 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -102,13 +43,28 @@ │ ════════════════════════ │ │ #particles: ………………………………………………… 80 │ │ #buffer_particles: ……………………………… 0 │ + │ #boundary_zones: …………………………………… 1 │ + │ fluid system: …………………………………………… FluidSystemMock2 │ + │ boundary model: ……………………………………… BoundaryModelTafuniMirroring │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", system) == show_box + + system = OpenBoundarySPHSystem(outflow, inflow; buffer_size=0, + boundary_model=BoundaryModelTafuniMirroring(), + fluid_system=FluidSystemMock2()) + + show_compact = "OpenBoundarySPHSystem{2}() with 160 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 160 │ + │ #buffer_particles: ……………………………… 0 │ + │ #boundary_zones: …………………………………… 2 │ │ fluid system: …………………………………………… FluidSystemMock2 │ - │ boundary model: ……………………………………… BoundaryModelLastiwka │ - │ boundary type: ………………………………………… OutFlow │ - │ prescribed velocity: ………………………… constant_vector │ - │ prescribed pressure: ………………………… constant_scalar │ - │ prescribed density: …………………………… constant_scalar │ - │ width: ……………………………………………………………… 0.2 │ + │ boundary model: ……………………………………… BoundaryModelTafuniMirroring │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box diff --git a/test/systems/packing_system.jl b/test/systems/packing_system.jl index d55afd559..b42e53427 100644 --- a/test/systems/packing_system.jl +++ b/test/systems/packing_system.jl @@ -19,7 +19,7 @@ │ neighborhood search: ………………………… GridNeighborhoodSearch │ │ #particles: ………………………………………………… 307 │ │ smoothing kernel: ………………………………… SchoenbergQuinticSplineKernel │ - │ tlsph: ……………………………………………………………… no │ + │ place_on_shell: ……………………………………… no │ │ boundary: ……………………………………………………… no │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box @@ -36,7 +36,7 @@ │ neighborhood search: ………………………… GridNeighborhoodSearch │ │ #particles: ………………………………………………… 307 │ │ smoothing kernel: ………………………………… SchoenbergQuinticSplineKernel │ - │ tlsph: ……………………………………………………………… no │ + │ place_on_shell: ……………………………………… no │ │ boundary: ……………………………………………………… yes │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box @@ -52,7 +52,7 @@ │ neighborhood search: ………………………… Nothing │ │ #particles: ………………………………………………… 307 │ │ smoothing kernel: ………………………………… SchoenbergQuinticSplineKernel │ - │ tlsph: ……………………………………………………………… no │ + │ place_on_shell: ……………………………………… no │ │ boundary: ……………………………………………………… no │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" end diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index cf6020d06..99f7152db 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -199,7 +199,7 @@ smoothing_length, density_diffusion=density_diffusion) - 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" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -211,8 +211,8 @@ │ state equation: ……………………………………… Val │ │ smoothing kernel: ………………………………… Val │ │ viscosity: …………………………………………………… nothing │ - │ tansport velocity formulation: Nothing │ │ density diffusion: ……………………………… Val{:density_diffusion}() │ + │ shifting technique: …………………………… nothing │ │ surface tension: …………………………………… nothing │ │ surface normal method: …………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │