Roe flux

[pjb236]

Add Roe flux capability.

Following chapter 11 of "Riemann Solvers and Numerical Methods" by Toro, the algorithm should

1. Compute Roe averaged velocities, total enthalpy, and speed of sound
2. Compute Roe averaged eigenvalues using an analytical expression, potentially with an "entropy fix".
3. Compute the Roe averaged right eigenvectors using analytical expressions for each entry
4. Project the jump in each conserved quantity onto the right eigenvectors. Again, there are analytical expressions for this in Toro's book or Roe's original paper.
5. Compute fluxes using equations 11.27-11.29 in Toro.

All of this should be possible following the same API used for the Rusanov flux in pressiodemoapps/impl_apps/eulerCommon/fluxes.hpp

Important note: there is a typo in Roe's original paper! The right hand side of Equation 22b should be \delta u_4 - w * \delta u_1.

[pjb236]

@fnrizzi I have 1d Roe flux working: it was able to solve the Sod Shock tube correctly. I've added a test for the flux function, but we don't have a test where we solve a 1D Euler problem with Roe Flux. Should I add a new test or modify one of the existing tests, such as the 1D sod cases, to use the Roe flux instead? Perhaps I could add separate sub-directories in 1d_sod for Roe flux tests like the weno3 and weno5 subdirectories?

I'll start adding 2D and 3D Roe fluxes tomorrow.

[pjb236]

@fnrizzi Good news: I've implemented 2D Roe flux and it works for the 2D normal shock and 2D double Mach reflection with WENO3. Bad news: I can't get it to work for WENO5. The branch is up to date if you want to try, it is set up to make a WENO5 test for the double Mach reflection that uses the Roe Flux. Looks like some of the boundaries are unstable after one timestep.

[fnrizzi]

template<class T, typename sc_t>
void eeRoeFluxThreeDof(T & F,
               const T & qL,
               const T & qR,
               const sc_t gamma)
{
  constexpr auto one  = static_cast<sc_t>(1);
  constexpr auto two  = static_cast<sc_t>(2);
  constexpr auto half = one/two;

  constexpr sc_t es_ef{0.1}; // Entropy fix parameter
  constexpr sc_t es{1.e-30};
  sc_t FL[3];
  sc_t FR[3];
  sc_t K1[3];
  sc_t K2[3];
  sc_t K3[3];

  // Compute Fluxes  
  // left
  const auto rL = qL(0);
  const auto uL = qL(1)/(rL + es);
  const auto pL = (gamma-1)*( qL(2) - half*rL*(uL*uL) );
  const auto HL = ( qL(2) + pL ) / rL;
  FL[0] = rL*uL;
  FL[1] = rL*uL*uL + pL;
  FL[2] = rL*uL*HL;

  // right
  const auto rR = qR(0);
  const auto uR = qR(1)/(rR + es);
  const auto pR = (gamma-1)*( qR(2) - half*rR*(uR*uR) );
  const auto HR = ( qR(2) + pR ) / rR;
  FR[0] = rR*uR;
  FR[1] = rR*uR*uR + pR;
  FR[2] = rR*uR*HR;

  // Compute Roe averaged variables
  const auto rLsqrt = std::sqrt(rL);
  const auto rRsqrt = std::sqrt(rR);
  const auto rSqrtSum = rLsqrt + rRsqrt;

  const auto uTilde = (rLsqrt * uL + rRsqrt * uR) / rSqrtSum;
  const auto HTilde = (rLsqrt * HL + rRsqrt * HR) / rSqrtSum;
  const auto aTilde = std::sqrt((gamma-1)*(HTilde - half * uTilde * uTilde));

  // Compute Eigenvalues
  auto lam1 = uTilde - aTilde;
  auto lam2 = uTilde;
  auto lam3 = uTilde + aTilde;

  // Entropy fix (San and Kara: Evaluation of Riemann flux solvers for WENO reconstruction
  // schemes: Kelvin–Helmholtz instability), Computers & Fluids, 117 (2015), 24-41. Eqn. 35
  const auto tol = 2 * es_ef * aTilde;
  if (std::abs(lam1) < tol) lam1 = lam1 * lam1 / (two * tol) + es_ef * aTilde;
  if (std::abs(lam2) < tol) lam2 = lam2 * lam2 / (two * tol) + es_ef * aTilde;
  if (std::abs(lam3) < tol) lam3 = lam3 * lam3 / (two * tol) + es_ef * aTilde;

  // Compute Eigenvectors
  K1[0] = one;
  K1[1] = uTilde - aTilde;
  K1[2] = HTilde - uTilde * aTilde;

  K2[0] = one;
  K2[1] = uTilde;
  K2[2] = half * uTilde * uTilde;

  K3[0] = one;
  K3[1] = uTilde + aTilde;
  K3[2] = HTilde + uTilde * aTilde;

  // Compute Eigenvector weights
  const auto dq1 = qR[0] - qL[0];
  const auto dq2 = qR[1] - qL[1];
  const auto dq3 = qR[2] - qL[2];

  const auto alpha2 = (gamma-1) / (aTilde * aTilde) * (dq1 * (HTilde - uTilde * uTilde) + uTilde * dq2 - dq3); 
  const auto alpha1 = (half / aTilde) * (dq1 * (uTilde + aTilde) - dq2 - aTilde * alpha2); 
  const auto alpha3 = dq1 - alpha1 - alpha2;

  // Compute Flux
  F(0) = half * (FL[0] + FR[0] - alpha1 * std::abs(lam1) * K1[0] - alpha2 * std::abs(lam2) * K2[0] - alpha3 * std::abs(lam3) * K3[0] );
  F(1) = half * (FL[1] + FR[1] - alpha1 * std::abs(lam1) * K1[1] - alpha2 * std::abs(lam2) * K2[1] - alpha3 * std::abs(lam3) * K3[1] );
  F(2) = half * (FL[2] + FR[2] - alpha1 * std::abs(lam1) * K1[2] - alpha2 * std::abs(lam2) * K2[2] - alpha3 * std::abs(lam3) * K3[2] );
}

