fast-poisson.h

Constant density projection method

Constant density projection method

For the constant density method we need to solve the constant coefficient Poisson equation: \displaystyle \nabla ^2 p^{n+1}= \nabla \cdot \left[\left(1- \frac{\rho_0}{\rho^{n+1}}\right)\nabla \hat{p}\right]+ \frac{\rho_0}{\Delta t} \nabla \cdot u^{*} After solving the Poisson equation, the velocity field will be updated based on the new pressure field \displaystyle u^{n+1} = u^{*}-\Delta t \left[1/\rho_0 \nabla p^{n+1}+ \left(1/\rho^{n+1}-1/\rho_0\right)\nabla \hat{p}\right] Notations are based on the paper of Dodd and Ferrante, 2014.

mgstats fast_project (int method, face vector u, scalar p, scalar pn,
		      (const) face vector alpha, double rho0, double dt)
{

We initialize \hat{p}=pn for fastpn method and \hat{p}=2 p^n - p^{n-1} for fastpstar method.

  scalar phat[];
  switch(method){
  case 1: {
    foreach()
      phat[] = p[];
    boundary ({phat});
    break;
  }
  case 2: {
    foreach() {
#if 1
      phat[] = 1.5*p[] - 0.5*pn[];
#else // this was suggested by Michael Dodd but does not seem to work
      phat[] = 2.*p[] - pn[];
#endif
      pn[] = p[];
    }
    boundary ({phat, pn});
    break;
  }
  }

We calculate the term \frac{\rho_0}{\Delta t} \nabla \cdot u^{*} and store it in div_ustar variable.

  scalar div_ustar[];
  foreach() {
    div_ustar[] = 0.;
    foreach_dimension()
      div_ustar[] += u.x[1] - u.x[];
    div_ustar[] /= dt*cm[]*Delta; 
    div_ustar[] *= rho0;
  }
  boundary ({div_ustar});

We calculate the term \nabla \cdot \left[\left(1-\frac{\rho_0}{\rho^{n+1}}\right)\nabla \hat{p}\right] and store it in div_phat variable.

  face vector f[], coef[];
  scalar div_phat[];
  foreach_face()
    coef.x[]= (1. - rho0*alpha.x[]);
#if TREE
  foreach_face()
    f.x[] = coef.x[]*(phat[] - phat[-1])/(Delta);
  boundary_flux ({f});
  foreach() {
    div_phat[] = 0.;
    foreach_dimension()
      div_phat[] += f.x[1] - f.x[];
    div_phat[] /= cm[]*Delta;
  }
#else // Cartesian
  foreach() {
    div_phat[] = 0.;
    foreach_dimension()
      div_phat[] += (coef.x[1]*(phat[1] - phat[]) -
		     coef.x[]*(phat[] - phat[-1]));
    div_phat[] /= cm[]*sq(Delta);
  }
#endif
  boundary ({div_phat});

We construct the right hand side of Poisson equation and store it in variable b.

  scalar b[];
  foreach()
    b[] = div_phat[] + div_ustar[];
  boundary ({b});

We solve the Poisson problem. The tolerance (set with TOLERANCE) is the maximum relative change in volume of a cell (due to the divergence of the flow) during one timestep i.e. the non-dimensional quantity \displaystyle |\nabla\cdot\mathbf{u}|\Delta t Given the scaling of the divergence above, this gives

  mgstats mgp = poisson (p, b, tolerance = TOLERANCE/sq(dt));

And compute \mathbf{u}^{n+1} using \mathbf{u}^*, p, $p^{hat}, $

  foreach_face()
    u.x[] -= dt*(1./rho0*(p[] - p[-1]) +
		 (alpha.x[] - 1./rho0)*(phat[] - phat[-1]))/Delta;
  boundary ((scalar *){u});   
 
  return mgp;
}