Implement initial density perturbations

.gitmodules

@@ -7,3 +7,6 @@
 [submodule "external/catch2"]
 	path = external/catch2
 	url = https://github.com/catchorg/Catch2.git
+[submodule "external/HighFive"]
+	path = external/HighFive
+	url = https://github.com/BlueBrain/HighFive.git

CMakeLists.txt

@@ -17,6 +17,9 @@ endif()
 
 project(athenaPK LANGUAGES CXX)
 
+# Adding HighFive library
+add_subdirectory(external/HighFive)
+
 set(CMAKE_CXX_EXTENSIONS OFF)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 set(CMAKE_CXX_STANDARD 17)
@@ -83,5 +86,4 @@ add_subdirectory(src)
 if (AthenaPK_ENABLE_TESTING)
   include(CTest)
   add_subdirectory(tst/regression)
-endif()
-
+endif()

src/eos/adiabatic_glmmhd.cpp

@@ -41,6 +41,7 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData<Real> *md) const {
   const auto nhydro = pkg->Param<int>("nhydro");
   const auto nscalars = pkg->Param<int>("nscalars");
 
+
   auto this_on_device = (*this);
 
   parthenon::par_for(
@@ -50,7 +51,23 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData<Real> *md) const {
         const auto &cons = cons_pack(b);
         auto &prim = prim_pack(b);
         // auto &nu = entropy_pack(b);
+
+        // Adding the gks gjs gis
+
+        // Getting the global indexing
+
+        auto pmb = md->GetBlockData(b)->GetBlockPointer();
+        auto pm = pmb->pmy_mesh;
+        auto hydro_pkg = pmb->packages.Get("Hydro");
+
+        const auto gis = pmb->loc.lx1 * pmb->block_size.nx1;
+        const auto gjs = pmb->loc.lx2 * pmb->block_size.nx2;
+        const auto gks = pmb->loc.lx3 * pmb->block_size.nx3;
+
+        // ...
+
+        //
 
-        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
+        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs,gis);
       });
 }

src/eos/adiabatic_glmmhd.hpp

@@ -30,10 +30,10 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
         gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
-
+  
   KOKKOS_INLINE_FUNCTION
   Real GetGamma() const { return gamma_; }
-
+  
   //----------------------------------------------------------------------------------------
   // \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO])
   // \brief returns adiabatic sound speed given vector of primitive variables
@@ -53,15 +53,15 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
     return std::sqrt(0.5 * (qsq + std::sqrt(tmp * tmp + 4.0 * asq * ct2)) / d);
   }
   //
-
+  
   //----------------------------------------------------------------------------------------
   // \!fn Real EquationOfState::ConsToPrim(View4D cons, View4D prim, const int& k, const
   // int& j, const int& i) \brief Fills an array of primitives given an array of
   // conserveds, potentially updating the conserved with floors
   template <typename View4D>
   KOKKOS_INLINE_FUNCTION void ConsToPrim(View4D cons, View4D prim, const int &nhydro,
                                          const int &nscalars, const int &k, const int &j,
-                                         const int &i) const {
+                                         const int &i, const int &gks, const int &gjs, const int &gis) const {
     auto gam = GetGamma();
     auto gm1 = gam - 1.0;
     auto density_floor_ = GetDensityFloor();

src/eos/adiabatic_hydro.cpp

@@ -30,26 +30,40 @@ using parthenon::ParArray4D;
 //           Container<Real> &rc,
 //           int il, int iu, int jl, int ju, int kl, int ku)
 // \brief Converts conserved into primitive variables in adiabatic hydro.
+
 void AdiabaticHydroEOS::ConservedToPrimitive(MeshData<Real> *md) const {
   auto const cons_pack = md->PackVariables(std::vector<std::string>{"cons"});
   auto prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
   auto ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire);
   auto jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire);
   auto kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire);
-
+  
   auto pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
   const auto nhydro = pkg->Param<int>("nhydro");
   const auto nscalars = pkg->Param<int>("nscalars");
-
+  
   auto this_on_device = (*this);
