Add Sedov blast input parameter file.

Also provide python to plot radial density profile and compare with analytical solution.

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.22)
 project(
   euler2d
   DESCRIPTION "euler2d_kokkos is mini application using C++/Kokkos library, used for teaching purpose"
-  LANGUAGES CXX C)
+  LANGUAGES CXX C Fortran)
 
 #
 # default local cmake macro repository

src/CMakeLists.txt

@@ -17,6 +17,8 @@ add_executable(euler2d
   SimpleTimer.h
   SimpleTimer.cpp
   CudaTimer.h
+  cnpy/cnpy.h
+  cnpy/cnpy.cpp
   main.cpp)
 
 target_include_directories(euler2d
@@ -34,3 +36,5 @@ 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)
+
+add_subdirectory(sedov_blast)

src/ComputeRadialProfileFunctor.h

@@ -0,0 +1,169 @@
+/**
+ * \file ComputeRadialProfile.h
+ */
+#ifndef EULER2D_COMPUTE_RADIAL_PROFILE_FUNCTOR_H_
+#define EULER2D_COMPUTE_RADIAL_PROFILE_FUNCTOR_H_
+
+#include "kokkos_shared.h"
+#include "HydroBaseFunctor.h"
+#include "HydroParams.h"
+#include "real_type.h"
+
+// other utilities
+#include "cnpy/cnpy_io.h"
+
+namespace euler2d
+{
+
+/*************************************************/
+/*************************************************/
+/*************************************************/
+/**
+ * Compute radial profile by averaging all direction.
+ *
+ * This functor is mainly aimed at checking numerical and analytical solution of the radial Sedov
+ * blast.
+ *
+ * \tparam device_t is the Kokkos device use for computation (CPU, GPU, ...)
+ */
+template <typename device_t>
+class ComputeRadialProfileFunctor : HydroBaseFunctor
+{
+
+public:
+  using DataArray_t = DataArray<device_t>;
+
+  //! our kokkos execution space
+  using exec_space = typename device_t::execution_space;
+
+private:
+  //! heavy data - conservative variables
+  DataArray_t m_Udata;
+
+  //! number of bins in radial direction
+  const int m_nbins;
+
+  //! max radial distance
+  const real_t m_max_radial_distance;
+
+  //! center of the box
+  Kokkos::Array<real_t, 2> m_box_center;
+
+  //! radial distance histogram
+  Kokkos::View<int *, device_t, Kokkos::MemoryTraits<Kokkos::Atomic>> m_radial_distance_histo;
+
+  //! averaged density radial profile
+  Kokkos::View<real_t *, device_t, Kokkos::MemoryTraits<Kokkos::Atomic>> m_density_profile;
+
+public:
+  ComputeRadialProfileFunctor(HydroParams              params,
+                              DataArray_t              Udata,
+                              real_t                   max_radial_distance,
+                              Kokkos::Array<real_t, 2> box_center)
+    : HydroBaseFunctor(params)
+    , m_Udata(Udata)
+    , m_nbins(params.blast_nbins)
+    , m_max_radial_distance(max_radial_distance)
+    , m_box_center(box_center)
+    , m_radial_distance_histo("radial_distance_histo", m_nbins)
+    , m_density_profile("density_profile", m_nbins)
+  {
+    Kokkos::deep_copy(m_radial_distance_histo, 0);
+    Kokkos::deep_copy(m_density_profile, ZERO_F);
+  };
+
+  // ====================================================================
+  // ====================================================================
+  //! static method which does it all: create and execute functor using range policy
+  //!
+  static void
+  apply(HydroParams params, DataArray_t Udata)
+  {
+
+    const real_t xmin = params.xmin;
+    const real_t ymin = params.ymin;
+    const real_t xmax = params.xmax;
+    const real_t ymax = params.ymax;
+    const real_t dx = params.dx;
+    const real_t dy = params.dy;
+
+    // get center of the box coordinates
+    Kokkos::Array<real_t, 2> box_center;
+
+    box_center[IX] = (xmin + xmax) / 2;
+    box_center[IY] = (ymin + ymax) / 2;
+
+    // compute maximum radial distance
+    const auto Dx = (xmax - xmin) / 2;
+    const auto Dy = (ymax - ymin) / 2;
+    const auto max_radial_distance = sqrt(Dx * Dx + Dy * Dy);
+
+    ComputeRadialProfileFunctor functor(params, Udata, max_radial_distance, box_center);
+
+    Kokkos::parallel_for(
+      "euler2d::ComputeRadialProfileFunctor",
+      Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>>({ 0, 0 }, { params.isize, params.jsize }),
+      functor);
+
+    // once the radial profile is computed, gather all the MPI pieces
+    const auto nbins = params.blast_nbins;
+    auto       total_radial_distance_histo =
+      Kokkos::View<int *, device_t>("total_radial_distance_histo", nbins);
+
+    const auto total_density_profile =
+      Kokkos::View<real_t *, device_t>("total_density_profile", nbins);
+
+    // then we compute the profile and output in numpy format
+    auto rad_dist =
+      Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, functor.m_radial_distance_histo);
+    auto dens_prof =
+      Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, functor.m_density_profile);
+
+    auto       distances = Kokkos::View<real_t *, Kokkos::HostSpace>("distances", nbins);
+    const auto dr = max_radial_distance / nbins;
+
+    for (int i = 0; i < nbins; ++i)
+    {
+      distances(i) = (i + 0.5) * dr;
+      dens_prof(i) /= rad_dist(i);
+    }
+
+    // now we can output results
+    save_cnpy(distances, "sedov_blast_radial_distances");
+    save_cnpy(dens_prof, "sedov_blast_density_profile");
+
+  } // apply
+
+  // ====================================================================
+  // ====================================================================
+  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;
+    const real_t x = xmin + dx / 2 + (i - ghostWidth) * dx;
+    const real_t y = ymin + dy / 2 + (j - ghostWidth) * dy;
+
+    const auto distance = sqrt((x - m_box_center[IX]) * (x - m_box_center[IX]) +
+                               (y - m_box_center[IY]) * (y - m_box_center[IY]));
+
+    // compute bin number from distance to center of the box
+    const auto bin = (int)(distance / m_max_radial_distance * m_nbins);
+
+    // update density profile
+    m_radial_distance_histo(bin) += 1;
+
+    m_density_profile(bin) += m_Udata(i, j, ID);
+
+  } // operator()
+
+}; // ComputeRadialProfileFunctor
+
+} // namespace euler2d
+
+#endif // EULER2D_COMPUTE_RADIAL_PROFILE_FUNCTOR_H_

