Reopen #221: Transport Velocity for EDAC #436
Conversation
| @inline function transport_velocity!(dv, system, ::TransportVelocityAdami, | ||
| rho_a, rho_b, m_a, m_b, grad_kernel, particle) | ||
| (; transport_velocity) = system | ||
| (; background_pressure) = transport_velocity | ||
|
|
||
| NDIMS = ndims(system) | ||
|
|
||
| volume_a = m_a / rho_a | ||
| volume_b = m_b / rho_b | ||
| volume_term = (volume_a^2 + volume_b^2) / m_a | ||
|
|
||
| for dim in 1:NDIMS | ||
| dv[NDIMS + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] | ||
| end | ||
|
|
||
| return dv | ||
| end |
There was a problem hiding this comment.
I'm pretty sure this doesn't work as intended.
Say we're using the Euler method for time integration.
Then we have only one stage. And at the beginning of this stage, velocity and transport velocity are the same, since you just did that in the update callback.
Maybe it's time to ask Adami if there is a way to implement this for multi-stage time integration methods?
There was a problem hiding this comment.
It is an explicit scheme so any suitable integrator can be used, no?
I mean, when using the euler method, it is just a normal evaluation without TVF (without correction). It's not incorrect then, no?
On the other hand, when using a multi-stage integration scheme, we correct the momentum at each stage as it is supposed to be.
Nevertheless you're right. We should contact Adami.
There was a problem hiding this comment.
The correction is applied during the stages. For momentum_velocity = advection_velocity (single stage Euler method) the convection term is zero. That is, no correction due to transport velocity.
There was a problem hiding this comment.
So far I understood, the TVF only takes affect for integration schemes with 2 stages and more.
There was a problem hiding this comment.
The kick-drift-kick form used in Adami's paper is basically only one stage. It never uses the non-transport velocity for advection, only the transport velocity. Here, the first stage would use the regular velocity.
There was a problem hiding this comment.
It never uses the non-transport velocity for advection, only the transport velocity.
I'm not sure what you mean, but, yes, the position update always uses the transport velocity.
Here, the first stage would use the regular velocity.
Yes, at first stage, no correction occurs, since transport velocity is equal to momentum velocity. But during the stages, the correction is applied due to the transport velocity.
So to speak, the particles are forced to be on the supposed trajectory during the stages according to the background pressure.
Ramachandran (2019) is using a two-stage Evaluate-Predict-Evaluate-Correct (EPEC) scheme.
There was a problem hiding this comment.
This is from Ramachandran 2019:
Note that in the first stage, the transport velocity is used to advect the particles, while in your implementation, the regular velocity is used in the first stage.
Thank you for pointing me to this article. I've been thinking a lot about how to correctly implement this for multi-stage methods. So the correct solution is to just compute the transport velocity from the regular velocity at the beginning of the kick! as
In the first screenshot, you can see that dt/2 is used instead in the first stage, which is the "stage-local" step size used for that stage.
So we have two options:
- Just use dt everywhere. It doesn't really matter that it is bigger than the stage-local step size, as the background pressure is arbitrary anyway.
- Use
dt = t_stage - t_laststep. The main difference to option 1 is that this dt is getting bigger for later stages, which is the same as in the first screenshot above, so this is closer (identical?) to what Ramachandran is doing. We could either fiddle with a stage limiter (which has the disadvantage that now all methods support it), or we could just store the time after the last full time step with a callback ast_laststep.
So then the transport velocity is not integrated, but locally computed at the beginning of kick! with the equation in the second screenshot and dt = t - t_laststep. This will also work for Euler.
There was a problem hiding this comment.
As discussed in zoom, the challenge is that
It is not as easy as it seems (or impossible?) with our time integrator approach. So far, the current implementation of this PR is the best solution or rather workaround.
However, @efaulhaber will implement a more flexible time integrator approach, right? Does an issue already exist or will you create one?
There was a problem hiding this comment.
So we wait with this PR till this is finished?
There was a problem hiding this comment.
No, I would say that the current implementation is sufficient (according to the validation results). That is, we also mark this feature as experimental implementation.
However, it would be interesting to compare the results with appropriated time integration schemes.
| pressure_average[particle] /= neighbor_counter[particle] | ||
| end | ||
| end | ||
| pressure_average ./= neighbor_counter |
There was a problem hiding this comment.
Why is it okay here if we divide by zero?
It's because if the particle doesn't have any neighbors, then the average pressure value is never used, right?
Can you please add a comment?
There was a problem hiding this comment.
We already discussed this in a resolve thread: #436 (comment)
There was a problem hiding this comment.
Its because the neighbor_counter is number of neighbors + 1
| fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, | ||
| sound_speed, | ||
| transport_velocity=TransportVelocityAdami(background_pressure), | ||
| viscosity=ViscosityAdami(; nu)) |
There was a problem hiding this comment.
| viscosity=ViscosityAdami(; nu)) | |
| viscosity=ViscosityAdami(; nu)) |
validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl
Outdated
Show resolved
Hide resolved
validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl
Outdated
Show resolved
Hide resolved
| @@ -0,0 +1,93 @@ | |||
| using TrixiParticles | |||
There was a problem hiding this comment.
Are these added to the tests? I didn't see them.
There was a problem hiding this comment.
I need #601 for this example
see #221
based on #440