Add some detail about making the MAC velocities

Creating the MAC velocities

~~~~~~~~~~~~~~~~~~~~~~~~~~~

To create the normal velocities on faces, we first extrapolate from the cell centers on each side and

upwind the face value to define  :math:`U^{pred}`

.  To compute the x-velocity on the x-faces of regular (ie not cut) cells, we call

   .. code:: shell

            AMREX_CUDA_HOST_DEVICE_FOR_3D(ubx, i, j, k,

             {

                 // X-faces

                 Real upls     = ccvel_fab(i  ,j,k,0) - 0.5 * xslopes_fab(i  ,j,k,0);

                 Real umns     = ccvel_fab(i-1,j,k,0) + 0.5 * xslopes_fab(i-1,j,k,0);

                 if ( umns < 0.0 && upls > 0.0 ) {

                    umac_fab(i,j,k) = 0.0;

                 } else {

                    Real avg = 0.5 * ( upls + umns );

                    if ( std::abs(avg) <  small_vel) { umac_fab(i,j,k) = 0.0;

                    } else if (avg >= 0)             { umac_fab(i,j,k) = umns;

                    } else                           { umac_fab(i,j,k) = upls;

                    }

                 }

             });

For cut cells we test on whether the area fraction is non-zero: 

   .. code:: shell

             AMREX_CUDA_HOST_DEVICE_FOR_3D(ubx, i, j, k,

             {

                 // X-faces

                 if (ax_fab(i,j,k) > 0.0)

                 {

                    Real upls     = ccvel_fab(i  ,j,k,0) - 0.5 * xslopes_fab(i  ,j,k,0);

                    Real umns     = ccvel_fab(i-1,j,k,0) + 0.5 * xslopes_fab(i-1,j,k,0);

                    if ( umns < 0.0 && upls > 0.0 ) {

                       umac_fab(i,j,k) = 0.0;

                    } else {

                       Real avg = 0.5 * ( upls + umns );

                       if ( std::abs(avg) <  small_vel) { umac_fab(i,j,k) = 0.0;

                       } else if (avg >= 0)             { umac_fab(i,j,k) = umns;

                       } else                           { umac_fab(i,j,k) = upls;

                       }

                    }

                 } else {

                       umac_fab(i,j,k) = huge_vel;

                 }

             });

We then perform a MAC projection on the face-centered velocities to enforce that they satisfy 

.. math:: \nabla \cdot (\varepsilon_g  U^{MAC})  = 0

We do this by solving 

.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi^{MAC} = \nabla \cdot \left( \varepsilon_g  U^{pred} \right)

then defining

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

Time Discretization

===================

In the absence of reactions, we assume that the fluid density is unchanged.

We compute the fluid volume fraction directly from the particle locations.

.. toctree::

   :maxdepth: 1

Thus here we focus on the discretization of the momentum equation

In the predictor

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

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

.. math:: (\varepsilon_g \rho_g U)^{\ast} &= (\varepsilon_g \rho_g U)^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

.. 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)

then defining

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

and 

.. math:: {p_g}^{n+1/2, \ast} = \phi

In the corrector

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

-  Define a new approximation to the new-time state, :math:`(\varepsilon_g \rho_g U)^{\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)

-  Project :math:`U^{\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)

then defining

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

and 

.. math:: {p_g}^{n+1/2} = \phi

   FluidTimeStep

   BuildingMacVelocities

Fluid Time Step

~~~~~~~~~~~~~~~

In the absence of reactions, we assume that the fluid density is unchanged.

We compute the fluid volume fraction directly from the particle locations.

Thus here we focus on the discretization of the momentum equation

In the predictor

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

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

.. math:: (\varepsilon_g \rho_g U)^{\ast} &= (\varepsilon_g \rho_g U)^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

.. 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)

then defining

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

and 

.. math:: {p_g}^{n+1/2, \ast} = \phi

In the corrector

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

-  Define a new approximation to the new-time state, :math:`(\varepsilon_g \rho_g U)^{\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)

-  Project :math:`U^{\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)

then defining

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

and 

.. math:: {p_g}^{n+1/2} = \phi

===========================

.. toctree::

   :maxdepth: 1

   FluidEquations

   FluidTimeDiscretization

the master branch at the beginning of each month.

.. toctree::

   :maxdepth: 1

   :maxdepth: 0

   :caption: Contents:

   Introduction