src/HydroParams.cpp

@@ -112,6 +112,8 @@ HydroParams::setup(ConfigMap & configMap)
   blast_density_out = configMap.getFloat("blast", "density_out", 1.2);
   blast_pressure_in = configMap.getFloat("blast", "pressure_in", 10.0);
   blast_pressure_out = configMap.getFloat("blast", "pressure_out", 0.1);
+  blast_total_energy_inside = configMap.getFloat("blast", "total_energy_inside", 0.0);
+  blast_nbins = configMap.getInteger("blast", "nbins", 100);
 
   implementationVersion = configMap.getFloat("OTHER", "implementationVersion", 0);
   if (implementationVersion != 0 and implementationVersion != 1 and implementationVersion != 2)

src/HydroParams.h

@@ -183,6 +183,8 @@ struct HydroParams
   real_t blast_density_out;
   real_t blast_pressure_in;
   real_t blast_pressure_out;
+  real_t blast_total_energy_inside;
+  int    blast_nbins;
 
   // other parameters
   int implementationVersion = 0; /*!< triggers which implementation to use (currently 3 versions)*/
@@ -224,6 +226,7 @@ struct HydroParams
     , blast_density_out(1.2)
     , blast_pressure_in(10.0)
     , blast_pressure_out(0.1)
+    , blast_total_energy_inside(0)
     , implementationVersion(0)
   {}
 

src/HydroRun.h

@@ -118,7 +118,7 @@ class HydroRun
    * Paraview/Visit to read these h5 files as a time series.
    */
   void
-  write_xdmf_time_series();
+  write_xdmf_time_series(int last_time_step);
 #endif
 }; // class HydroRun
 
