WIP: added 2D inviscid Roe flux. #28

Works for:
2D Normal shock weno5
2D Double Mach Reflection weno3

Does not work for:
2D Double Mach Reflection weno5

include/pressiodemoapps/euler2d.hpp

@@ -43,12 +43,12 @@ T createEe2dImpl(const mesh_t & meshObj,
       ("Stencil size in the mesh object not compatible with desired inviscid flux reconstruction.");
   }
 
-  if (probEnum == ::pressiodemoapps::Euler2d::DoubleMachReflection and
-      recEnum == ::pressiodemoapps::InviscidFluxReconstruction::Weno5)
-  {
-    throw std::runtime_error
-      ("Double mach reflection does NOT currently suppot Weno5.");
-  }
+  //if (probEnum == ::pressiodemoapps::Euler2d::DoubleMachReflection and
+  //    recEnum == ::pressiodemoapps::InviscidFluxReconstruction::Weno5)
+  //{
+  //  throw std::runtime_error
+  //    ("Double mach reflection does NOT currently suppot Weno5.");
+  //}
 
 
   if (probEnum == ::pressiodemoapps::Euler2d::PeriodicSmooth)

include/pressiodemoapps/impl_apps/euler2d/euler2d_app_impl.hpp

@@ -248,6 +248,10 @@ class Euler2dAppT
 	  eeRusanovFluxFourDof(FL, uMinusHalfNeg, uMinusHalfPos, normalX_, m_gamma);
 	  eeRusanovFluxFourDof(FR, uPlusHalfNeg,  uPlusHalfPos,  normalX_, m_gamma);
 	  break;
+	case ::pressiodemoapps::InviscidFluxScheme::Roe:
+	  eeRoeFluxFourDof(FL, uMinusHalfNeg, uMinusHalfPos, normalX_, m_gamma);
+	  eeRoeFluxFourDof(FR, uPlusHalfNeg,  uPlusHalfPos,  normalX_, m_gamma);
+	  break;
 	}
 
       // Y
@@ -259,6 +263,10 @@ class Euler2dAppT
 	  eeRusanovFluxFourDof(FD, uMinusHalfNeg, uMinusHalfPos, normalY_, m_gamma);
 	  eeRusanovFluxFourDof(FU, uPlusHalfNeg,  uPlusHalfPos,  normalY_, m_gamma);
 	  break;
+	case ::pressiodemoapps::InviscidFluxScheme::Roe:
+	  eeRoeFluxFourDof(FD, uMinusHalfNeg, uMinusHalfPos, normalY_, m_gamma);
+	  eeRoeFluxFourDof(FU, uPlusHalfNeg,  uPlusHalfPos,  normalY_, m_gamma);
+	  break;
 	}
 
       const auto vIndex = smPt*numDofPerCell;
@@ -319,6 +327,10 @@ class Euler2dAppT
 	  eeRusanovFluxFourDof(FL, uMinusHalfNeg, uMinusHalfPos, normalX_, m_gamma);
 	  eeRusanovFluxFourDof(FR, uPlusHalfNeg,  uPlusHalfPos,  normalX_, m_gamma);
 	  break;
+	case ::pressiodemoapps::InviscidFluxScheme::Roe:
+	  eeRoeFluxFourDof(FL, uMinusHalfNeg, uMinusHalfPos, normalX_, m_gamma);
+	  eeRoeFluxFourDof(FR, uPlusHalfNeg,  uPlusHalfPos,  normalX_, m_gamma);
+	  break;
 	}
 
       // Y
@@ -329,7 +341,11 @@ class Euler2dAppT
 	  eeRusanovFluxFourDof(FD, uMinusHalfNeg, uMinusHalfPos, normalY_, m_gamma);
 	  eeRusanovFluxFourDof(FU, uPlusHalfNeg,  uPlusHalfPos,  normalY_, m_gamma);
 	  break;
-	}
+	case ::pressiodemoapps::InviscidFluxScheme::Roe:
+	  eeRoeFluxFourDof(FD, uMinusHalfNeg, uMinusHalfPos, normalY_, m_gamma);
+	  eeRoeFluxFourDof(FU, uPlusHalfNeg,  uPlusHalfPos,  normalY_, m_gamma);
+	  break;
+    }
 
       const auto vIndex = smPt*numDofPerCell;
       vEval(vIndex);

include/pressiodemoapps/impl_apps/eulerCommon/fluxes.hpp

@@ -134,7 +134,7 @@ void eeRusanovFluxThreeDof(T & F,
   F(2) = half*( FL[2] + FR[2] + smax*(qL(2) - qR(2)) );
 }
 
-template<class T, typename sc_t>
+template<class T, class normal_t, typename sc_t>
 void eeRoeFluxFourDof(T & F,
                const T & qL,
                const T & qR,
@@ -189,13 +189,14 @@ void eeRoeFluxFourDof(T & F,
   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 * VTilde));
+  const auto aTilde  = std::sqrt((gamma-1)*(HTilde - half * VTilde2));
 
 
   // Compute Eigenvalues
   auto lam1 = unTilde - aTilde;
   auto lam2 = unTilde;
-  auto lam3 = unTilde + aTilde;
+  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
@@ -232,8 +233,8 @@ void eeRoeFluxFourDof(T & F,
   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]*u - 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 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];

tests_cpp/2d_double_mach_reflection/CMakeLists.txt

@@ -1,3 +1,4 @@
 
 add_subdirectory(firstorder)
-add_subdirectory(weno3)
+add_subdirectory(weno3)
+add_subdirectory(weno5)

tests_cpp/2d_double_mach_reflection/main.cc

@@ -8,14 +8,19 @@ int main(int argc, char *argv[])
 
   namespace pda = pressiodemoapps;
   const auto meshObj = pda::loadCellCenterUniformMeshEigen(".");
-#if defined USE_WENO3
+#ifdef USE_WENO5
+  const auto order   = pda::InviscidFluxReconstruction::Weno5;
+  const auto iFlux   = pda::InviscidFluxScheme::Roe;
+#elif defined USE_WENO3
   const auto order   = pda::InviscidFluxReconstruction::Weno3;
+  const auto iFlux   = pda::InviscidFluxScheme::Rusanov;
 #else
   const auto order   = pda::InviscidFluxReconstruction::FirstOrder;
+  const auto iFlux   = pda::InviscidFluxScheme::Rusanov;
 #endif
 
   const auto probId  = pda::Euler2d::DoubleMachReflection;
-  auto appObj      = pda::createProblemEigen(meshObj, probId, order);
+  auto appObj      = pda::createProblemEigen(meshObj, probId, order, iFlux, 1);
   using app_t = decltype(appObj);
   using app_state_t = typename app_t::state_type;
   using ode_state_t = pressio::containers::Vector<app_state_t>;
@@ -28,7 +33,7 @@ int main(int argc, char *argv[])
   const auto dt = 0.0005;
   const auto Nsteps = 0.25/dt;
   auto time = 0.0;
-  FomObserver<ode_state_t> Obs("doubleMach2d_solution.bin", 50);
+  FomObserver<ode_state_t> Obs("doubleMach2d_solution.bin", 1);
   app_state_t state1(stateSize);
   app_state_t state2(stateSize);
   app_state_t veloTmp(stateSize);