Restructure TVF

[efaulhaber]

This PR consists of several changes that I did in steps.

Step 1

The current approach is integrating the transport velocity through stages. So, between time steps, \tilde{v} is set to v, and then it is integrated for one time step with d\tilde{v}/dt set to the term inside the parentheses here:

Note that this approach is currently flawed, as the dv/dt term in the RHS is missing, so only the physical velocity from the last time step is used at every stage.

In the paper, with the different time integration, this transport velocity is only computed once per step:

While we can interpret this as "transport velocity is integrated within the step starting from v", which leads to the approach explained above, we could also interpret this as "dr/dt = v + δv, where δv is dt * v_c and v_c is the second term inside the parentheses above".
This interpretation leads to a δv = dt * v_c that is either constant throughout the time step and computed between steps, or v_c is recomputed at every stage. The latter is the one I implemented for the following comparisons, but I don't think the difference between computing this only between steps is significant, so we might be able to save some significant computation here.

Results

As can be seen here in the left column, directly behind the cylinder, the old TVF implementation (top row) is less effective at preventing void regions than the new implementation (top row). The old implementation with p_0 = 100k (which is rho_0 * c^2) is about as effective in preventing void regions as the new version with p_0 = 50k, whereas the new version with p_0 = 100k is about as effective as the old version with p_0 = 200k. This indicates that the new implementation might just be a rescaling of the old implementation. We can see here, similar to the results before, that the top left and bottom right show similar particle distribution uniformity (mostly in the red area above the wake, but also in the size of the "reordered" region in front of the cylinder), as well as the top right and bottom left. However, we can now see in the left column that the old TVF implementation changes the solution, producing a slightly wider wake region and lower velocity in the center of the wake. The higher background pressure changes the results even more (bottom right), whereas the top row solution is not affected by the different background pressure. This is a comparison between TVF in this PR (top), consistent shifting from #843 (middle) and inconsistent shifting (bottom).

Step 2

The original 2013 paper was using kernel summation to compute the density. When we use this implementation with the continuity equation, we get the same inconsistency as the inconsistent PST, as particles are now moving with the transport velocity, but the physical velocity is used in the continuity equation.
This can lead to an uncontrolled increase of density in a closed simulation.

In "A generalized transport-velocity formulation for smoothed particle hydrodynamics" (Zhang et al. 2017), the density is integrated with the continuity equation and the transport velocity is used here:

[LasNikas]

very nice!

[codecov]

Codecov Report

❌ Patch coverage is 41.60000% with 73 lines in your changes missing coverage. Please review.
✅ Project coverage is 71.03%. Comparing base (995f182) to head (92e9e7a).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/schemes/fluid/transport_velocity.jl	18.36%	40 Missing ⚠️
src/general/semidiscretization.jl	52.38%	10 Missing ⚠️
src/preprocessing/particle_packing/system.jl	0.00%	6 Missing ⚠️
src/schemes/boundary/open_boundary/system.jl	33.33%	4 Missing ⚠️
...rc/schemes/fluid/entropically_damped_sph/system.jl	62.50%	3 Missing ⚠️
src/callbacks/update.jl	0.00%	2 Missing ⚠️
src/general/system.jl	80.00%	2 Missing ⚠️
src/callbacks/density_reinit.jl	0.00%	1 Missing ⚠️
src/callbacks/particle_shifting.jl	0.00%	1 Missing ⚠️
src/callbacks/post_process.jl	0.00%	1 Missing ⚠️
... and 3 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #864      +/-   ##
==========================================
+ Coverage   70.98%   71.03%   +0.04%     
==========================================
  Files         107      107              
  Lines        7045     7018      -27     
==========================================
- Hits         5001     4985      -16     
+ Misses       2044     2033      -11     

Flag	Coverage Δ
unit	71.03% <41.60%> (+0.04%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[efaulhaber]

Update: I compared the integrated implementation on main (with the bug fixed) with the new implementation. Here is the kinetic energy in the Taylor-Green vortex.

As expected, the new implementation yields a constant background pressure throughout different time integration methods and step sizes.
Also, the implementation on main with a CFL of 0.25 with the Leapfrog method (the red line shows an adaptive CFL condition, the purple line a fixed time step), which is what the method was designed for in the paper, is close to the gold and grey line, which is the new implementation. The turquoise line is slightly different, probably due to time integration errors.

[svchb]

Can you add the continuity equation based on the advection velocity? At least in my testing this was much better!

[efaulhaber]

This is step 2 that I described above. It will be the next PR. This PR doesn't change features, it only refactors the code and fixes a bug.

[efaulhaber]

/run-gpu-tests

[svchb]

/run-gpu-tests

[efaulhaber]

/run-gpu-tests

[svchb]

/run-gpu-tests