@@ -696,7 +696,8 @@ HydroRun<device_t>::godunov_unsplit_impl(DataArray_t data_in,
    * point to data file : xdmf2d.h5
    */
   template <typename device_t>
-  void HydroRun<device_t>::write_xdmf_time_series()
+  void
+  HydroRun<device_t>::write_xdmf_time_series(int last_time_step)
   {
     FILE *      xdmf = 0;
     const int & nx = params.nx;
@@ -722,8 +723,15 @@ HydroRun<device_t>::godunov_unsplit_impl(DataArray_t data_in,
             "    <Grid Name=\"TimeSeries\" GridType=\"Collection\" CollectionType=\"Temporal\">\n");
 
     // for each time step write a <grid> </grid> item
-    for (int iStep = 0; iStep <= params.nStepmax; iStep += params.nOutput)
+    const int nbOutput = (params.nStepmax + params.nOutput - 1) / params.nOutput + 1;
+    for (int iOut = 0; iOut < nbOutput; ++iOut)
     {
+      int iStep = iOut * params.nOutput;
+
+      if (iOut == (nbOutput - 1))
+      {
+        iStep = last_time_step;
+      }
 
       std::ostringstream outNum;
       outNum.width(7);
@@ -786,7 +794,7 @@ HydroRun<device_t>::godunov_unsplit_impl(DataArray_t data_in,
     fprintf(xdmf, "</Xdmf>\n");
     fclose(xdmf);
 
-  }    // HydroRun<device_t>::write_xdmf_xml
+  } // HydroRun<device_t>::write_xdmf_xml
 #endif // USE_HDF5
 
 } // namespace euler2d

src/HydroRunFunctors.h

@@ -1384,32 +1384,63 @@ class InitImplodeFunctor : public HydroBaseFunctor
 /*************************************************/
 /*************************************************/
 /*************************************************/
-template<typename device_t>
+/**
+ * Initialize blast test case.
+ *
+ * Two types of initialization:
+ *
+ * - if parameter total_energy_inside is positive, we initialize in total energy
+ * - if parameter total_energy_inside is negative or null, we initialize using pressure in /
+ * pressure out
+ *
+ * If you want to do the well know sedov blast test, you need to initialize using total energy, not
+ * pressure.
+ */
+template <typename device_t>
 class InitBlastFunctor : public HydroBaseFunctor
 {
 
 public:
+  struct TagRegularInit
+  {};
+  struct TagCorrectTotalEnergy
+  {};
+
   using DataArray_t = DataArray<device_t>;
   using exec_space = typename device_t::execution_space;
 
   InitBlastFunctor(HydroParams params, DataArray_t Udata)
     : HydroBaseFunctor(params)
-    , Udata(Udata){};
+    , Udata(Udata) {};
 
   // static method which does it all: create and execute functor
   static void
   apply(HydroParams params, DataArray_t Udata)
   {
+    real_t              volume_inside = ZERO_F;
+    Kokkos::Sum<real_t> reducer(volume_inside);
+
     InitBlastFunctor functor(params, Udata);
-    Kokkos::parallel_for(
-      "InitBlast",
-      Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>>({ 0, 0 }, { params.isize, params.jsize }),
-      functor);
+    Kokkos::parallel_reduce("InitBlast - regular",
+                            Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>, TagRegularInit>(
+                              { 0, 0 }, { params.isize, params.jsize }),
+                            functor,
+                            reducer);
+
+    if (params.blast_total_energy_inside > 0)
+    {
+      functor.m_volume_inside = volume_inside;
+      Kokkos::parallel_for(
+        "InitBlast - correct total energy",
+        Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>, TagCorrectTotalEnergy>(
+          { 0, 0 }, { params.isize, params.jsize }),
+        functor);
+    }
   }
 
   KOKKOS_INLINE_FUNCTION
   void
-  operator()(const int & i, const int & j) const
+  operator()(TagRegularInit const &, const int & i, const int & j, real_t & volume) const
   {
 
     const int ghostWidth = params.ghostWidth;
@@ -1443,6 +1474,8 @@ class InitBlastFunctor : public HydroBaseFunctor
       Udata(i, j, IP) = blast_pressure_in / (gamma0 - 1.0);
       Udata(i, j, IU) = 0.0;
       Udata(i, j, IV) = 0.0;
+
+      volume += dx * dy;
     }
     else
     {
@@ -1452,9 +1485,44 @@ class InitBlastFunctor : public HydroBaseFunctor
       Udata(i, j, IV) = 0.0;
     }
 
-  } // end operator ()
+  } // end operator () - TagRegularInit
+
+  KOKKOS_INLINE_FUNCTION
+  void
+  operator()(TagCorrectTotalEnergy const &, 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;
+
+    const real_t gamma0 = params.settings.gamma0;
+
+    // blast problem parameters
+    const real_t blast_radius = params.blast_radius;
+    const real_t radius2 = blast_radius * blast_radius;
+    const real_t blast_center_x = params.blast_center_x;
+    const real_t blast_center_y = params.blast_center_y;
+    const real_t blast_total_energy_inside = params.blast_total_energy_inside;
+
+    real_t x = xmin + dx / 2 + (i - ghostWidth) * dx;
+    real_t y = ymin + dy / 2 + (j - ghostWidth) * dy;
+
+    real_t d2 =
+      (x - blast_center_x) * (x - blast_center_x) + (y - blast_center_y) * (y - blast_center_y);
+
+    if (d2 < radius2)
+    {
+      Udata(i, j, IP) = blast_total_energy_inside / m_volume_inside;
+    }
+
+  } // end operator () - TagCorrectTotalEnergy
 
   DataArray_t Udata;
+  real_t      m_volume_inside = 0.0;
 
 }; // InitBlastFunctor