Fixed parthenon typo

parthenon-hpc-lab · Oct 12, 2023 · 8a8ba8f · 8a8ba8f
1 parent f188f33
commit 8a8ba8f
Show file tree

Hide file tree

Showing 9 changed files with 72 additions and 145 deletions.
diff --git a/src/eos/adiabatic_glmmhd.cpp b/src/eos/adiabatic_glmmhd.cpp
@@ -52,22 +52,14 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData<Real> *md) const {
         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, gks, gjs,gis);
+
+        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
       });
 }
diff --git a/src/eos/adiabatic_glmmhd.hpp b/src/eos/adiabatic_glmmhd.hpp
@@ -61,7 +61,7 @@ class AdiabaticGLMMHDEOS : 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 int &gks, const int &gjs, const int &gis) const {
+                                         const int &i) const {
     auto gam = GetGamma();
     auto gm1 = gam - 1.0;
     auto density_floor_ = GetDensityFloor();

diff --git a/src/eos/adiabatic_hydro.cpp b/src/eos/adiabatic_hydro.cpp
@@ -55,15 +55,9 @@ void AdiabaticHydroEOS::ConservedToPrimitive(MeshData<Real> *md) const {
         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, gks, gjs ,gis);
+        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
       });
 }
diff --git a/src/eos/adiabatic_hydro.hpp b/src/eos/adiabatic_hydro.hpp
@@ -51,7 +51,7 @@ 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 int &gks, const int &gjs, const int &gis) const {
+                                         const int &i) const {
 
     Real gm1 = GetGamma() - 1.0;
     auto density_floor_ = GetDensityFloor();
@@ -72,10 +72,6 @@ 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.
@@ -113,21 +109,7 @@ class AdiabaticHydroEOS : public EquationOfState {
 
     // 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.

diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp
@@ -371,7 +371,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
   }
 
   /************************************************************
-   * Read Velocity perturbation
+   * Read Density perturbation
    ************************************************************/
 
   const auto mu_rho = pin->GetOrAddReal("problem/cluster/init_perturb", "mu_rho", 0.0); // Mean density of perturbations
@@ -401,7 +401,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
     // this way is easier for now)
     Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy},
                std::vector<int>({3}));
-    hydro_pkg->AddField("tmp_perturb_rho", m);
+    hydro_pkg->AddField("tmp_perturb", m);
 
   }
 
@@ -711,7 +711,6 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
   auto const &cons = md->PackVariables(std::vector<std::string>{"cons"});
   const auto num_blocks = md->NumBlocks();    
 
-
   /************************************************************
    * Set initial density perturbations (read from HDF5 file)
    ************************************************************/
@@ -721,53 +720,11 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
   const Real thickness_ism    = pin->GetOrAddReal("problem/cluster/init_perturb", "thickness_ism", 0.0);
   const bool overpressure_ring = pin->GetOrAddBoolean("problem/cluster/init_perturb", "overpressure_ring", false);
   const bool spherical_collapse  = pin->GetOrAddBoolean("problem/cluster/init_perturb", "spherical_collapse", false);
-  const bool numerical_diffusion = pin->GetOrAddBoolean("problem/cluster/init_perturb", "numerical_diffusion", false);
-
+
   hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho);
 
   Real passive_scalar = 0.0; // Not useful here
 
-  // Numerical diffusion problem
-
-  if (numerical_diffusion == true) {
-
-    pmb->par_reduce(
-        "Init density field", 0, num_blocks - 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, Real &lsum) {
-
-          auto pmbb  = md->GetBlockData(b)->GetBlockPointer(); // Meshblock b
-
-          const auto gis = pmbb->loc.lx1 * pmb->block_size.nx1;
-          const auto gjs = pmbb->loc.lx2 * pmb->block_size.nx2;
-          const auto gks = pmbb->loc.lx3 * pmb->block_size.nx3;
-
-          const auto &coords = cons.GetCoords(b);
-          const auto &u = cons(b);
-
-          const Real x = coords.Xc<1>(i);
-          const Real y = coords.Xc<2>(j);
-          const Real z = coords.Xc<3>(k);
-          const Real background_density  = pin->GetOrAddReal("problem/cluster/init_perturb", "background_density", 3);
-          const Real foreground_density  = pin->GetOrAddReal("problem/cluster/init_perturb", "foreground_density", 150);
-
-          // Setting density for the left side of the box
-          if (x <= 0) {
-
-            u(IDN, k, j, i) = background_density;
-
-          }
-
-          else {
-
-            u(IDN, k, j, i) = foreground_density;  
-
-          }
-
-        },
-        passive_scalar);
-
-  }
-
   // Spherical collapse test with an initial overdensity
 
   if (spherical_collapse == true) {
@@ -913,55 +870,50 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
         passive_scalar);
 
   }
-    
+
   /************************************************************
    * Setting up a perturbed density field (hardcoded version)
    ************************************************************/
 
   const auto mu_rho = pin->GetOrAddReal("problem/cluster/init_perturb", "mu_rho", 0.0);
