diff --git a/src/eos/adiabatic_glmmhd.hpp b/src/eos/adiabatic_glmmhd.hpp index 0a4cbc9f..ae073d97 100644 --- a/src/eos/adiabatic_glmmhd.hpp +++ b/src/eos/adiabatic_glmmhd.hpp @@ -93,9 +93,9 @@ class AdiabaticGLMMHDEOS : public EquationOfState { // 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, - "Got negative density. Consider enabling first-order flux " - "correction or setting a reasonable density floor."); + //PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0, + // "Got negative density. Consider enabling first-order flux " + // "correction or setting a reasonable density floor."); // apply density floor, without changing momentum or energy u_d = (u_d > density_floor_) ? u_d : density_floor_; @@ -134,9 +134,9 @@ class AdiabaticGLMMHDEOS : 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. - PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0, - "Got negative pressure. Consider enabling first-order flux " - "correction or setting a reasonable pressure or temperature floor."); + //PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0, + // "Got negative pressure. Consider enabling first-order flux " + // "correction or setting a reasonable pressure or temperature floor."); // Pressure floor (if present) takes precedence over temperature floor if ((pressure_floor_ > 0.0) && (w_p < pressure_floor_)) { diff --git a/src/hydro/srcterms/gravitational_field.hpp b/src/hydro/srcterms/gravitational_field.hpp index 31feb1fe..fd3bb3ca 100644 --- a/src/hydro/srcterms/gravitational_field.hpp +++ b/src/hydro/srcterms/gravitational_field.hpp @@ -36,7 +36,18 @@ void GravitationalFieldSrcTerm(parthenon::MeshData *md, IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); 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/cloud.cpp b/src/pgen/cloud.cpp index bc9fed25..086b9b38 100644 --- a/src/pgen/cloud.cpp +++ b/src/pgen/cloud.cpp @@ -143,6 +143,7 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) { // \brief Problem Generator for the cloud in wind setup void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { + auto hydro_pkg = pmb->packages.Get("Hydro"); auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); @@ -156,9 +157,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_FAIL("Requested to initialize magnetic fields by `cloud/plasma_beta > 0`, " "but `hydro/fluid` is not supporting MHD."); } - + auto steepness = pin->GetOrAddReal("problem/cloud", "cloud_steepness", 10); - + // initialize conserved variables auto &mbd = pmb->meshblock_data.Get(); auto &u_dev = mbd->Get("cons").data; @@ -191,13 +192,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { u(IM2, k, j, i) = mom; // Can use rhoe_wind here as simulation is setup in pressure equil. u(IEN, k, j, i) = rhoe_wind + 0.5 * mom * mom / rho; - + if (mhd_enabled) { u(IB1, k, j, i) = Bx; u(IB2, k, j, i) = By; u(IEN, k, j, i) += 0.5 * (Bx * Bx + By * By); } - + // Init passive scalars for (auto n = nhydro; n < nhydro + nscalars; n++) { if (rad <= r_cloud) { @@ -240,13 +241,14 @@ void InflowWindX2(std::shared_ptr> &mbd, bool coarse) { } parthenon::AmrTag ProblemCheckRefinementBlock(MeshBlockData *mbd) { + auto pmb = mbd->GetBlockPointer(); auto w = mbd->Get("prim").data; - + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); - + auto hydro_pkg = pmb->packages.Get("Hydro"); const auto nhydro = hydro_pkg->Param("nhydro"); diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 58c188a5..3f03755f 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -117,22 +117,22 @@ Real ClusterEstimateTimestep(MeshData *md) { //======================================================================================== void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg) { - + /************************************************************ * Read Uniform Gas ************************************************************/ - + const bool init_uniform_gas = pin->GetOrAddBoolean("problem/cluster/uniform_gas", "init_uniform_gas", false); hydro_pkg->AddParam<>("init_uniform_gas", init_uniform_gas); - + if (init_uniform_gas) { const Real uniform_gas_rho = pin->GetReal("problem/cluster/uniform_gas", "rho"); const Real uniform_gas_ux = pin->GetReal("problem/cluster/uniform_gas", "ux"); const Real uniform_gas_uy = pin->GetReal("problem/cluster/uniform_gas", "uy"); const Real uniform_gas_uz = pin->GetReal("problem/cluster/uniform_gas", "uz"); const Real uniform_gas_pres = pin->GetReal("problem/cluster/uniform_gas", "pres"); - + hydro_pkg->AddParam<>("uniform_gas_rho", uniform_gas_rho); hydro_pkg->AddParam<>("uniform_gas_ux", uniform_gas_ux); hydro_pkg->AddParam<>("uniform_gas_uy", uniform_gas_uy); @@ -238,15 +238,15 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd /************************************************************ * Read SNIA Feedback ************************************************************/ - + SNIAFeedback snia_feedback(pin, hydro_pkg); - + /************************************************************ * Read Stellar Feedback ************************************************************/ - + StellarFeedback stellar_feedback(pin, hydro_pkg); - + /************************************************************ * Read Clips (ceilings and floors) ************************************************************/ @@ -297,11 +297,11 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd std::string reduction_strs[] = {"stellar_mass", "added_dfloor_mass", "removed_eceil_energy", "removed_vceil_energy", "added_vAceil_mass"}; - + // Add a param for each reduction, then add it as a summation reduction for // history outputs auto hst_vars = hydro_pkg->Param(parthenon::hist_param_key); - + for (auto reduction_str : reduction_strs) { hydro_pkg->AddParam(reduction_str, 0.0, true); hst_vars.emplace_back(parthenon::HistoryOutputVar( @@ -541,7 +541,8 @@ 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_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"); @@ -565,7 +566,48 @@ 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 { /************************************************************ @@ -722,7 +764,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); - auto hydro_pkg = pmb->packages.Get("Hydro"); + auto hydro_pkg = pmb->packages.Get("Hydro"); const auto fluid = hydro_pkg->Param("fluid"); auto const &cons = md->PackVariables(std::vector{"cons"}); const auto num_blocks = md->NumBlocks(); @@ -732,11 +774,75 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { ************************************************************/ const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); - - hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + const bool spherical_cloud = pin->GetOrAddBoolean("problem/cluster/init_perturb", "spherical_cloud", false); + const bool cluster_cloud = pin->GetOrAddBoolean("problem/cluster/init_perturb", "cluster_cloud", false); Real passive_scalar = 0.0; // Not useful here + if (cluster_cloud == true) { + + const Real x_cloud = pin->GetOrAddReal("problem/cluster/init_perturb", "x_cloud", 0.0); // 0.0 pc + const Real y_cloud = pin->GetOrAddReal("problem/cluster/init_perturb", "y_cloud", 5e-2); // 100 pc + const Real r_cloud = pin->GetOrAddReal("problem/cluster/init_perturb", "r_cloud", 1e-4); // 100 pc + const Real rho_cloud = pin->GetOrAddReal("problem/cluster/init_perturb", "rho_cloud", 150.0); // 1e-24 g/cm^3 + const Real steepness = pin->GetOrAddReal("problem/cluster/init_perturb", "steep_cloud", 10.0); + + Real passive_scalar = 0.0; // Useless + + /* + for (int b = 0; b < md->NumBlocks(); b++) { + + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + + // Initialize the conserved variables + auto &u = pmb->meshblock_data.Get()->Get("cons").data; + auto &coords = pmb->coords; + + Real r_min = 1e6; + + for (int k = kb.s; k <= kb.e; k++) { + for (int j = jb.s; j <= jb.e; j++) { + for (int i = ib.s; i <= ib.e; i++) { + + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); + + if (r < r_min) {r_min = r;} + + } + } + } + if (r_min < 20 * r_cloud) {std::cout << "Should be refining now" << std::endl;} + if (r_min < 20 * r_cloud) { + parthenon::AmrTag::refine; + parthenon::AmrTag::same;} + } + */ + + 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 &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 r = std::sqrt(SQR(x - x_cloud) + SQR(y - y_cloud) + SQR(z)); + + Real rho = rho_cloud * (1.0 - std::tanh(steepness * (r / r_cloud - 1.0))); + u(IDN, k, j, i) += rho; + + }, + passive_scalar); + + } + /* -------------- Setting up a clumpy atmosphere -------------- 1) Extract the values of the density from an input hdf5 file using H5Easy @@ -750,7 +856,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { auto filename_rho = pin->GetOrAddString("problem/cluster/init_perturb", "init_perturb_rho_file","none"); const Real perturb_amplitude = pin->GetOrAddReal("problem/cluster/init_perturb", "perturb_amplitude", 1); - const Real box_size_over_two = pin->GetOrAddReal("parthenon/mesh", "x1max", 0.250); + const Real box_size_over_two = pin->GetOrAddReal("problem/cluster/init_perturb", "xmax", 0.250); + + // Loading files std::string keys_rho = "data"; H5Easy::File file(filename_rho, HighFive::File::ReadOnly); @@ -762,6 +870,8 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { std::cout << "Entering initialisation of rho field"; + // Setting up the perturbations + 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) { @@ -775,14 +885,23 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real y = coords.Xc<2>(j); const Real z = coords.Xc<3>(k); + const Real r = sqrt(SQR(x) + SQR(y) + SQR(z)); + // Getting the corresponding index in the const Real rho_init_index_x = floor((x + box_size_over_two) / (2 * box_size_over_two) * (rho_init_size - 1)); const Real rho_init_index_y = floor((y + box_size_over_two) / (2 * box_size_over_two) * (rho_init_size - 1)); const Real rho_init_index_z = floor((z + box_size_over_two) / (2 * box_size_over_two) * (rho_init_size - 1)); + //const Real damping = 1 - std::exp(r - box_size_over_two); + + if ((rho_init_index_x >= 0) && (rho_init_index_x <= rho_init_size - 1) && (rho_init_index_y >= 0) && (rho_init_index_y <= rho_init_size - 1) && (rho_init_index_z >= 0) && (rho_init_index_z <= rho_init_size - 1)) { + // Case where the box is filled with perturbations of equal mean amplitude - u(IDN, k, j, i) += perturb_amplitude * rho_init[rho_init_index_x][rho_init_index_y][rho_init_index_z] * (u(IDN, k, j, i) / 29.6); + //u(IDN, k, j, i) += perturb_amplitude * rho_init[rho_init_index_x][rho_init_index_y][rho_init_index_z] * (u(IDN, k, j, i) / 29.6) * std::max(0.0,damping); + u(IDN, k, j, i) += perturb_amplitude * rho_init[rho_init_index_x][rho_init_index_y][rho_init_index_z]; + + } }, passive_scalar); @@ -825,10 +944,11 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { PARTHENON_REQUIRE( u(IM1, k, j, i) == 0.0 && u(IM2, k, j, i) == 0.0 && u(IM3, k, j, i) == 0.0, "Found existing non-zero velocity when setting velocity perturbations."); - + u(IM1, k, j, i) = rho * perturb_pack(b, 0, k, j, i); u(IM2, k, j, i) = rho * perturb_pack(b, 1, k, j, i); u(IM3, k, j, i) = rho * perturb_pack(b, 2, k, j, i); + // No need to touch the energy yet as we'll normalize later lsum += (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) * @@ -887,29 +1007,32 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { auto perturb_pack = md->PackVariables(std::vector{"tmp_perturb"}); // Modifying the magnetic potential so that it follows the rho profile + + /* pmb->par_for( "Init sigma_b", 0, num_blocks - 1, kb.s-1, kb.e+1, jb.s-1, jb.e+1, ib.s-1, ib.e+1, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { 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 r = sqrt(SQR(x) + SQR(y) + SQR(z)); auto pmb = md->GetBlockData(b)->GetBlockPointer(); 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; - + perturb_pack(b, 0, k, j, i) *= 1 / (1 + (r / r_scale)); perturb_pack(b, 1, k, j, i) *= 1 / (1 + (r / r_scale)); perturb_pack(b, 2, k, j, i) *= 1 / (1 + (r / r_scale)); }); + */ pmb->par_reduce( "Init sigma_b", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -947,12 +1070,12 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &b2_sum, 1, MPI_PARTHENON_REAL, MPI_SUM, MPI_COMM_WORLD)); #endif // MPI_PARALLEL - + const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min; const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min; const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min; auto b_norm = std::sqrt(b2_sum / (Lx * Ly * Lz) / (SQR(sigma_b))); - + pmb->par_for( "Norm sigma_b", 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) { diff --git a/src/pgen/cluster/cluster_gravity.hpp b/src/pgen/cluster/cluster_gravity.hpp index 43cd4c69..dd7c20a0 100644 --- a/src/pgen/cluster/cluster_gravity.hpp +++ b/src/pgen/cluster/cluster_gravity.hpp @@ -17,7 +17,7 @@ namespace cluster { // Types of BCG's -enum class BCG {NONE, HERNQUIST, ISOTHERMAL}; +enum class BCG {NONE, HERNQUIST, VAUCOULEURS, ISOTHERMAL}; // Hernquiest BCG: Hernquist 1990 DOI:10.1086/168845 /************************************************************ @@ -132,6 +132,8 @@ class ClusterGravity { which_bcg_g_ = BCG::HERNQUIST; } else if (which_bcg_g_str == "ISOTHERMAL") { which_bcg_g_ = BCG::ISOTHERMAL; + } else if (which_bcg_g_str == "VAUCOULEURS") { + which_bcg_g_ = BCG::VAUCOULEURS; } else { @@ -172,6 +174,12 @@ class ClusterGravity { g_const_bcg_ = 2 * units.k_boltzmann() * T_bcg_s_ / (mu * units.mh()); } + + if (which_bcg_g_str == "VAUCOULEURS") { + + g_const_bcg_ = 1; + + } else { @@ -191,16 +199,16 @@ class ClusterGravity { : ClusterGravity(pin) { hydro_pkg->AddParam<>("cluster_gravity", *this); } - + // Inline functions to compute gravitational acceleration KOKKOS_INLINE_FUNCTION parthenon::Real g_from_r(const parthenon::Real r_in) const __attribute__((always_inline)) { - + const parthenon::Real r = std::max(r_in, smoothing_r_); const parthenon::Real r2 = r * r; - + parthenon::Real g_r = 0; - + // Add NFW gravity if (include_nfw_g_) { g_r += g_const_nfw_ * (log(1 + r / r_nfw_s_) - r / (r + r_nfw_s_)) / r2; @@ -214,6 +222,16 @@ class ClusterGravity { g_r += g_const_bcg_ / ((1 + r / r_bcg_s_) * (1 + r / r_bcg_s_)); case BCG::ISOTHERMAL: g_r += g_const_bcg_ / r; + case BCG::VAUCOULEURS: + + // From Matthews et al. 2005 + parthenon::Real vaucouleurs_K1 = std::pow(1e3 * r, 0.5975) / (3.206e-7); + parthenon::Real vaucouleurs_K2 = std::pow(1e3 * r, 1.849) / (1.861e-6); + parthenon::Real vaucouleurs_K = std::pow(std::pow(vaucouleurs_K1, 0.9) + std::pow(vaucouleurs_K2, 0.9),-1/0.9); + + vaucouleurs_K *= (3.15576e+16 * 3.15576e+16 / 3.085677580962325e+24); + g_r += 0.8 * vaucouleurs_K; + break; } @@ -231,7 +249,7 @@ class ClusterGravity { const parthenon::Real r = std::max(r_in, smoothing_r_); parthenon::Real rho = 0; - + // Add NFW gravity if (include_nfw_g_) { rho += rho_const_nfw_ / (r * pow(r + r_nfw_s_, 2)); diff --git a/src/pgen/cluster/cluster_gravity_working.hpp b/src/pgen/cluster/cluster_gravity_working.hpp new file mode 100644 index 00000000..ef6b6fb3 --- /dev/null +++ b/src/pgen/cluster/cluster_gravity_working.hpp @@ -0,0 +1,265 @@ +#ifndef CLUSTER_CLUSTER_GRAVITY_HPP_ +#define CLUSTER_CLUSTER_GRAVITY_HPP_ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file cluster_gravity.hpp +// \brief Class for defining gravitational acceleration for a cluster+bcg+smbh + +// Parthenon headers +#include + +// AthenaPK headers +#include "../../units.hpp" + +namespace cluster { + +// Types of BCG's +enum class BCG {NONE, HERNQUIST, ISOTHERMAL}; +// Hernquiest BCG: Hernquist 1990 DOI:10.1086/168845 + +/************************************************************ + * Cluster Gravity Class, for computing gravitational acceleration + * Lightweight object for inlined computation within kernels + ************************************************************/ +class ClusterGravity { + + // Parameters for which gravity sources to include + bool include_nfw_g_; + BCG which_bcg_g_; + bool include_smbh_g_; + + // NFW Parameters + parthenon::Real r_nfw_s_; + // G , Mass, and Constants rolled into one + parthenon::Real gravitationnal_const_; + parthenon::Real g_const_nfw_; + parthenon::Real rho_const_nfw_; + + // BCG Parameters + parthenon::Real alpha_bcg_s_; + parthenon::Real beta_bcg_s_; + parthenon::Real r_bcg_s_; + parthenon::Real T_bcg_s_; + // G , Mass, and Constants rolled into one + parthenon::Real g_const_bcg_; + parthenon::Real rho_const_bcg_; + + // SMBH Parameters + // G , Mass, and Constants rolled into one + parthenon::Real g_const_smbh_; + + // Radius underwhich to truncate + parthenon::Real smoothing_r_; + + // Static Helper functions to calculate constants to minimize in-kernel work + // ie. compute the constants characterizing the profile (M_nfw, etc...) + + static parthenon::Real calc_R_nfw_s(const parthenon::Real rho_crit, + const parthenon::Real m_nfw_200, + const parthenon::Real c_nfw) { + const parthenon::Real rho_nfw_0 = + 200 / 3. * rho_crit * pow(c_nfw, 3.) / (log(1 + c_nfw) - c_nfw / (1 + c_nfw)); + const parthenon::Real R_nfw_s = + pow(m_nfw_200 / (4 * M_PI * rho_nfw_0 * (log(1 + c_nfw) - c_nfw / (1 + c_nfw))), + 1. / 3.); + return R_nfw_s; + } + static parthenon::Real calc_g_const_nfw(const parthenon::Real gravitational_constant, + const parthenon::Real m_nfw_200, + const parthenon::Real c_nfw) { + return gravitational_constant * m_nfw_200 / (log(1 + c_nfw) - c_nfw / (1 + c_nfw)); + } + static parthenon::Real calc_rho_const_nfw(const parthenon::Real gravitational_constant, + const parthenon::Real m_nfw_200, + const parthenon::Real c_nfw) { + return m_nfw_200 / (4 * M_PI * (log(1 + c_nfw) - c_nfw / (1 + c_nfw))); + } + static parthenon::Real calc_g_const_bcg(const parthenon::Real gravitational_constant, + BCG which_bcg_g, const parthenon::Real m_bcg_s, + const parthenon::Real r_bcg_s, + const parthenon::Real alpha_bcg_s, + const parthenon::Real beta_bcg_s) { + switch (which_bcg_g) { + case BCG::NONE: + return 0; + case BCG::HERNQUIST: + return gravitational_constant * m_bcg_s / (r_bcg_s * r_bcg_s); + } + return NAN; + } + static parthenon::Real calc_rho_const_bcg(const parthenon::Real gravitational_constant, + BCG which_bcg_g, + const parthenon::Real m_bcg_s, + const parthenon::Real r_bcg_s, + const parthenon::Real alpha_bcg_s, + const parthenon::Real beta_bcg_s) { + switch (which_bcg_g) { + case BCG::NONE: + return 0; + case BCG::HERNQUIST: + return m_bcg_s * r_bcg_s / (2 * M_PI); + } + return NAN; + } + static KOKKOS_INLINE_FUNCTION parthenon::Real + calc_g_const_smbh(const parthenon::Real gravitational_constant, + const parthenon::Real m_smbh) { + return gravitational_constant * m_smbh; + } + + public: + // ClusterGravity(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg) + // is called from cluster.cpp to add the ClusterGravity object to hydro_pkg + // + // ClusterGravity(parthenon::ParameterInput *pin) is used in SNIAFeedback to + // calculate the BCG density profile + ClusterGravity(parthenon::ParameterInput *pin) { + Units units(pin); + + gravitationnal_const_ = units.gravitational_constant(); + + // Determine which element to include + include_nfw_g_ = + pin->GetOrAddBoolean("problem/cluster/gravity", "include_nfw_g", false); + const std::string which_bcg_g_str = + pin->GetOrAddString("problem/cluster/gravity", "which_bcg_g", "NONE"); + if (which_bcg_g_str == "NONE") { + which_bcg_g_ = BCG::NONE; + } else if (which_bcg_g_str == "HERNQUIST") { + which_bcg_g_ = BCG::HERNQUIST; + } else if (which_bcg_g_str == "ISOTHERMAL") { + which_bcg_g_ = BCG::ISOTHERMAL; + } + + else { + std::stringstream msg; + msg << "### FATAL ERROR in function [InitUserMeshData]" << std::endl + << "Unknown BCG type " << which_bcg_g_str << std::endl; + PARTHENON_FAIL(msg); + } + + include_smbh_g_ = + pin->GetOrAddBoolean("problem/cluster/gravity", "include_smbh_g", false); + + // Initialize the NFW Profile + const parthenon::Real hubble_parameter = pin->GetOrAddReal( + "problem/cluster", "hubble_parameter", 70 * units.km_s() / units.mpc()); + const parthenon::Real rho_crit = 3 * hubble_parameter * hubble_parameter / + (8 * M_PI * units.gravitational_constant()); + + const auto He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction"); + const auto mu = 1 / (He_mass_fraction * 3. / 4. + (1 - He_mass_fraction) * 2); + const parthenon::Real M_nfw_200 = + pin->GetOrAddReal("problem/cluster/gravity", "m_nfw_200", 8.5e14 * units.msun()); + const parthenon::Real c_nfw = + pin->GetOrAddReal("problem/cluster/gravity", "c_nfw", 6.81); + r_nfw_s_ = calc_R_nfw_s(rho_crit, M_nfw_200, c_nfw); + g_const_nfw_ = calc_g_const_nfw(units.gravitational_constant(), M_nfw_200, c_nfw); + + // Initialize the BCG Profile + alpha_bcg_s_ = pin->GetOrAddReal("problem/cluster/gravity", "alpha_bcg_s", 0.1); + beta_bcg_s_ = pin->GetOrAddReal("problem/cluster/gravity", "beta_bcg_s", 1.43); + const parthenon::Real M_bcg_s = + pin->GetOrAddReal("problem/cluster/gravity", "m_bcg_s", 7.5e10 * units.msun()); + r_bcg_s_ = pin->GetOrAddReal("problem/cluster/gravity", "r_bcg_s", 4 * units.kpc()); + T_bcg_s_ = pin->GetOrAddReal("problem/cluster/gravity", "T_bcg_s", 10000); // Temperature in the isothermal case + + if (which_bcg_g_str == "ISOTHERMAL") { + + g_const_bcg_ = 2 * units.k_boltzmann() * T_bcg_s_ / (mu * units.mh()); + + } + + else { + + g_const_bcg_ = calc_g_const_bcg(units.gravitational_constant(), which_bcg_g_, M_bcg_s, + r_bcg_s_, alpha_bcg_s_, beta_bcg_s_); + } + + const parthenon::Real m_smbh = + pin->GetOrAddReal("problem/cluster/gravity", "m_smbh", 3.4e8 * units.msun()); + g_const_smbh_ = calc_g_const_smbh(units.gravitational_constant(), m_smbh), + + smoothing_r_ = + pin->GetOrAddReal("problem/cluster/gravity", "g_smoothing_radius", 0.0); + } + + ClusterGravity(parthenon::ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg) + : ClusterGravity(pin) { + hydro_pkg->AddParam<>("cluster_gravity", *this); + } + + // Inline functions to compute gravitational acceleration + KOKKOS_INLINE_FUNCTION parthenon::Real g_from_r(const parthenon::Real r_in) const + __attribute__((always_inline)) { + + const parthenon::Real r = std::max(r_in, smoothing_r_); + const parthenon::Real r2 = r * r; + + parthenon::Real g_r = 0; + + // Add NFW gravity + if (include_nfw_g_) { + g_r += g_const_nfw_ * (log(1 + r / r_nfw_s_) - r / (r + r_nfw_s_)) / r2; + } + + // Add BCG gravity + switch (which_bcg_g_) { + case BCG::NONE: + break; + case BCG::HERNQUIST: + g_r += g_const_bcg_ / ((1 + r / r_bcg_s_) * (1 + r / r_bcg_s_)); + case BCG::ISOTHERMAL: + g_r += g_const_bcg_ / r; + break; + } + + // Add SMBH, point mass gravity + if (include_smbh_g_) { + g_r += g_const_smbh_ / r2; + } + + return g_r; + } + // Inline functions to compute density + KOKKOS_INLINE_FUNCTION parthenon::Real rho_from_r(const parthenon::Real r_in) const + __attribute__((always_inline)) { + + const parthenon::Real r = std::max(r_in, smoothing_r_); + + parthenon::Real rho = 0; + + // Add NFW gravity + if (include_nfw_g_) { + rho += rho_const_nfw_ / (r * pow(r + r_nfw_s_, 2)); + } + + // Add BCG gravity + switch (which_bcg_g_) { + case BCG::NONE: + break; + case BCG::HERNQUIST: + rho += rho_const_bcg_ / (r * pow(r + r_bcg_s_, 3)); + break; + case BCG::ISOTHERMAL: + rho += g_const_bcg_ * 1 / (4 * M_PI * gravitationnal_const_ * r * r); + } + + // SMBH, point mass gravity -- density is not defined. Throw an error + if (include_smbh_g_ && r <= smoothing_r_) { + Kokkos::abort("ClusterGravity::SMBH density is not defined"); + } + + return rho; + } + + // SNIAFeedback needs to be a friend to disable the SMBH and NFW + friend class SNIAFeedback; +}; + +} // namespace cluster + +#endif // CLUSTER_CLUSTER_GRAVITY_HPP_ diff --git a/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp b/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp index 933182dc..3d82258b 100644 --- a/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp +++ b/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp @@ -77,7 +77,7 @@ class HydrostaticEquilibriumSphere { const parthenon::Real n = rho / (mu_ * mh_); return n; } - + // Get electron number density from density KOKKOS_INLINE_FUNCTION parthenon::Real ne_from_rho(const parthenon::Real rho) const { const parthenon::Real ne = mu_ / mu_e_ * n_from_rho(rho); diff --git a/src/refinement/other.cpp b/src/refinement/other.cpp index 586521f9..fa1ffc94 100644 --- a/src/refinement/other.cpp +++ b/src/refinement/other.cpp @@ -16,17 +16,19 @@ using parthenon::IndexRange; // refinement condition: check max density parthenon::AmrTag MaxDensity(MeshBlockData *rc) { - auto pmb = rc->GetBlockPointer(); - auto w = rc->Get("prim").data; + auto pmb = rc->GetBlockPointer(); + auto w = rc->Get("prim").data; + //auto &coords = pmb->coords; + const auto deref_below = pmb->packages.Get("Hydro")->Param("refinement/maxdensity_deref_below"); const auto refine_above = pmb->packages.Get("Hydro")->Param("refinement/maxdensity_refine_above"); - + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); - + Real maxrho = 0.0; pmb->par_reduce( "overdens check refinement", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e + 1, @@ -34,9 +36,12 @@ parthenon::AmrTag MaxDensity(MeshBlockData *rc) { lmaxrho = std::max(lmaxrho, w(IDN, k, j, i)); }, Kokkos::Max(maxrho)); - + if (maxrho > refine_above) return parthenon::AmrTag::refine; if (maxrho < deref_below) return parthenon::AmrTag::derefine; + + + return parthenon::AmrTag::same; }