创建 MAC 速度

为了在面上创建法向速度，我们首先使用之前计算的斜率从每一侧的单元中心进行外推，并使用迎风法定义面值 \(U^{pred}\)。

为了在常规（即未切割）单元的 x 面上计算 x 方向速度，我们调用

AMREX_CUDA_HOST_DEVICE_FOR_3D(ubx, i, j, k,
 {
     // X-面
     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;
        }
     }
 });

对于切割单元，我们测试面积分数是否非零：

AMREX_CUDA_HOST_DEVICE_FOR_3D(ubx, i, j, k,
{
    // X-面
    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;
    }
});

然后我们对面心速度执行 MAC 投影，以确保它们满足

\[\nabla \cdot (\varepsilon_g U^{MAC}) = 0\]

我们通过求解以下方程实现这一点

\[\nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi^{MAC} = \nabla \cdot \left( \varepsilon_g U^{pred} \right)\]

然后定义

\[U^{MAC} = U^{pred} - \frac{1}{\rho_g} \nabla \phi^{MAC}\]