-
+  const auto background_rho = pin->GetOrAddReal("problem/cluster/init_perturb", "background_rho", 0.0);
+
   if (mu_rho != 0.0) {
 
-    std::cout << "Entering rho perturbation mode \n" ;
-
     auto few_modes_ft_rho = hydro_pkg->Param<FewModesFTLog>("cluster/few_modes_ft_rho");
-
-    std::cout << "few_modes_ft_rho defined \n" ;
 
     // Init phases on all blocks
     for (int b = 0; b < md->NumBlocks(); b++) {
       auto pmb = md->GetBlockData(b)->GetBlockPointer();
       few_modes_ft_rho.SetPhases(pmb.get(), pin);
 
     }
-
-    std::cout << "Phase is set on each block \n" ;
 
     // As for t_corr in few_modes_ft, the choice for dt is
     // in principle arbitrary because the inital v_hat is 0 and the v_hat_new will contain
     // the perturbation (and is normalized in the following to get the desired sigma_v)
-    //const Real dt = 1.0;
-    //few_modes_ft_rho.Generate(md, dt, "tmp_perturb_rho");
+    const Real dt = 1.0;
+    few_modes_ft_rho.Generate(md, dt, "tmp_perturb");
 
-    //Real v2_sum_rho = 0.0; // used for normalization
+    Real v2_sum_rho = 0.0; // used for normalization
 
-    //auto perturb_pack_rho = md->PackVariables(std::vector<std::string>{"tmp_perturb_rho"});
+    auto perturb_pack_rho = md->PackVariables(std::vector<std::string>{"tmp_perturb"});
 
-    //pmb->par_reduce(
-    //    "Init sigma_v", 0, num_blocks - 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, Real &lsum) {
-    //      const auto &coords = cons.GetCoords(b);
-    //      const auto &u = cons(b);
-            
-    //      u(IDN, k, j, i) = 1000 + perturb_pack_rho(b, 0, k, j, i);
-    //      std::cout << "Rho value:" << perturb_pack_rho(b, 0, k, j, i) << "\n" ;
-    //    },
-    //    v2_sum_rho);
-
-//#ifdef MPI_PARALLEL
-    //PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &v2_sum_rho, 1, MPI_PARTHENON_REAL,
-    //                                  MPI_SUM, MPI_COMM_WORLD));
-//#endif // MPI_PARALLEL
+    pmb->par_reduce(
+        "Init sigma_v", 0, num_blocks - 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, Real &lsum) {
+          const auto &coords = cons.GetCoords(b);
+          const auto &u = cons(b);
+
+          u(IDN, k, j, i) = background_rho + mu_rho * std::abs(perturb_pack_rho(b, 0, k, j, i));
+          std::cout << "Rho value:" << perturb_pack_rho(b, 0, k, j, i) << "\n" ;
+        },
+        v2_sum_rho);
+
+#ifdef MPI_PARALLEL
+    PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &v2_sum_rho, 1, MPI_PARTHENON_REAL,
+                                      MPI_SUM, MPI_COMM_WORLD));
+#endif // MPI_PARALLEL
 
   }
 

diff --git a/src/pgen/cluster/agn_triggering.cpp b/src/pgen/cluster/agn_triggering.cpp
@@ -188,7 +188,7 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass,
               AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(),
                                              k, j, i);
               // Update the Primitives
-              eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0);
+              eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
             }
           }
         }
@@ -309,7 +309,7 @@ void AGNTriggering::RemoveBondiAccretedGas(parthenon::MeshData<parthenon::Real>
                                          i);
 
           // Update the Primitives
-          eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0);
+          eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
         }
       });
 }

diff --git a/src/pgen/cluster/snia_feedback.cpp b/src/pgen/cluster/snia_feedback.cpp
@@ -114,7 +114,7 @@ void SNIAFeedback::FeedbackSrcTerm(parthenon::MeshData<parthenon::Real> *md,
         cons(IEN, k, j, i) += snia_energy_density;
         AddDensityToConsAtFixedVel(snia_mass_density, cons, prim, k, j, i);
 
-        eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0);
+        eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
       });
 }
 

diff --git a/src/utils/few_modes_ft.cpp b/src/utils/few_modes_ft.cpp
@@ -115,9 +115,9 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) {
   // below take the logical grid size into account. For example, the local phases at level
   // 1 should be calculated assuming a grid that is twice as large as the root grid.
   const auto root_level = pm->GetRootLevel();
-  auto gnx1 = pm->mesh_size.nx1 * std::pow(2, pmb->loc.level() - root_level);
-  auto gnx2 = pm->mesh_size.nx2 * std::pow(2, pmb->loc.level() - root_level);
-  auto gnx3 = pm->mesh_size.nx3 * std::pow(2, pmb->loc.level() - root_level);
+  auto gnx1 = pm->mesh_size.nx1 * std::pow(2, pmb->loc.level - root_level);
+  auto gnx2 = pm->mesh_size.nx2 * std::pow(2, pmb->loc.level - root_level);
+  auto gnx3 = pm->mesh_size.nx3 * std::pow(2, pmb->loc.level - root_level);
 
   // Restriction should also be easily fixed, just need to double check transforms and
   // volume weighting everywhere
@@ -130,9 +130,9 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) {
   const auto nx2 = pmb->block_size.nx2;
   const auto nx3 = pmb->block_size.nx3;
 
-  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 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;
 
   // make local ref to capure in lambda
   const auto num_modes = num_modes_;