From 8a8ba8f859573a3da570660eca2798749440d546 Mon Sep 17 00:00:00 2001 From: HPC user Date: Thu, 12 Oct 2023 11:45:36 +0200 Subject: [PATCH] Fixed parthenon typo --- src/eos/adiabatic_glmmhd.cpp | 12 +--- src/eos/adiabatic_glmmhd.hpp | 2 +- src/eos/adiabatic_hydro.cpp | 8 +-- src/eos/adiabatic_hydro.hpp | 22 +------ src/pgen/cluster.cpp | 98 +++++++--------------------- src/pgen/cluster/agn_triggering.cpp | 4 +- src/pgen/cluster/snia_feedback.cpp | 2 +- src/utils/few_modes_ft.cpp | 12 ++-- src/utils/few_modes_ft_lognormal.cpp | 57 +++++++++------- 9 files changed, 72 insertions(+), 145 deletions(-) diff --git a/src/eos/adiabatic_glmmhd.cpp b/src/eos/adiabatic_glmmhd.cpp index 124fd2a7..28c0984a 100644 --- a/src/eos/adiabatic_glmmhd.cpp +++ b/src/eos/adiabatic_glmmhd.cpp @@ -52,22 +52,14 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData *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 index 9b59f3e2..13b60373 100644 --- a/src/eos/adiabatic_glmmhd.hpp +++ b/src/eos/adiabatic_glmmhd.hpp @@ -61,7 +61,7 @@ class AdiabaticGLMMHDEOS : public EquationOfState { template 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 index 475a754c..651000e4 100644 --- a/src/eos/adiabatic_hydro.cpp +++ b/src/eos/adiabatic_hydro.cpp @@ -55,15 +55,9 @@ void AdiabaticHydroEOS::ConservedToPrimitive(MeshData *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 index 0addb9ee..d51adf4b 100644 --- a/src/eos/adiabatic_hydro.hpp +++ b/src/eos/adiabatic_hydro.hpp @@ -51,7 +51,7 @@ class AdiabaticHydroEOS : public EquationOfState { template 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 index 20eaf903..06150d21 100644 --- 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({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 *md) { auto const &cons = md->PackVariables(std::vector{"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 *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,20 +870,17 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *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("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++) { @@ -934,34 +888,32 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { 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{"tmp_perturb_rho"}); + auto perturb_pack_rho = md->PackVariables(std::vector{"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 index 5f39957e..242de4fa 100644 --- 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 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 index 25087235..ed41657f 100644 --- a/src/pgen/cluster/snia_feedback.cpp +++ b/src/pgen/cluster/snia_feedback.cpp @@ -114,7 +114,7 @@ void SNIAFeedback::FeedbackSrcTerm(parthenon::MeshData *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 index f477b211..f80fd103 100644 --- 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_; diff --git a/src/utils/few_modes_ft_lognormal.cpp b/src/utils/few_modes_ft_lognormal.cpp index d47b30c9..ffdbc17d 100644 --- a/src/utils/few_modes_ft_lognormal.cpp +++ b/src/utils/few_modes_ft_lognormal.cpp @@ -33,9 +33,11 @@ FewModesFTLog::FewModesFTLog(parthenon::ParameterInput *pin, parthenon::StateDes std::string prefix, int num_modes, ParArray2D k_vec, Real k_min, Real k_max, Real sol_weight, Real t_corr, uint32_t rseed, bool fill_ghosts) - : prefix_(prefix), num_modes_(num_modes), k_min_(k_min), k_max_(k_max), + : prefix_(prefix), num_modes_(num_modes), k_vec_(k_vec), k_min_(k_min), k_max_(k_max), t_corr_(t_corr), fill_ghosts_(fill_ghosts) { + //std::cout << "FewModesFTLog: k_vec_(0,0)=" << + if ((num_modes > 100) && (parthenon::Globals::my_rank == 0)) { std::cout << "### WARNING using more than 100 explicit modes will significantly " << "increase the runtime." << std::endl @@ -49,11 +51,13 @@ FewModesFTLog::FewModesFTLog(parthenon::ParameterInput *pin, parthenon::StateDes // Need to make this comparison on the host as (for some reason) an extended cuda device // lambda cannot live in the constructor of an object. auto k_vec_host = k_vec.GetHostMirrorAndCopy(); + for (int i = 0; i < num_modes; i++) { + PARTHENON_REQUIRE(std::abs(k_vec_host(0, i)) <= gnx1 / 2, "k_vec x1 mode too large"); PARTHENON_REQUIRE(std::abs(k_vec_host(1, i)) <= gnx2 / 2, "k_vec x2 mode too large"); PARTHENON_REQUIRE(std::abs(k_vec_host(2, i)) <= gnx3 / 2, "k_vec x3 mode too large"); - } + } const auto nx1 = pin->GetInteger("parthenon/meshblock", "nx1"); const auto nx2 = pin->GetInteger("parthenon/meshblock", "nx2"); @@ -88,6 +92,7 @@ FewModesFTLog::FewModesFTLog(parthenon::ParameterInput *pin, parthenon::StateDes } void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { + auto pm = pmb->pmy_mesh; auto hydro_pkg = pmb->packages.Get("Hydro"); @@ -102,8 +107,6 @@ void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { // separate and not touch the main driver at the expense of using one pack per rank -- // which is typically fastest on devices anyway. - std::cout << "Entering SetPhases routine \n" ; - const auto pack_size = pin->GetInteger("parthenon/mesh", "pack_size"); PARTHENON_REQUIRE_THROWS(pack_size == -1, "Few modes FT currently needs parthenon/mesh/pack_size=-1 " @@ -140,8 +143,6 @@ void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { // make local ref to capure in lambda const auto num_modes = num_modes_; auto &k_vec = k_vec_; - - std::cout << "SetPhases: k_vec(0,0) = " << k_vec(0,0) ; Complex I(0.0, 1.0); @@ -150,8 +151,6 @@ void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { auto &phases_j = base->Get(prefix_ + "_phases_j").data; auto &phases_k = base->Get(prefix_ + "_phases_k").data; - std::cout << "SetPhases: all constants defined, starting the real shit \n" ; - const auto ng = fill_ghosts_ ? parthenon::Globals::nghost : 0; pmb->par_for( "FMFT: calc phases_i", 0, nx1 - 1 + 2 * ng, KOKKOS_LAMBDA(int i) { @@ -159,22 +158,13 @@ void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { Real w_kx; Complex phase; - std::cout << "SetPhases: first loop \n" ; - for (int m = 0; m < num_modes; m++) { - std::cout << "SetPhases: entering for loop on modes \n" ; - - std::cout << "SetPhases: defining w_kx \n" ; - std::cout << "SetPhases: k_vec(0, m)" << k_vec(0, m) << "\n" ; - std::cout << "gnx1:" << static_cast(gnx1) << " \n" ; w_kx = k_vec(0, m) * 2. * M_PI / static_cast(gnx1); - std::cout << "SetPhases: w_kx defined \n" ; - // adjust phase factor to Complex->Real IFT: u_hat*(k) = u_hat(-k) - std::cout << "SetPhases: about to compute phase \n" ; + if (k_vec(0, m) == 0.0) { phase = 0.5 * Kokkos::exp(I * w_kx * gi); } else { @@ -183,7 +173,6 @@ void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { phases_i(i, m, 0) = phase.real(); phases_i(i, m, 1) = phase.imag(); - std::cout << "SetPhases: phases computed \n" ; } }); @@ -382,19 +371,37 @@ ParArray2D MakeRandomModesLog(const int num_modes, const Real k_min, const std::mt19937 rng; rng.seed(rseed); - std::uniform_int_distribution<> dist(-k_high, k_high); // Function to be used to define size of the + std::uniform_real_distribution<> dist(-k_high, k_high); + std::uniform_real_distribution<> distlog(std::log10(k_low), std::log10(k_high)); int n_mode = 0; int n_attempt = 0; constexpr int max_attempts = 1000000; - Real kx1, kx2, kx3, k_mag, ampl; + Real kx1, kx2, kx3, skx1, skx2, skx3, lkx1, lkx2, lkx3, k_mag, ampl; bool mode_exists = false; while (n_mode < num_modes && n_attempt < max_attempts) { n_attempt += 1; - - kx1 = dist(rng); - kx2 = dist(rng); - kx3 = dist(rng); + + skx1 = dist(rng); + skx2 = dist(rng); + skx3 = dist(rng); + + skx1 = skx1 / std::abs(skx1); + skx2 = skx2 / std::abs(skx2); + skx3 = skx3 / std::abs(skx3); + + lkx1 = distlog(rng); + lkx2 = distlog(rng); + lkx3 = distlog(rng); + + //kx1 = skx1 * std::floor(std::pow(10,lkx1)); + //kx2 = skx2 * std::floor(std::pow(10,lkx2)); + //kx3 = skx3 * std::floor(std::pow(10,lkx3)); + + kx1 = std::floor(std::pow(10,lkx1)); + kx2 = std::floor(std::pow(10,lkx2)); + kx3 = std::floor(std::pow(10,lkx3)); + k_mag = std::sqrt(SQR(kx1) + SQR(kx2) + SQR(kx3)); // Expected amplitude of the spectral function. If this is changed, it also needs to