Add init conditions for shocked bubble test case.

src/CMakeLists.txt

@@ -36,5 +36,6 @@ 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)
 configure_file(test_discontinuity.ini test_discontinuity.ini COPYONLY)
+configure_file(test_shocked_bubble.ini test_shocked_bubble.ini COPYONLY)
 
 add_subdirectory(sedov_blast)

src/HydroParams.cpp

@@ -12,6 +12,26 @@ namespace euler2d
 
 const char * varNames[4] = { "rho", "E", "mx", "my" };
 
+ShockedBubbleParams
+readShockedBubbleParams(const ConfigMap & configMap)
+{
+  ShockedBubbleParams par;
+
+  par.bubble_radius = configMap.getFloat("shocked_bubble", "bubble_radius", 0.025);
+  par.bubble_center_x = configMap.getFloat("shocked_bubble", "bubble_center_x", 0.225);
+  par.bubble_center_y = configMap.getFloat("shocked_bubble", "bubble_center_y", 0.0445);
+  par.bubble_density = configMap.getFloat("shocked_bubble", "bubble_density", 3.863);
+  par.bubble_pressure = configMap.getFloat("shocked_bubble", "bubble_pressure", 1.0132e5);
+  par.preshock_density = configMap.getFloat("shocked_bubble", "preshock_density", 1.225);
+  par.preshock_pressure = configMap.getFloat("shocked_bubble", "preshock_pressure", 1.0132e5);
+  par.postshock_density = configMap.getFloat("shocked_bubble", "postshock_density", 1.686);
+  par.postshock_pressure = configMap.getFloat("shocked_bubble", "postshock_pressure", 1.59e5);
+  par.postshock_velocity = configMap.getFloat("shocked_bubble", "postshock_velocity", 113.5);
+  par.shock_loc = configMap.getFloat("shocked_bubble", "shock_loc", 0.170);
+
+  return par;
+}
+
 /*
  * Hydro Parameters (read parameter file)
  */
@@ -97,6 +117,11 @@ HydroParams::setup(ConfigMap & configMap)
   {
     problemType = PROBLEM_DISCONTINUITY;
   }
+  else if (!problemStr.compare("shocked_bubble"))
+  {
+    problemType = PROBLEM_SHOCKED_BUBBLE;
+    shock_par = readShockedBubbleParams(configMap);
+  }
   else
   {
     std::cout << "Problem is invalid\n";

src/HydroParams.h

@@ -24,6 +24,7 @@ enum ComponentIndex
 {
   ID = 0, /*!< ID Density field index */
   IP = 1, /*!< IP Pressure/Energy field index */
+  IE = 1,
   IU = 2, /*!< X velocity / momentum index */
   IV = 3  /*!< Y velocity / momentum index */
 };
@@ -92,7 +93,8 @@ enum ProblemType
   PROBLEM_IMPLODE,
   PROBLEM_BLAST,
   PROBLEM_FOUR_QUADRANT,
-  PROBLEM_DISCONTINUITY
+  PROBLEM_DISCONTINUITY,
+  PROBLEM_SHOCKED_BUBBLE
 };
 
 // variable names in the order as in component index
@@ -125,6 +127,24 @@ struct HydroSettings
 
 }; // struct HydroSettings
 
+struct ShockedBubbleParams
+{
+  real_t bubble_radius;
+  real_t bubble_center_x;
+  real_t bubble_center_y;
+  real_t bubble_density;
+  real_t bubble_pressure;
+  real_t preshock_density;
+  real_t preshock_pressure;
+  real_t postshock_density;
+  real_t postshock_pressure;
+  real_t postshock_velocity;
+  real_t shock_loc;
+};
+
+ShockedBubbleParams
+readShockedBubbleParams(const ConfigMap & configMap);
+
 /*
  * Hydro Parameters (declaration)
  */
@@ -184,6 +204,9 @@ struct HydroParams
   real_t blast_total_energy_inside;
   int    blast_nbins;
 
+  // shocked bubble parameters
+  ShockedBubbleParams shock_par;
+
   // other parameters
   int implementationVersion = 0; /*!< triggers which implementation to use (currently 3 versions)*/
 

src/HydroRun.h

@@ -100,6 +100,8 @@ class HydroRun
   init_four_quadrant(DataArray_t Udata);
   void
   init_discontinuity(DataArray_t Udata);
+  void
+  init_shocked_bubble(DataArray_t Udata);
 
   // host routines (save data to file, device data are copied into host
   // inside this routine)
@@ -187,6 +189,10 @@ HydroRun<device_t>::HydroRun(HydroParams & params, ConfigMap & configMap)
   {
     init_discontinuity(U);
   }
