Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Roe flux #28

Open
pjb236 opened this issue Jul 1, 2021 · 3 comments
Open

Roe flux #28

pjb236 opened this issue Jul 1, 2021 · 3 comments
Assignees

Comments

@pjb236
Copy link
Contributor

pjb236 commented Jul 1, 2021

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 pjb236 self-assigned this Jul 1, 2021
pjb236 added a commit that referenced this issue Jul 26, 2021
pjb236 added a commit that referenced this issue Jul 29, 2021
pjb236 added a commit that referenced this issue Jul 29, 2021
… entropy fix is activated for one eigenvalue. #28
pjb236 added a commit that referenced this issue Jul 29, 2021
…sing Rusanov flux to be overwritten by Roe flux. #28
@pjb236
Copy link
Contributor Author

pjb236 commented Jul 29, 2021

@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 added a commit that referenced this issue Jul 30, 2021
pjb236 added a commit that referenced this issue Aug 5, 2021
Works for:
2D Normal shock weno5
2D Double Mach Reflection weno3

Does not work for:
2D Double Mach Reflection weno5
@pjb236
Copy link
Contributor Author

pjb236 commented Aug 5, 2021

@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
Copy link
Member

fnrizzi commented Dec 22, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants