Possible bug in TransportVelocityStep

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^2_j)\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])

Hi Stephen,

The stepper is correct.

See that the auhat is computed with a negative sign

pysph/pysph/sph/wc/transport_velocity.py

Line 311 in e66505c

tmp = -self.pb * mi1 * (Vi2 + Vj2)

here d_u is updated

pysph/pysph/sph/integrator_step.py

Line 273 in e66505c

d_u[d_idx] += dtb2*d_au[d_idx]

So when we add auhat in

pysph/pysph/sph/integrator_step.py

Line 278 in e66505c

d_uhat[d_idx] = d_u[d_idx] + dtb2*d_auhat[d_idx]

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.

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