+  else if (params.problemType == PROBLEM_SHOCKED_BUBBLE)
+  {
+    init_shocked_bubble(U);
+  }
   else
   {
     std::cout << "Problem : " << params.problemType
@@ -412,6 +418,21 @@ HydroRun<device_t>::godunov_unsplit_impl(DataArray_t data_in,
 
   } // HydroRun<device_t>::init_discontinuity
 
+  // =======================================================
+  // =======================================================
+  /**
+   * Shocked bubble.
+   *
+   */
+  template <typename device_t>
+  void
+  HydroRun<device_t>::init_shocked_bubble(DataArray_t Udata)
+  {
+
+    InitShockedBubbleFunctor<device_t>::apply(params, Udata);
+
+  } // HydroRun<device_t>::init_shocked_bubble
+
   // =======================================================
   // =======================================================
   /**
@@ -426,7 +447,8 @@ HydroRun<device_t>::godunov_unsplit_impl(DataArray_t data_in,
    * list of all 19 initial conditions.
    */
   template <typename device_t>
-  void HydroRun<device_t>::init_four_quadrant(DataArray_t Udata)
+  void
+  HydroRun<device_t>::init_four_quadrant(DataArray_t Udata)
   {
 
     InitFourQuadrantFunctor<device_t>::apply(params, Udata);

src/HydroRunFunctors.h

@@ -1709,6 +1709,98 @@ class InitDiscontinuityFunctor : public HydroBaseFunctor
 
 }; // InitDiscontinuityFunctor
 
+/*************************************************/
+/*************************************************/
+/*************************************************/
+template <typename device_t>
+class InitShockedBubbleFunctor : public HydroBaseFunctor
+{
+
+public:
+  using DataArray_t = DataArray<device_t>;
+  using exec_space = typename device_t::execution_space;
+
+  InitShockedBubbleFunctor(HydroParams params, DataArray_t Udata)
+    : HydroBaseFunctor(params)
+    , Udata(Udata)
+    , shock_par(params.shock_par){};
+
+  // static method which does it all: create and execute functor
+  static void
+  apply(HydroParams params, DataArray_t Udata)
+  {
+    InitShockedBubbleFunctor functor(params, Udata);
+    Kokkos::parallel_for(
+      "InitShockedBubble",
+      Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>>({ 0, 0 }, { params.isize, params.jsize }),
+      functor);
+  }
+
+  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;
+
+    real_t x = xmin + dx / 2 + (i - ghostWidth) * dx;
+    real_t y = ymin + dy / 2 + (j - ghostWidth) * dy;
+
+    const real_t gamma0 = params.settings.gamma0;
+
+    // post shock region (coming from the left)
+    if (x < shock_par.shock_loc)
+    {
+      Udata(i, j, ID) = shock_par.postshock_density;
+      Udata(i, j, IU) = shock_par.postshock_density * shock_par.postshock_velocity;
+      Udata(i, j, IV) = 0.0;
+      const real_t rho_eint = shock_par.postshock_pressure / (gamma0 - 1);
+      Udata(i, j, IE) =
+        rho_eint + 0.5 * (Udata(i, j, IU) * Udata(i, j, IU) + Udata(i, j, IV) * Udata(i, j, IV)) /
+                     Udata(i, j, ID);
+    }
+    else
+    {
+      // inside or outside bubble ?
+      const auto radius = sqrt((x - shock_par.bubble_center_x) * (x - shock_par.bubble_center_x) +
+                               (y - shock_par.bubble_center_y) * (y - shock_par.bubble_center_y));
+
+      if (radius < shock_par.bubble_radius)
+      {
+        // bubble at rest
+        Udata(i, j, ID) = shock_par.bubble_density;
+        Udata(i, j, IU) = 0.0;
+        Udata(i, j, IV) = 0.0;
+        const real_t rho_eint = shock_par.bubble_pressure / (gamma0 - 1);
+        Udata(i, j, IE) =
+          rho_eint + 0.5 * (Udata(i, j, IU) * Udata(i, j, IU) + Udata(i, j, IV) * Udata(i, j, IV)) /
+                       Udata(i, j, ID);
+      }
+      else
+      {
+        // air at rest
+        Udata(i, j, ID) = shock_par.preshock_density;
+        Udata(i, j, IU) = 0.0;
+        Udata(i, j, IV) = 0.0;
+        const real_t rho_eint = shock_par.preshock_pressure / (gamma0 - 1);
+        Udata(i, j, IE) =
+          rho_eint + 0.5 * (Udata(i, j, IU) * Udata(i, j, IU) + Udata(i, j, IV) * Udata(i, j, IV)) /
+                       Udata(i, j, ID);
+      }
+    }
+
+  } // end operator ()
+
+  DataArray_t         Udata;
+  ShockedBubbleParams shock_par;
+
+}; // InitShockedBubbleFunctor
+
 /*************************************************/
 /*************************************************/
 /*************************************************/

src/test_shocked_bubble.ini

@@ -0,0 +1,39 @@
+[run]
+tEnd=10.0
+nStepmax=500
+nOutput=10
+
+[mesh]
+nx=445
+ny=89
+
+xmin=0.0
+xmax=0.445
+
+ymin=0.0
+ymax=0.089
+
+boundary_type_xmin=2
+boundary_type_xmax=2
+
+boundary_type_ymin=1
+boundary_type_ymax=1
+
+[hydro]
+gamma0=1.2
+cfl=0.5
+niter_riemann=10
+slope_type=2
+problem=shocked_bubble
+riemann=hllc
+
+#[shocked_bubble]
+#bubble_pressure=1.0132
+#preshock_pressure=1.0132
+#postshock_pressure=1.59
+
+[output]
+outputPrefix=test_shocked_bubble
+
+[other]
+implementationVersion=0