Possible bug in TransportVelocityStep

[Qcellaris]

Dear developers,

When we look at equation 13 in "A transport-velocity formulation for smoothed particle hydrodynamics", by Adami, Hu and Adams 2013, we see the term: $$\frac{\tilde d\mathbf{v}_i}{dt}-\frac{p_b}{m_i}\sum_j(V^2_i+V^2j)\frac{\partial W}{\partial r_{ij}}\mathbf{e}\{ij}$$ of which the first term is computed in equation 8. In the TVFScheme of pysph the result is stored in au, av and aw. The second term is stored in pysph as auhat, avhat and awhat. Then if we look at the time integration:

class TransportVelocityStep(IntegratorStep):
    def initialize(self):
        pass

    def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_uhat, d_auhat, d_vhat, d_avhat, d_what, d_awhat, d_x, d_y, d_z, dt):
        dtb2 = 0.5*dt

        # velocity update eqn (14)
        d_u[d_idx] += dtb2*d_au[d_idx]
        d_v[d_idx] += dtb2*d_av[d_idx]
        d_w[d_idx] += dtb2*d_aw[d_idx]

        # advection velocity update eqn (15)
        d_uhat[d_idx] = d_u[d_idx] + dtb2*d_auhat[d_idx]
        d_vhat[d_idx] = d_v[d_idx] + dtb2*d_avhat[d_idx]
        d_what[d_idx] = d_w[d_idx] + dtb2*d_awhat[d_idx]

        # position update eqn (16)
        d_x[d_idx] += dt*d_uhat[d_idx]
        d_y[d_idx] += dt*d_vhat[d_idx]
        d_z[d_idx] += dt*d_what[d_idx]

    def stage2(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_vmag2, dt):
        dtb2 = 0.5*dt

        # corrector update eqn (17)
        d_u[d_idx] += dtb2*d_au[d_idx]
        d_v[d_idx] += dtb2*d_av[d_idx]
        d_w[d_idx] += dtb2*d_aw[d_idx]

        # magnitude of velocity squared
        d_vmag2[d_idx] = (d_u[d_idx]*d_u[d_idx] + d_v[d_idx]*d_v[d_idx] + d_w[d_idx]*d_w[d_idx])

We see that in the update of the advection velocity auhat, avhat and awhat are used. If I look at the paper correctly the background pressure force of equation 15 consists of both terms so I think it should be (au-auhat), (av-avhat) and (aw-awhat), i.e.:

 # advection velocity update eqn (15)
        d_uhat[d_idx] = d_u[d_idx] + dtb2*(d_au[d_idx]-d_auhat[d_idx])
        d_vhat[d_idx] = d_v[d_idx] + dtb2*(d_av[d_idx]-d_avhat[d_idx])
        d_what[d_idx] = d_w[d_idx] + dtb2*(d_aw[d_idx]-d_awhat[d_idx])

[pawansnegi]

Hi Stephen,

The stepper is correct.

See that the auhat is computed with a negative sign

https://github.com/pypr/pysph/blob/e66505c4c4f1eba1b66fa2579a123b820f9f6781/pysph/sph/wc/transport_velocity.py#L311

here d_u is updated

https://github.com/pypr/pysph/blob/e66505c4c4f1eba1b66fa2579a123b820f9f6781/pysph/sph/integrator_step.py#L273

So when we add auhat in

https://github.com/pypr/pysph/blob/e66505c4c4f1eba1b66fa2579a123b820f9f6781/pysph/sph/integrator_step.py#L278

we have

d_u (old) + dtb2 * (au - auhat)

and the negative sign is already taken care in the equation.

I hope this clears the doubt.

[Qcellaris]

You are completely right, I didn't look carefully enough. Thanks for your reply!