diff --git a/src/eos/adiabatic_hydro.cpp b/src/eos/adiabatic_hydro.cpp index fa0af7a7..a59c8c90 100644 --- a/src/eos/adiabatic_hydro.cpp +++ b/src/eos/adiabatic_hydro.cpp @@ -48,12 +48,6 @@ void AdiabaticHydroEOS::ConservedToPrimitive(MeshData *md) const { 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 &cons = cons_pack(b); auto &prim = prim_pack(b); diff --git a/src/eos/adiabatic_hydro.hpp b/src/eos/adiabatic_hydro.hpp index f4d714f4..d5a0debf 100644 --- a/src/eos/adiabatic_hydro.hpp +++ b/src/eos/adiabatic_hydro.hpp @@ -114,7 +114,7 @@ class AdiabaticHydroEOS : public EquationOfState { // 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, - "Wow. Got negative pressure. Consider enabling first-order flux " + "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 diff --git a/src/hydro/srcterms/gravitational_field.hpp b/src/hydro/srcterms/gravitational_field.hpp index e0ae98af..31feb1fe 100644 --- a/src/hydro/srcterms/gravitational_field.hpp +++ b/src/hydro/srcterms/gravitational_field.hpp @@ -37,21 +37,6 @@ void GravitationalFieldSrcTerm(parthenon::MeshData *md, IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); - /* - std::cout << "Values of gravitationnal field" << std::endl; - std::cout - << gravitationalField.g_from_r(1e-4) << ", " << gravitationalField.g_from_r(3.16e-4) - << ", " << gravitationalField.g_from_r(1e-3) << ", " << - gravitationalField.g_from_r(3.16e-3) - << ", " << gravitationalField.g_from_r(1e-2) << ", " << - gravitationalField.g_from_r(3.16e-2) - << ", " << gravitationalField.g_from_r(1e-1) << ", " << - gravitationalField.g_from_r(3.16e-1) - << ", " << gravitationalField.g_from_r(1.0) << ", " << - gravitationalField.g_from_r(3.16) - << ", " << gravitationalField.g_from_r(10.0) << std::endl; - */ - parthenon::par_for( DEFAULT_LOOP_PATTERN, "GravitationalFieldSrcTerm", parthenon::DevExecSpace(), 0, cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 798a7f22..34809aa9 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -44,7 +44,6 @@ #include "../hydro/srcterms/tabular_cooling.hpp" #include "../main.hpp" #include "../utils/few_modes_ft.hpp" -#include "../utils/few_modes_ft_lognormal.hpp" // Cluster headers #include "cluster/agn_feedback.hpp" @@ -62,7 +61,6 @@ namespace cluster { using namespace parthenon::driver::prelude; using namespace parthenon::package::prelude; using utils::few_modes_ft::FewModesFT; -using utils::few_modes_ft_log::FewModesFTLog; void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt) { @@ -365,46 +363,6 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddField("plasma_beta", m); } - /************************************************************ - * Read Density perturbation - ************************************************************/ - - const auto mu_rho = pin->GetOrAddReal("problem/cluster/init_perturb", "mu_rho", - 0.0); // Mean density of perturbations - - if (mu_rho != 0.0) { - - auto k_min_rho = pin->GetReal("problem/cluster/init_perturb", - "k_min_rho"); // Minimum wavenumber of perturbation - auto num_modes_rho = - pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_rho", 40); - auto sol_weight_rho = - pin->GetOrAddReal("problem/cluster/init_perturb", "sol_weight_rho", 1.0); - uint32_t rseed_rho = - pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_rho", 1); - - // Computing the kmax ie. the Nyquist limit - auto grid_ni = pin->GetOrAddInteger( - "parthenon/mesh", "nx1", 64); // Assuming cubic grid with equal size in each axis - auto k_max_rho = grid_ni / 2; - - const auto t_corr_rho = 1e-10; - auto k_vec_rho = utils::few_modes_ft_log::MakeRandomModesLog( - num_modes_rho, k_min_rho, k_max_rho, rseed_rho); // Generating random modes - - auto few_modes_ft_rho = - FewModesFTLog(pin, hydro_pkg, "cluster_perturb_rho", num_modes_rho, k_vec_rho, - k_min_rho, k_max_rho, sol_weight_rho, t_corr_rho, rseed_rho); - - hydro_pkg->AddParam<>("cluster/few_modes_ft_rho", few_modes_ft_rho); - - // Add field for initial perturation (must not need to be consistent but defining it - // this way is easier for now) - Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, - std::vector({3})); - hydro_pkg->AddField("tmp_perturb", m); - } - /************************************************************ * Read Velocity perturbation ************************************************************/ @@ -545,8 +503,6 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { * Initialize the initial hydro state ************************************************************/ const auto &init_uniform_gas = hydro_pkg->Param("init_uniform_gas"); - const auto init_perseus = pin->GetOrAddBoolean( - "problem/cluster/hydrostatic_equilibrium", "init_perseus", false); if (init_uniform_gas) { const Real rho = hydro_pkg->Param("uniform_gas_rho"); @@ -570,49 +526,6 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { u(IM3, k, j, i) = Mz; u(IEN, k, j, i) = E; }); - } - if (init_perseus) { - - // initialize conserved variables - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::UniformGas", - parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { - // Units - const Real mh = units.mh(); - const Real k_boltzmann = units.k_boltzmann(); - - const Real mu = hydro_pkg->Param("mu"); - const Real mu_e = hydro_pkg->Param("mu_e"); - - // Calculate radius - const Real r = sqrt(coords.Xc<1>(i) * coords.Xc<1>(i) + - coords.Xc<2>(j) * coords.Xc<2>(j) + - coords.Xc<3>(k) * coords.Xc<3>(k)); - - const Real re3 = r * 1e3; // Mpc to kpc - - // Get pressure and density from generated profile - const Real ne_r = - 0.0192 * 1 / (1 + std::pow(re3 / 18, 3)) + - 0.046 * 1 / std::pow((1 + std::pow(re3 / 57, 2)), 1.8) + - 0.0048 * 1 / - std::pow((1 + std::pow(re3 / 200, 2)), 0.87); // cgs units (cm-3) - const Real T_r = - 8.12e7 * (1 + std::pow(re3 / 71, 3)) / (2.3 + std::pow(re3 / 71, 3)); // K - - const Real ne_r_cu = ne_r * std::pow(units.cm(), -3); // Code units, Mpc^-3 - - const Real rho_r = mu_e * mh * ne_r_cu; // From ne_r (conversion) - const Real P_r = - k_boltzmann * (mu_e / mu) * ne_r_cu * T_r; // From T_r (conversion) - - u(IDN, k, j, i) = rho_r; - u(IM1, k, j, i) = 0.0; - u(IM2, k, j, i) = 0.0; - u(IM3, k, j, i) = 0.0; - u(IEN, k, j, i) = P_r / gm1; - }); } else { @@ -640,7 +553,6 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real P_r = P_rho_profile.P_from_r(r); const Real rho_r = P_rho_profile.rho_from_r(r); - // u(IDN, k, j, i) = rho_r; u(IDN, k, j, i) = rho_r; u(IM1, k, j, i) = 0.0; u(IM2, k, j, i) = 0.0; @@ -1001,9 +913,6 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { * Set initial magnetic field perturbations (resets magnetic field field) ************************************************************/ const auto sigma_b = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_b", 0.0); - const auto alpha_b = - pin->GetOrAddReal("problem/cluster/init_perturb", "alpha_b", 2.0 / 3.0); - const auto r_scale = pin->GetOrAddReal("problem/cluster/init_perturb", "r_scale", 0.1); if (sigma_b != 0.0) { auto few_modes_ft = hydro_pkg->Param("cluster/few_modes_ft_b");