Add four quadrant problem (for code comparison).

src/CMakeLists.txt

@@ -32,3 +32,4 @@ configure_file(test_blast.ini test_blast.ini COPYONLY)
 configure_file(test_blast_large.ini test_blast_large.ini COPYONLY)
 configure_file(test_implode.ini test_implode.ini COPYONLY)
 configure_file(test_implode_big.ini test_implode_big.ini COPYONLY)
+configure_file(test_four_quadrant.ini test_four_quadrant.ini COPYONLY)

src/HydroParams.cpp

@@ -90,6 +90,10 @@ HydroParams::setup(ConfigMap & configMap)
   {
     problemType = PROBLEM_BLAST;
   }
+  else if (!problemStr.compare("four_quadrant"))
+  {
+    problemType = PROBLEM_FOUR_QUADRANT;
+  }
   else
   {
     std::cout << "Problem is invalid\n";

src/HydroParams.h

@@ -90,7 +90,8 @@ enum ImplementationVersion
 enum ProblemType
 {
   PROBLEM_IMPLODE,
-  PROBLEM_BLAST
+  PROBLEM_BLAST,
+  PROBLEM_FOUR_QUADRANT
 };
 
 // variable names in the order as in component index

src/HydroRun.h

@@ -36,7 +36,7 @@ isBigEndian()
 /**
  * Main hydrodynamics data structure.
  */
-template<typename device_t>
+template <typename device_t>
 class HydroRun
 {
 
@@ -95,6 +95,8 @@ class HydroRun
   init_implode(DataArray_t Udata);
   void
   init_blast(DataArray_t Udata);
+  void
+  init_four_quadrant(DataArray_t Udata);
 
   // host routines (save data to file, device data are copied into host
   // inside this routine)
@@ -123,7 +125,7 @@ class HydroRun
 /**
  *
  */
-template<typename device_t>
+template <typename device_t>
 HydroRun<device_t>::HydroRun(HydroParams & params, ConfigMap & configMap)
   : params(params)
   , configMap(configMap)
@@ -174,6 +176,10 @@ HydroRun<device_t>::HydroRun(HydroParams & params, ConfigMap & configMap)
   {
     init_blast(U);
   }
+  else if (params.problemType == PROBLEM_FOUR_QUADRANT)
+  {
+    init_four_quadrant(U);
+  }
   else
   {
     std::cout << "Problem : " << params.problemType
@@ -381,7 +387,29 @@ HydroRun<device_t>::init_blast(DataArray_t Udata)
 
 // =======================================================
 // =======================================================
-template<typename device_t>
+/**
+ * Hydrodynamical four quadrant Test.
+ *
+ * In the 2D case, there are 19 different possible configurations (see
+ * article by Lax and Liu, "Solution of two-dimensional riemann
+ * problems of gas dynamics by positive schemes",SIAM journal on
+ * scientific computing, 1998, vol. 19, no2, pp. 319-340).
+ *
+ * Here only problem #3 is implemented. See https://github.com/pkestene/euler_kokkos for complete
+ * list of all 19 initial conditions.
+ */
+template <typename device_t>
+void
+HydroRun<device_t>::init_four_quadrant(DataArray_t Udata)
+{
+
+  InitFourQuadrantFunctor<device_t>::apply(params, Udata);
+
+} // HydroRun<device_t>::init_four_quadrant
+
+// =======================================================
+// =======================================================
+template <typename device_t>
 void
 HydroRun<device_t>::saveData(DataArray_t Udata, int iStep, std::string name)
 {

src/HydroRunFunctors.h

@@ -1211,7 +1211,128 @@ class InitBlastFunctor : public HydroBaseFunctor
 /*************************************************/
 /*************************************************/
 /*************************************************/
-template<typename device_t, FaceIdType faceId>
+template <typename device_t>
+class InitFourQuadrantFunctor : public HydroBaseFunctor
+{
+
+public:
+  using DataArray_t = DataArray<device_t>;
+  using exec_space = typename device_t::execution_space;
+
+  InitFourQuadrantFunctor(HydroParams params, DataArray_t Udata)
+    : HydroBaseFunctor(params)
+    , Udata(Udata){};
+
+  // static method which does it all: create and execute functor
+  static void
+  apply(HydroParams params, DataArray_t Udata)
+  {
+    InitFourQuadrantFunctor functor(params, Udata);
+    Kokkos::parallel_for(
+      "InitFourQuadrant",
+      Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>>({ 0, 0 }, { params.isize, params.jsize }),
+      functor);
+  }
+
+  KOKKOS_FUNCTION void
+  primToCons(Kokkos::Array<real_t, 4> & U, real_t gamma0) const
+  {
+
+    real_t rho = U[ID];
+    real_t p = U[IP];
+    real_t u = U[IU];
+    real_t v = U[IV];
+
+    U[IU] *= rho; // rho*u
+    U[IV] *= rho; // rho*v
+    U[IP] = p / (gamma0 - 1.0) + rho * (u * u + v * v) * 0.5;
+
+  } // primToCons
+
+
+  KOKKOS_INLINE_FUNCTION
+  void
+  operator()(const int & i, const int & j) const
+  {
+
+    const int ghostWidth = params.ghostWidth;
+
+    const real_t xmin = params.xmin;
+    const real_t ymin = params.ymin;
+    const real_t dx = params.dx;
+    const real_t dy = params.dy;
+
+    // four quadrant problem parameters
+    const real_t xt = 0.8;
+    const real_t yt = 0.8;
+
+    const real_t gamma0 = params.settings.gamma0;
+
+    // clang-format off
+    // primitive variables
+    Kokkos::Array<real_t, 4> U0{ 1.5   , 1.5  , 0.0  , 0.0   };
+    Kokkos::Array<real_t, 4> U1{ 0.5323, 0.3  , 1.206, 0.0   };
+    Kokkos::Array<real_t, 4> U2{ 0.138 , 0.029, 1.206, 1.206 };
+    Kokkos::Array<real_t, 4> U3{ 0.5323, 0.3  , 0.0  , 1.206 };
+    // clang-format on
+
+    primToCons(U0, gamma0);
+    primToCons(U1, gamma0);
+    primToCons(U2, gamma0);
+    primToCons(U3, gamma0);
+
+    real_t x = xmin + dx / 2 + (i - ghostWidth) * dx;
+    real_t y = ymin + dy / 2 + (j - ghostWidth) * dy;
+
+    if (x < xt)
+    {
+      if (y < yt)
+      {
+        // region 2
+        Udata(i, j, ID) = U2[ID];
+        Udata(i, j, IP) = U2[IP];
+        Udata(i, j, IU) = U2[IU];
+        Udata(i, j, IV) = U2[IV];
+      }
+      else
+      {
+        // region 1
+        Udata(i, j, ID) = U1[ID];
+        Udata(i, j, IP) = U1[IP];
+        Udata(i, j, IU) = U1[IU];
+        Udata(i, j, IV) = U1[IV];
+      }
+    }
+    else
+    {
+      if (y < yt)
+      {
+        // region 3
+        Udata(i, j, ID) = U3[ID];
+        Udata(i, j, IP) = U3[IP];
+        Udata(i, j, IU) = U3[IU];
+        Udata(i, j, IV) = U3[IV];
+      }
+      else
+      {
+        // region 0
+        Udata(i, j, ID) = U0[ID];
+        Udata(i, j, IP) = U0[IP];
+        Udata(i, j, IU) = U0[IU];
+        Udata(i, j, IV) = U0[IV];
+      }
+    }
+
+  } // end operator ()
+
+  DataArray_t Udata;
+
+}; // InitFourQuadrantFunctor
+
+/*************************************************/
+/*************************************************/
+/*************************************************/
+template <typename device_t, FaceIdType faceId>
 class MakeBoundariesFunctor : public HydroBaseFunctor
 {
 

src/test_four_quadrant.ini

@@ -0,0 +1,33 @@
+[run]
+tEnd=0.6
+nStepmax=1200
+nOutput=100
+
+[mesh]
+nx=256
+ny=256
+boundary_type_xmin=2
+boundary_type_xmax=2
+boundary_type_ymin=2
+boundary_type_ymax=2
+
+xmin=0.0
+xmax=1.0
+
+ymin=0.0
+ymax=1.0
+
+[hydro]
+gamma0=1.666
+cfl=0.8
+niter_riemann=10
+iorder=2
+slope_type=2
+problem=four_quadrant
+riemann=hllc
+
+[output]
+outputPrefix=test_four_quadrant
+
+[other]
+implementationVersion=0