Time step in pressure projection + velocity correction

Issue #8 resolved
Antonio Alguacil created an issue

Hello, I'm working on deep learning applied to fluid mechanics following the approach of Tompson et al. - Accelerating Eulerian Fluid Simulation With Convolutional Networks. They generated their dataset using Mantaflow using a timestep of 0.1 for their solver and I have one question about the implementation of your code. In your plume_2d example you choose your timestep = 1, but what if I modify this value (e.g dt = 0.1). I've seen that this time step is naturally taken into account in the Semi-Lagrangian advection but NOT in the resolution of the Poisson equation (e.g. coefficient in the diagonal are 1 whereas they should be dt/ (rhodxdx) and knowing that dx=1 and rho=1 in all your code, the dt should still remain there.), NEITHER in the velocity correction (source/plugin/pressure.cpp, CorrectVelocity): there, you got

//! Kernel: Apply velocity update from poisson equation
KERNEL(bnd=1)
void CorrectVelocity(FlagGrid& flags, MACGrid& vel, Grid<Real>& pressure)
{
    int idx = flags.index(i,j,k);
    if (flags.isEmpty(idx)) {
      // This was originally a debug print, but now I'll leave it in here as a
      // way to make sure we're never marking flags as empty (otherwise our own
      // velocity update code will break). (i.e. We DO NOT handle the empty
      // cell case).
      printf("isEmpty (i, j, k) = (%d, %d, %d) - NOT EXPECTED!\n", i, j, k);
    }
    if (flags.isFluid(idx))
    {
        if (flags.isFluid(i-1,j,k)) vel[idx].x -= (pressure[idx] - pressure(i-1,j,k));
        if (flags.isFluid(i,j-1,k)) vel[idx].y -= (pressure[idx] - pressure(i,j-1,k));
        if (flags.is3D() && flags.isFluid(i,j,k-1)) vel[idx].z -= (pressure[idx] - pressure(i,j,k-1));

        if (flags.isEmpty(i-1,j,k)) vel[idx].x -= pressure[idx];
        if (flags.isEmpty(i,j-1,k)) vel[idx].y -= pressure[idx];
        if (flags.is3D() && flags.isEmpty(i,j,k-1)) vel[idx].z -= pressure[idx];
    }
    else if (flags.isEmpty(idx)&&!flags.isOutflow(idx)) // don't change velocities in outflow cells
    {
        if (flags.isFluid(i-1,j,k)) vel[idx].x += pressure(i-1,j,k);
        else                        vel[idx].x  = 0.f;
        if (flags.isFluid(i,j-1,k)) vel[idx].y += pressure(i,j-1,k);
        else                        vel[idx].y  = 0.f;
        if (flags.is3D() ) {
        if (flags.isFluid(i,j,k-1)) vel[idx].z += pressure(i,j,k-1);
        else                        vel[idx].z  = 0.f;
        }
    }
}

and I would have expected a dt multiplying the discrete pressure gradient (p(i,j,k) - p(i-1,j,k)). Could you confirm that you have make the assumption of dt=1? This timestep is used in the advection and buoyancy/gravity algorithms but not in pressure projection/velocity correction.

Thank you very much for your time! Antonio

Comments (1)

  1. Antonio Alguacil reporter

    I'm going to leave the response that I kindly got from Nils Thuerey: no time step multiplication is needed in the projection step as if you divide the rhs of the system by dt, you also multiply the gradient of pressure by the same amount, effectively cancelling both contributions. Therefore, it is better to leave it as in Mantaflow in order to obtain a well conditionned Laplacian matrix, Thank you very much for your work!!

  2. Log in to comment