template<class T, class normal_t, typename sc_t>
void eeRoeFluxFourDof(T & F,
               const T & qL,
               const T & qR,
               const normal_t & n,
               const sc_t gamma)
{
  constexpr auto one  = static_cast<sc_t>(1);
  constexpr auto two  = static_cast<sc_t>(2);
  constexpr auto half = one/two;

  constexpr sc_t es_ef{0.1}; // Entropy fix parameter
  constexpr sc_t es{1.e-30};
  sc_t FL[4];
  sc_t FR[4];
  sc_t K1[4];
  sc_t K2[4];
  sc_t K3[4];
  sc_t K4[4];

  // Compute Fluxes  
  // left
  const auto rL = qL(0);
  const auto uL = qL(1)/(rL + es);
  const auto vL = qL(2)/(rL + es);
  const auto unL = uL*n[0] + vL*n[1];
  const auto pL = (gamma-1)*( qL(3) - half*rL*(uL*uL+vL*vL) );
  const auto HL = ( qL(3) + pL ) / rL;
  FL[0] = rL*unL;
  FL[1] = rL*unL*uL + pL*n[0];
  FL[2] = rL*unL*vL + pL*n[1];
  FL[3] = rL*unL*HL;

  // right
  const auto rR = qR(0);
  const auto uR = qR(1)/(rR + es);
  const auto vR = qR(2)/(rR + es);
  const auto unR = uR*n[0] + vR*n[1];
  const auto pR = (gamma-1)*( qR(3) - half*rR*(uR*uR+vR*vR) );
  const auto HR = ( qR(3) + pR ) / rR;
  FR[0] = rR*unR;
  FR[1] = rR*unR*uR + pR*n[0];
  FR[2] = rR*unR*vR + pR*n[1];
  FR[3] = rR*unR*HR;

  // Compute Roe averaged variables
  const auto rLsqrt = std::sqrt(rL);
  const auto rRsqrt = std::sqrt(rR);
  const auto rSqrtSum = rLsqrt + rRsqrt;

  const auto uTilde  = (rLsqrt * uL + rRsqrt * uR) / rSqrtSum;
  const auto vTilde  = (rLsqrt * vL + rRsqrt * vR) / rSqrtSum;
  const auto unTilde = (rLsqrt * unL + rRsqrt * unR) / rSqrtSum;
  const auto HTilde  = (rLsqrt * HL + rRsqrt * HR) / rSqrtSum;
  const auto VTilde2 = uTilde * uTilde + vTilde * vTilde;
  const auto aTilde  = std::sqrt((gamma-1)*(HTilde - half * VTilde2));

  // Compute Eigenvalues
  auto lam1 = unTilde - aTilde;
  auto lam2 = unTilde;
  auto lam3 = unTilde;
  auto lam4 = unTilde + aTilde;

  // Entropy fix (San and Kara: Evaluation of Riemann flux solvers for WENO reconstruction
  // schemes: Kelvin–Helmholtz instability), Computers & Fluids, 117 (2015), 24-41. Eqn. 35
  const auto tol = two * es_ef * aTilde;
  if (std::abs(lam1) < tol) lam1 = lam1 * lam1 / (two * tol) + es_ef * aTilde;
  if (std::abs(lam2) < tol) lam2 = lam2 * lam2 / (two * tol) + es_ef * aTilde;
  if (std::abs(lam3) < tol) lam3 = lam3 * lam3 / (two * tol) + es_ef * aTilde;
  if (std::abs(lam4) < tol) lam4 = lam4 * lam4 / (two * tol) + es_ef * aTilde;

  // Compute Eigenvectors
  K1[0] = one;
  K1[1] = uTilde - aTilde*n[0];
  K1[2] = vTilde - aTilde*n[1];
  K1[3] = HTilde - unTilde * aTilde;

  K2[0] = one;
  K2[1] = uTilde;
  K2[2] = vTilde;
  K2[3] = half * VTilde2;

  K3[0] = 0.0;
  K3[1] = n[1];
  K3[2] = n[0];
  K3[3] = n[0]*vTilde + n[1]*uTilde;

  K4[0] = one;
  K4[1] = uTilde + aTilde*n[0];
  K4[2] = vTilde + aTilde*n[1];
  K4[3] = HTilde + unTilde * aTilde;

  // Compute Eigenvector weights
  const auto dq1 = qR[0] - qL[0];
  const auto dq2 = qR[1] - qL[1];
  const auto dq3 = qR[2] - qL[2];
  const auto dq4 = qR[3] - qL[3];

  const auto alpha3 = (dq3*n[0] - dq2*n[1] + dq1*n[1]*uTilde - dq1*n[0]*vTilde)/(n[0]*n[0] - n[1]*n[1]);
  const auto alpha2 = (gamma-1) / (aTilde * aTilde) * (dq1 * (HTilde - VTilde2) + uTilde * dq2 -dq4 +dq3*vTilde);

  const auto b1 = dq1-alpha2;
  const auto b2 = dq2-alpha2*uTilde-alpha3*n[1];
  const auto b3 = dq3 - alpha2*vTilde-alpha3*n[0];
  const auto alpha1 = (half / aTilde) * ((b1*(uTilde+aTilde*n[0])-b2)*n[0] + (b1*(vTilde+aTilde*n[1])-b3)*n[1]);
  const auto alpha4 = dq1 - alpha1 - alpha2;

  // Compute Flux
  F(0) = half * (FL[0] + FR[0] - alpha1 * std::abs(lam1) * K1[0] - alpha2 * std::abs(lam2) * K2[0] - alpha3 * std::abs(lam3) * K3[0] - alpha4 * std::abs(lam4) * K4[0]);
  F(1) = half * (FL[1] + FR[1] - alpha1 * std::abs(lam1) * K1[1] - alpha2 * std::abs(lam2) * K2[1] - alpha3 * std::abs(lam3) * K3[1] - alpha4 * std::abs(lam4) * K4[1]);
  F(2) = half * (FL[2] + FR[2] - alpha1 * std::abs(lam1) * K1[2] - alpha2 * std::abs(lam2) * K2[2] - alpha3 * std::abs(lam3) * K3[2] - alpha4 * std::abs(lam4) * K4[2]);
  F(3) = half * (FL[3] + FR[3] - alpha1 * std::abs(lam1) * K1[3] - alpha2 * std::abs(lam2) * K2[3] - alpha3 * std::abs(lam3) * K3[3] - alpha4 * std::abs(lam4) * K4[3]);
}

since code has changed quite a bit, putting these here and we also need at some point to add Jacobians