Update U --> U_g and tau^*** instead of tau^** in corrector.

Conservation of fluid momentum:

.. math:: \frac{ \partial (\varepsilon_g \rho_g U)}{\partial t} + \nabla \cdot (\varepsilon_g \rho_g U_g U_g) + \varepsilon_g \nabla p_g = \nabla \cdot \tau

.. math:: \frac{ \partial (\varepsilon_g \rho_g U_g)}{\partial t} + \nabla \cdot (\varepsilon_g \rho_g U_g U_g) + \varepsilon_g \nabla p_g = \nabla \cdot \tau

           + \sum_p \beta_p (V_p - U_g) + \rho_g g

where :math:`\sum_p \beta_p (V_p - U_g)` is the drag term in which :math:`V_p` represents the particle velocity and :math:`\beta_p` is the drag coefficient associated with that particle

In the predictor

-  Define :math:`U^{MAC,n}`, the face-centered (staggered) MAC velocity which is used for advection, using :math:`U^n`

-  Define :math:`U^{MAC,n}`, the face-centered (staggered) MAC velocity which is used for advection, using :math:`U_g^n`

-  Define an approximation to the new-time state, :math:`(\varepsilon_g \rho_g U)^{\ast}` by setting 

-  Define an approximation to the new-time state, :math:`(\varepsilon_g \rho_g U_g)^{\ast}` by setting 

.. math:: (\varepsilon_g \rho_g U)^{\ast} &= (\varepsilon_g \rho_g U)^n -  

.. math:: (\varepsilon_g \rho_g U_g)^{\ast} &= (\varepsilon_g \rho_g U_g)^n -  

           \Delta t \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g) + \varepsilon_g \nabla {p_g}^{n-1/2} \right) \\ &+ 

           \Delta t \left( \nabla \cdot \tau^n + \sum_p \beta_p (V_p - {U_g}^{\ast}) + \rho_g \varepsilon_g g \right)

-  Project :math:`U^{\ast}` by solving

-  Project :math:`U_g^{\ast}` by solving

.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g  U)^{\ast}+ {\varepsilon_g}{\rho_g} \nabla {p_g}^{n-1/2} \right)

.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g  U_g)^{\ast}+ {\varepsilon_g}{\rho_g} \nabla {p_g}^{n-1/2} \right)

then defining

.. math:: U^{\ast \ast} = U^{\ast} - \frac{1}{\rho_g} \nabla \phi

.. math:: U_g^{\ast \ast} = U_g^{\ast} - \frac{1}{\rho_g} \nabla \phi

and 

In the corrector

-  Define :math:`U^{MAC,\ast \ast}` at the "new" time using :math:`U^{\ast \ast}`

-  Define :math:`U^{MAC,\ast \ast}` at the "new" time using :math:`U_g^{\ast \ast}`

-  Define a new approximation to the new-time state, :math:`(\varepsilon_g \rho_g U)^{\ast \ast \ast}` by setting  

-  Define a new approximation to the new-time state, :math:`(\varepsilon_g \rho_g U_g)^{\ast \ast \ast}` by setting  

.. math:: (\varepsilon_g \rho_g U)^{\ast \ast \ast} &= (\varepsilon_g \rho_g U)^n - \frac{\Delta t}{2} \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^n + \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^{\ast \ast}\right) + \\ &+ \frac{\Delta t}{2} \left( \nabla \cdot \tau^n + \nabla \cdot \tau^{\ast \ast} \right) + \Delta t \left( - \varepsilon_g \nabla {p_g}^{n+1/2,\ast} + \sum_p \beta_p (V_p - {U_g}^{\ast \ast \ast}) + \varepsilon_g \rho_g g \right)

.. math:: (\varepsilon_g \rho_g U_g)^{\ast \ast \ast} &= (\varepsilon_g \rho_g U_g)^n - \frac{\Delta t}{2} \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^n + \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^{\ast \ast}\right) + \\ &+ \frac{\Delta t}{2} \left( \nabla \cdot \tau^n + \nabla \cdot \tau^{\ast \ast \ast} \right) + \Delta t \left( - \varepsilon_g \nabla {p_g}^{n+1/2,\ast} + \sum_p \beta_p (V_p - {U_g}^{\ast \ast \ast}) + \varepsilon_g \rho_g g \right)

-  Project :math:`U^{\ast \ast \ast}` by solving

-  Project :math:`U_g^{\ast \ast \ast}` by solving

.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g  U)^{\ast \ast \ast} + \frac{\varepsilon_g}{\rho_g} \nabla {p_g}^{n+1/2,\ast} \right)

.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g  U_g)^{\ast \ast \ast} + \frac{\varepsilon_g}{\rho_g} \nabla {p_g}^{n+1/2,\ast} \right)

then defining

.. math:: U^{n+1} = U^{\ast \ast \ast} - \frac{1}{\rho_g} \nabla \phi

.. math:: U_g^{n+1} = U_g^{\ast \ast \ast} - \frac{1}{\rho_g} \nabla \phi

and