Implement consistent shifting by Sun et al. 2019 (#888)

* Implement consistent shifting by Sun et al. 2019

* Fix tests

* Add option to switch between dynamic and static v_max

* Fix docs

* Add more constructor tests

* Fix `dv_shifting`

* Add reference

* Remove doctest output

* Fix ugly spaces

* Fix tests

* Fix tests

* Update NEWS.md

* Reformat

* Add #879 to NEWS.md

* Fix tests

* Revise docstring

* Make API more flexible

* Fix tests

Author: efaulhaber

NEWS.md

@@ -11,14 +11,28 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - API for `OpenBoundarySPHSystem` and `BoundaryZone` changed.
   It is now possible to pass multiple `BoundaryZone`s to a single `OpenBoundarySPHSystem`.
   Reference values are now assigned individually to each `BoundaryZone`. (#866)
+
+- The argument of `TransportVelocityAdami` is now a keyword argument.
+  `TransportVelocityAdami(1000.0)` now becomes
+  `TransportVelocityAdami(background_pressure=1000.0)` (#884).
+
 - 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.
+  instead, and pass `shifting_technique=ParticleShiftingTechniqueSun2017()` to the system (#884).
 
 - Renamed the keyword argument `tlsph` to `place_on_shell` for `ParticlePackingSystem`,
-  `sample_boundary`, `extrude_geometry`, `RectangularShape`, and `SphereShape`.
+  `sample_boundary`, `extrude_geometry`, `RectangularShape`, and `SphereShape` (#814).
+
+- Custom quantity functions passed to `SolutionSavingCallback` or `PostprocessCallback`
+  that were not using the documented API but were functions of
+  `(system, v_ode, u_ode, semi, t)` now need to be functions
+  of `(system, dv_ode, du_ode, v_ode, u_ode, semi, t)` (#879).
+
+### Features
+
+- Added consistent particle shifting by Sun et al. (2019) as `ConsistentShiftingSun2019` (#888).
 
 ## Version 0.3.1
 

docs/src/refs.bib

@@ -707,6 +707,17 @@ @article{Sun2018
 	doi       = {10.1016/j.cpc.2017.11.016},
 }
 
+@article{Sun2019,
+	author    = {Sun, P. N. and Colagrossi, A. and Marrone, S. and Antuono, M. and Zhang, A. -M.},
+	title     = {A consistent approach to particle shifting in the δPlus-{SPH} model},
+	journal   = {Computer Methods in Applied Mechanics and Engineering},
+	volume    = {348},
+	year      = {2019},
+	issn      = {0045-7825},
+	pages     = {912--934},
+	doi       = {10.1016/j.cma.2019.01.045},
+}
+
 @Article{Tafuni2018,
   author    = {A. Tafuni and J.M. Dom{\'{\i}}nguez and R. Vacondio and A.J.C. Crespo},
   journal   = {Computer Methods in Applied Mechanics and Engineering},

examples/fluid/lid_driven_cavity_2d.jl

@@ -65,15 +65,15 @@ if wcsph
                                                state_equation, smoothing_kernel,
                                                pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure,
                                                smoothing_length, viscosity=viscosity,
-                                               shifting_technique=TransportVelocityAdami(pressure))
+                                               shifting_technique=TransportVelocityAdami(background_pressure=pressure))
 else
     state_equation = nothing
     density_calculator = ContinuityDensity()
     fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel,
                                                smoothing_length,
                                                density_calculator=density_calculator,
                                                sound_speed, viscosity=viscosity,
-                                               shifting_technique=TransportVelocityAdami(pressure))
+                                               shifting_technique=TransportVelocityAdami(background_pressure=pressure))
 end
 
 # ==========================================================================================

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),
-                                           shifting_technique=TransportVelocityAdami(pressure),
+                                           shifting_technique=TransportVelocityAdami(background_pressure=pressure),
                                            acceleration=(acceleration_x, 0.0))
 
 # ==========================================================================================

examples/fluid/pipe_flow_2d.jl

@@ -83,7 +83,7 @@ if wcsph
                                                state_equation, smoothing_kernel,
                                                density_diffusion=density_diffusion,
                                                smoothing_length, viscosity=viscosity,
-                                               shifting_technique=ParticleShiftingTechnique(),
+                                               shifting_technique=ParticleShiftingTechnique(v_max_factor=1.5),
                                                buffer_size=n_buffer_particles)
 else
     # Alternatively the EDAC scheme can be used

examples/fluid/taylor_green_vortex_2d.jl

@@ -86,13 +86,15 @@ if wcsph
                                                pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure,
                                                smoothing_length,
                                                viscosity=ViscosityAdami(; nu),
-                                               shifting_technique=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,
-                                               shifting_technique=TransportVelocityAdami(background_pressure),
+                                               shifting_technique=TransportVelocityAdami(;
+                                                                                         background_pressure),
                                                viscosity=ViscosityAdami(; nu))
 end
 

src/TrixiParticles.jl

@@ -68,7 +68,9 @@ export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback
 export ContinuityDensity, SummationDensity
-export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique
+export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique,
+       ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019,
+       ContinuityEquationTermSun2019, MomentumEquationTermSun2019
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel

src/schemes/boundary/open_boundary/boundary_zones.jl

@@ -80,7 +80,7 @@ There are three ways to specify the actual shape of the boundary zone:
     or evolved using the characteristic flow variables ([BoundaryModelCharacteristicsLastiwka](@ref BoundaryModelCharacteristicsLastiwka)).
 
 # Examples
-```jldoctest; output = false
+```jldoctest; output=false
 # 2D
 plane_points = ([0.0, 0.0], [0.0, 1.0])
 plane_normal = [1.0, 0.0]

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -53,10 +53,10 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         # 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,
+                             particle_system, neighbor_system,
                              v_particle_system, v_neighbor_system,
-                             m_a, m_b, rho_a, rho_b, pos_diff, distance,
-                             grad_kernel, correction)
+                             particle, neighbor, 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,