-
+  
   parthenon::par_for(
       DEFAULT_LOOP_PATTERN, "ConservedToPrimitive", parthenon::DevExecSpace(), 0,
       cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+
+        // Getting the global indexing
+
+        auto pmb = md->GetBlockData(b)->GetBlockPointer();
+        auto pm = pmb->pmy_mesh;
+        auto hydro_pkg = pmb->packages.Get("Hydro");
+
+        const auto gis = pmb->loc.lx1 * pmb->block_size.nx1;
+        const auto gjs = pmb->loc.lx2 * pmb->block_size.nx2;
+        const auto gks = pmb->loc.lx3 * pmb->block_size.nx3;
+
+        // ...
+
         const auto &cons = cons_pack(b);
         auto &prim = prim_pack(b);
-
-        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
+        
+        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs ,gis);
       });
 }

src/eos/adiabatic_hydro.hpp

@@ -51,15 +51,16 @@ class AdiabaticHydroEOS : public EquationOfState {
   template <typename View4D>
   KOKKOS_INLINE_FUNCTION void ConsToPrim(View4D cons, View4D prim, const int &nhydro,
                                          const int &nscalars, const int &k, const int &j,
-                                         const int &i) const {
+                                         const int &i, const int &gks, const int &gjs, const int &gis) const {
+
     Real gm1 = GetGamma() - 1.0;
     auto density_floor_ = GetDensityFloor();
     auto pressure_floor_ = GetPressureFloor();
     auto e_floor_ = GetInternalEFloor();
-
+    
     auto velocity_ceiling_ = GetVelocityCeiling();
     auto e_ceiling_ = GetInternalECeiling();
-
+    
     Real &u_d = cons(IDN, k, j, i);
     Real &u_m1 = cons(IM1, k, j, i);
     Real &u_m2 = cons(IM2, k, j, i);
@@ -71,7 +72,11 @@ class AdiabaticHydroEOS : public EquationOfState {
     Real &w_vy = prim(IV2, k, j, i);
     Real &w_vz = prim(IV3, k, j, i);
     Real &w_p = prim(IPR, k, j, i);
-
+
+    if (u_d <= 0.0) {std::cout << "(k,j,i)=(" << k << "," << j << "," << i << ") \n" ;
+                     std::cout << "(gks,gjs,gis)=(" << gks << "," << gjs << "," << gis << ") \n" ;
+                   }
+
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative density is encountered.
     PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0,
@@ -88,7 +93,7 @@ class AdiabaticHydroEOS : public EquationOfState {
 
     Real e_k = 0.5 * di * (SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
     w_p = gm1 * (u_e - e_k);
-
+    
     // apply velocity ceiling. By default ceiling is std::numeric_limits<Real>::infinity()
     const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz);
     if (w_v2 > SQR(velocity_ceiling_)) {
@@ -105,11 +110,29 @@ class AdiabaticHydroEOS : public EquationOfState {
       u_e -= e_k - e_k_new;
       e_k = e_k_new;
     }
-
+    
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative pressure is encountered.
+
+    // Trying to understand why negative pressure appear in tests
+
+    if (w_p <= 0.0) {std::cout << "w_p = gm1 * (u_e - e_k)";
+                    std::cout << "w_p=" << w_p << "\n" ;
+                    std::cout << "gm1=" << gm1 << "\n" ;
+                    std::cout << "u_e=" << u_e << "\n" ;
+                    std::cout << "e_k=" << e_k << "\n" ;
+                    std::cout << "(k,j,i)=(" << k << "," << j << "," << i << ") \n" ;
+                    std::cout << "(gks,gjs,gis)=(" << gks << "," << gjs << "," << gis << ") \n" ;
+                   }
+
+    //if (pressure_floor_ < 0.0) {std::cout << "pressure_floor_" << pressure_floor_ << "\n" ;}
+    //if (e_floor_ < 0.0) {std::cout << "e_floor_" << e_floor_ << "\n" ;}
+
+    // The first argument check whether one of the conditions is true
+    // By default, floors are deactivated, ie. pressure floor and e_floor
+    // are -1. The code will eventually crash when w_p turns < 0.
     PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
-                      "Got negative pressure. Consider enabling first-order flux "
+                      "Wow. Got negative pressure. Consider enabling first-order flux "
                       "correction or setting a reasonble pressure or temperature floor.");
 
     // Pressure floor (if present) takes precedence over temperature floor