diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 076fbb9f..06ee3df4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,6 +18,7 @@ add_executable( refinement/gradient.cpp refinement/other.cpp utils/few_modes_ft.cpp + utils/few_modes_ft_lognormal.cpp ) add_subdirectory(pgen) diff --git a/src/hydro/srcterms/tabular_cooling.hpp b/src/hydro/srcterms/tabular_cooling.hpp index bea2499b..0e9f9c24 100644 --- a/src/hydro/srcterms/tabular_cooling.hpp +++ b/src/hydro/srcterms/tabular_cooling.hpp @@ -155,19 +155,19 @@ class CoolingTableObj { log_lambda = 0.5 * log_temp - 0.5 * log_temp_final_ + log_lambdas_(n_temp_ - 1); } else { // Inside table, interpolate assuming log spaced temperatures - + // Determine where temp is in the table const unsigned int i_temp = static_cast((log_temp - log_temp_start_) / d_log_temp_); const Real log_temp_i = log_temp_start_ + d_log_temp_ * i_temp; - + // log_temp should be between log_temps[i_temp] and log_temps[i_temp+1] PARTHENON_REQUIRE(log_temp >= log_temp_i && log_temp <= log_temp_i + d_log_temp_, "FATAL ERROR in [CoolingTable::DeDt]: Failed to find log_temp"); - + const Real log_lambda_i = log_lambdas_(i_temp); const Real log_lambda_ip1 = log_lambdas_(i_temp + 1); - + // Linearly interpolate lambda at log_temp log_lambda = log_lambda_i + (log_temp - log_temp_i) * (log_lambda_ip1 - log_lambda_i) / d_log_temp_; @@ -181,7 +181,7 @@ class CoolingTableObj { return de_dt; } - + KOKKOS_INLINE_FUNCTION parthenon::Real DeDt(const parthenon::Real &e, const parthenon::Real &rho) const { bool is_valid = true; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index e2b4e7ca..38c149ea 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -23,6 +23,7 @@ #include // stringstream #include // runtime_error #include // c_str() + // #include // HDF5 #include "../../external/HighFive/include/highfive/H5Easy.hpp" @@ -41,6 +42,7 @@ #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" @@ -57,6 +59,7 @@ 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; template void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, @@ -282,8 +285,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddParam<>("dipole_b_field_my", dipole_b_field_my); hydro_pkg->AddParam<>("dipole_b_field_mz", dipole_b_field_mz); } - - + /************************************************************ * Read Cluster Gravity Parameters ************************************************************/ @@ -310,8 +312,9 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd HydrostaticEquilibriumSphere hse_sphere(pin, hydro_pkg, cluster_gravity, entropy_profile); - - /************************************************************ + + + /************************************************************ * Read Precessing Jet Coordinate system ************************************************************/ @@ -326,6 +329,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd /************************************************************ * Read AGN Triggering ************************************************************/ + AGNTriggering agn_triggering(pin, hydro_pkg); /************************************************************ @@ -334,6 +338,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd // Build Magnetic Tower MagneticTower magnetic_tower(pin, hydro_pkg); + // Determine if magnetic_tower_power_scaling is needed // Is AGN Power and Magnetic fraction non-zero? @@ -398,7 +403,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddField("mach_sonic", m); // temperature hydro_pkg->AddField("temperature", m); - + if (hydro_pkg->Param("enable_cooling") == Cooling::tabular) { // cooling time hydro_pkg->AddField("cooling_time", m); @@ -411,25 +416,47 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd // plasma beta hydro_pkg->AddField("plasma_beta", m); } + + /************************************************************ + * Read Velocity 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_rho", m); - /************************************************************ - * Read Density perturbation - ************************************************************/ - - //const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); - //hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); - + } /************************************************************ * Read Velocity perturbation ************************************************************/ - + + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); if (sigma_v != 0.0) { // peak of init vel perturb @@ -455,6 +482,11 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd std::vector({3})); hydro_pkg->AddField("tmp_perturb", m); } + + /************************************************************ + * Read Magnetic field perturbation + ************************************************************/ + const auto sigma_b = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_b", 0.0); if (sigma_b != 0.0) { PARTHENON_REQUIRE_THROWS(hydro_pkg->Param("fluid") == Fluid::glmmhd, @@ -487,6 +519,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddField("tmp_perturb", m); } } + } //======================================================================================== @@ -533,8 +566,6 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real My = rho * uy; const Real Mz = rho * uz; const Real E = rho * (0.5 * (ux * ux + uy * uy + uz * uz) + pres / (gm1 * rho)); - - std::cout << "Setting uniform gas"; parthenon::par_for( DEFAULT_LOOP_PATTERN, "Cluster::ProblemGenerator::UniformGas", @@ -581,7 +612,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { u(IEN, k, j, i) = P_r / gm1; }); } - + if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { /************************************************************ * Initialize the initial magnetic field state via a vector potential @@ -687,8 +718,8 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { } // END if(hydro_pkg->Param("fluid") == Fluid::glmmhd) } - - /************************************************************ + + /************************************************************ * Initial parameters ************************************************************/ @@ -702,22 +733,67 @@ 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 + + /************************************************************ + * Set initial density perturbations (read from HDF5 file) ************************************************************/ const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + const bool full_box = pin->GetOrAddBoolean("problem/cluster/init_perturb", "full_box", true); + 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 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) { + if (spherical_collapse == 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) { @@ -727,7 +803,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { 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); @@ -747,12 +823,14 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { u(IDN, k, j, i) = foreground_density; } + //u(IDN, k, j, i) += foreground_density * std::exp(- 2 * SQR(r) / SQR(overdensity_radius) ) ; + }, passive_scalar); - + } - + /* -------------- Setting up a clumpy atmosphere -------------- 1) Extract the values of the density from an input hdf5 file using H5Easy @@ -764,23 +842,19 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { if (init_perturb_rho == true) { - auto init_perturb_rho_file = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_file"); - auto init_perturb_rho_keys = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_keys"); - - hydro_pkg->AddParam<>("cluster/init_perturb_rho_file", init_perturb_rho_file); - hydro_pkg->AddParam<>("cluster/init_perturb_rho_keys", init_perturb_rho_keys); + auto filename_rho = pin->GetOrAddString("problem/cluster/init_perturb", "init_perturb_rho_file","none"); + auto grid_size = pin->GetOrAddInteger("parthenon/mesh","nx1",256); - std::cout << "Setting density perturbation"; - - // Read HDF5 file containing the density - std::string filename_rho = "/work/bbd0833/test/rho.h5"; std::string keys_rho = "data"; H5Easy::File file(filename_rho, HighFive::File::ReadOnly); - auto rho_init = H5Easy::load, 64>, 64>>(file, keys_rho); + auto rho_init = H5Easy::load, 256>, 256>>(file, keys_rho); + + //std::vector rho_init(grid_size * grid_size * grid_size); + //std::vector rho_init(grid_size * grid_size * grid_size) = H5Easy::load(file, keys_rho); Real passive_scalar = 0.0; // Useless - std::cout << "entering initialisation of rho field"; + std::cout << "Entering initialisation of rho field"; pmb->par_reduce( "Init density field", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -792,14 +866,26 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const auto gjs = pmbb->loc.lx2 * pmb->block_size.nx2; const auto gks = pmbb->loc.lx3 * pmb->block_size.nx3; - //std::cout << "gis=" << gis << "\n"; - const auto &coords = cons.GetCoords(b); const auto &u = cons(b); // Adding the value of the density field + if (full_box){ + + u(IDN, k, j, i) += rho_init[gks + k - 2][gjs + j - 2][gis + i - 2]; + + } + + else { + + const Real z = coords.Xc<3>(k); + + // if (z > -thickness_ism / 2 && z < thickness_ism / 2){ + + u(IDN, k, j, i) += rho_init[gks + k - 2][gjs + j - 2][gis + i - 2] * std::exp(- SQR(z) / SQR(thickness_ism)) ; + - u(IDN, k, j, i) = rho_init[gks + k - 2][gjs + j - 2][gis + i - 2]; + } // Ring of overpressure // compute radius @@ -831,7 +917,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { u(IEN, k, j, i) += 0.5 * u(IDN, k, j, i) * (pow(velocity_blast,2)); } } - + if (spherical_collapse){ const Real x = coords.Xc<1>(i); @@ -841,12 +927,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); u(IDN, k, j, i) = 3; - - if (r > 0 && r < overdensity_radius){ - - u(IDN, k, j, i) *= 50; - - } + u(IDN, k, j, i) *= 50 * std::exp(- 2 * SQR(r) / SQR(overdensity_radius)); } @@ -854,6 +935,57 @@ 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); + + 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++) { + 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"); + + //Real v2_sum_rho = 0.0; // used for normalization + + //auto perturb_pack_rho = md->PackVariables(std::vector{"tmp_perturb_rho"}); + + //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 + + } /************************************************************ * Set initial velocity perturbations (requires no other velocities for now) @@ -925,7 +1057,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { u(IDN, k, j, i); }); } - + /************************************************************ * Set initial magnetic field perturbations (resets magnetic field field) ************************************************************/ diff --git a/src/pgen/cluster_working.cpp b/src/pgen/cluster_working.cpp new file mode 100644 index 00000000..cbdb6dd3 --- /dev/null +++ b/src/pgen/cluster_working.cpp @@ -0,0 +1,1125 @@ +//======================================================================================== +// 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.cpp +// \brief Idealized galaxy cluster problem generator +// +// Setups up an idealized galaxy cluster with an ACCEPT-like entropy profile in +// hydrostatic equilbrium with an NFW+BCG+SMBH gravitational profile, +// optionally with an initial magnetic tower field. Includes AGN feedback, AGN +// triggering via cold gas, simple SNIA Feedback(TODO) +//======================================================================================== + +// C headers + +// C++ headers +#include // min, max +#include // sqrt() +#include // fopen(), fprintf(), freopen() +#include // endl +#include +#include // stringstream +#include // runtime_error +#include // c_str() + +// #include // HDF5 +#include "../../external/HighFive/include/highfive/H5Easy.hpp" + +// Parthenon headers +#include "kokkos_abstraction.hpp" +#include "mesh/domain.hpp" +#include "mesh/mesh.hpp" +#include +#include + +// AthenaPK headers +#include "../eos/adiabatic_glmmhd.hpp" +#include "../eos/adiabatic_hydro.hpp" +#include "../hydro/hydro.hpp" +#include "../hydro/srcterms/gravitational_field.hpp" +#include "../hydro/srcterms/tabular_cooling.hpp" +#include "../main.hpp" +#include "../utils/few_modes_ft.hpp" + +// Cluster headers +#include "cluster/agn_feedback.hpp" +#include "cluster/agn_triggering.hpp" +#include "cluster/cluster_gravity.hpp" +#include "cluster/entropy_profiles.hpp" +#include "cluster/hydrostatic_equilibrium_sphere.hpp" +#include "cluster/magnetic_tower.hpp" +#include "cluster/snia_feedback.hpp" +#include "parthenon_array_generic.hpp" +#include "utils/error_checking.hpp" + +namespace cluster { +using namespace parthenon::driver::prelude; +using namespace parthenon::package::prelude; +using utils::few_modes_ft::FewModesFT; + +template +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt, const EOS eos) { + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + // Apply clips -- ceilings on temperature, velocity, alfven velocity, and + // density floor -- within a radius of the AGN + const auto &dfloor = hydro_pkg->Param("cluster_dfloor"); + const auto &eceil = hydro_pkg->Param("cluster_eceil"); + const auto &vceil = hydro_pkg->Param("cluster_vceil"); + const auto &vAceil = hydro_pkg->Param("cluster_vAceil"); + const auto &clip_r = hydro_pkg->Param("cluster_clip_r"); + + if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits::infinity() || + vceil < std::numeric_limits::infinity() || + vAceil < std::numeric_limits::infinity())) { + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + const auto nhydro = hydro_pkg->Param("nhydro"); + const auto nscalars = hydro_pkg->Param("nscalars"); + + const Real clip_r2 = SQR(clip_r); + const Real vceil2 = SQR(vceil); + const Real vAceil2 = SQR(vAceil); + const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Cluster::ApplyClusterClips", 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) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + const auto &coords = cons_pack.GetCoords(b); + + const Real r2 = + SQR(coords.Xc<1>(i)) + SQR(coords.Xc<2>(j)) + SQR(coords.Xc<3>(k)); + + if (r2 < clip_r2) { + // Cell falls within clipping radius + + const int gks = 0; + const int gjs = 0; + const int gis = 0; + + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs, gis); // Three last parameters are just passive here + + if (dfloor > 0) { + const Real rho = prim(IDN, k, j, i); + if (rho < dfloor) { + cons(IDN, k, j, i) = dfloor; + prim(IDN, k, j, i) = dfloor; + } + } + + if (vceil < std::numeric_limits::infinity()) { + // Apply velocity ceiling + const Real v2 = SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) + + SQR(prim(IV3, k, j, i)); + if (v2 > vceil2) { + // Fix the velocity to the velocity ceiling + const Real v = sqrt(v2); + cons(IM1, k, j, i) *= vceil / v; + cons(IM2, k, j, i) *= vceil / v; + cons(IM3, k, j, i) *= vceil / v; + prim(IV1, k, j, i) *= vceil / v; + prim(IV2, k, j, i) *= vceil / v; + prim(IV3, k, j, i) *= vceil / v; + + // Remove kinetic energy + cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); + } + } + + if (vAceil2 < std::numeric_limits::infinity()) { + // Apply Alfven velocity ceiling by raising density + const Real rho = prim(IDN, k, j, i); + const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) + + SQR(prim(IB3, k, j, i))); + + // compute Alfven mach number + const Real va2 = (B2 / rho); + + if (va2 > vAceil2) { + // Increase the density to match the alfven velocity ceiling + const Real rho_new = std::sqrt(B2 / vAceil2); + cons(IDN, k, j, i) = rho_new; + prim(IDN, k, j, i) = rho_new; + } + } + + if (eceil < std::numeric_limits::infinity()) { + // Apply internal energy ceiling as a pressure ceiling + const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i)); + if (internal_e > eceil) { + cons(IEN, k, j, i) -= prim(IDN, k, j, i) * (internal_e - eceil); + prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil; + } + } + } + }); + } +} + +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + auto fluid = hydro_pkg->Param("fluid"); + if (fluid == Fluid::euler) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else if (fluid == Fluid::glmmhd) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else { + PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS"); + } +} + +void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const bool &gravity_srcterm = hydro_pkg->Param("gravity_srcterm"); + + if (gravity_srcterm) { + const ClusterGravity &cluster_gravity = + hydro_pkg->Param("cluster_gravity"); + + GravitationalFieldSrcTerm(md, beta_dt, cluster_gravity); + } + + const auto &agn_feedback = hydro_pkg->Param("agn_feedback"); + agn_feedback.FeedbackSrcTerm(md, beta_dt, tm); + + const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); + magnetic_tower.FixedFieldSrcTerm(md, beta_dt, tm); + + const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); + snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + + ApplyClusterClips(md, tm, beta_dt); +}; + +Real ClusterEstimateTimestep(MeshData *md) { + Real min_dt = std::numeric_limits::max(); + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + // TODO time constraints imposed by thermal AGN feedback, jet velocity, + // magnetic tower + const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + const Real agn_triggering_min_dt = agn_triggering.EstimateTimeStep(md); + min_dt = std::min(min_dt, agn_triggering_min_dt); + + return min_dt; +} + +//======================================================================================== +//! \fn void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor +//! *hydro_pkg) \brief Init package data from parameter input +//======================================================================================== + +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); + hydro_pkg->AddParam<>("uniform_gas_uz", uniform_gas_uz); + hydro_pkg->AddParam<>("uniform_gas_pres", uniform_gas_pres); + } + + /************************************************************ + * Read Uniform Magnetic Field + ************************************************************/ + + const bool init_uniform_b_field = pin->GetOrAddBoolean( + "problem/cluster/uniform_b_field", "init_uniform_b_field", false); + hydro_pkg->AddParam<>("init_uniform_b_field", init_uniform_b_field); + + if (init_uniform_b_field) { + const Real uniform_b_field_bx = pin->GetReal("problem/cluster/uniform_b_field", "bx"); + const Real uniform_b_field_by = pin->GetReal("problem/cluster/uniform_b_field", "by"); + const Real uniform_b_field_bz = pin->GetReal("problem/cluster/uniform_b_field", "bz"); + + hydro_pkg->AddParam<>("uniform_b_field_bx", uniform_b_field_bx); + hydro_pkg->AddParam<>("uniform_b_field_by", uniform_b_field_by); + hydro_pkg->AddParam<>("uniform_b_field_bz", uniform_b_field_bz); + } + + /************************************************************ + * Read Uniform Magnetic Field + ************************************************************/ + + const bool init_dipole_b_field = pin->GetOrAddBoolean("problem/cluster/dipole_b_field", + "init_dipole_b_field", false); + hydro_pkg->AddParam<>("init_dipole_b_field", init_dipole_b_field); + + if (init_dipole_b_field) { + const Real dipole_b_field_mx = pin->GetReal("problem/cluster/dipole_b_field", "mx"); + const Real dipole_b_field_my = pin->GetReal("problem/cluster/dipole_b_field", "my"); + const Real dipole_b_field_mz = pin->GetReal("problem/cluster/dipole_b_field", "mz"); + + hydro_pkg->AddParam<>("dipole_b_field_mx", dipole_b_field_mx); + hydro_pkg->AddParam<>("dipole_b_field_my", dipole_b_field_my); + hydro_pkg->AddParam<>("dipole_b_field_mz", dipole_b_field_mz); + } + + /************************************************************ + * Read Cluster Gravity Parameters + ************************************************************/ + + // Build cluster_gravity object + ClusterGravity cluster_gravity(pin, hydro_pkg); + // hydro_pkg->AddParam<>("cluster_gravity", cluster_gravity); + + // Include gravity as a source term during evolution + const bool gravity_srcterm = + pin->GetBoolean("problem/cluster/gravity", "gravity_srcterm"); + hydro_pkg->AddParam<>("gravity_srcterm", gravity_srcterm); + + /************************************************************ + * Read Initial Entropy Profile + ************************************************************/ + + // Build entropy_profile object + ACCEPTEntropyProfile entropy_profile(pin); + + /************************************************************ + * Build Hydrostatic Equilibrium Sphere + ************************************************************/ + + HydrostaticEquilibriumSphere hse_sphere(pin, hydro_pkg, cluster_gravity, + entropy_profile); + + + /************************************************************ + * Read Precessing Jet Coordinate system + ************************************************************/ + + JetCoordsFactory jet_coords_factory(pin, hydro_pkg); + + /************************************************************ + * Read AGN Feedback + ************************************************************/ + + AGNFeedback agn_feedback(pin, hydro_pkg); + + /************************************************************ + * Read AGN Triggering + ************************************************************/ + + AGNTriggering agn_triggering(pin, hydro_pkg); + + /************************************************************ + * Read Magnetic Tower + ************************************************************/ + + // Build Magnetic Tower + MagneticTower magnetic_tower(pin, hydro_pkg); + + + // Determine if magnetic_tower_power_scaling is needed + // Is AGN Power and Magnetic fraction non-zero? + bool magnetic_tower_power_scaling = + (agn_feedback.magnetic_fraction_ != 0 && + (agn_feedback.fixed_power_ != 0 || + agn_triggering.triggering_mode_ != AGNTriggeringMode::NONE)); + hydro_pkg->AddParam("magnetic_tower_power_scaling", magnetic_tower_power_scaling); + + /************************************************************ + * Read SNIA Feedback + ************************************************************/ + + SNIAFeedback snia_feedback(pin, hydro_pkg); + + /************************************************************ + * Read Clips (ceilings and floors) + ************************************************************/ + + // Disable all clips by default with a negative radius clip + Real clip_r = pin->GetOrAddReal("problem/cluster/clips", "clip_r", -1.0); + + // By default disable floors by setting a negative value + Real dfloor = pin->GetOrAddReal("problem/cluster/clips", "dfloor", -1.0); + + // By default disable ceilings by setting to infinity + Real vceil = pin->GetOrAddReal("problem/cluster/clips", "vceil", + std::numeric_limits::infinity()); + Real vAceil = pin->GetOrAddReal("problem/cluster/clips", "vAceil", + std::numeric_limits::infinity()); + Real Tceil = pin->GetOrAddReal("problem/cluster/clips", "Tceil", + std::numeric_limits::infinity()); + Real eceil = Tceil; + if (eceil < std::numeric_limits::infinity()) { + if (!hydro_pkg->AllParams().hasKey("mbar_over_kb")) { + PARTHENON_FAIL("Temperature ceiling requires units and gas composition. " + "Either set a 'units' block and the 'hydro/He_mass_fraction' in " + "input file or use a pressure floor " + "(defined code units) instead."); + } + auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + eceil = Tceil / mbar_over_kb / (hydro_pkg->Param("AdiabaticIndex") - 1.0); + } + hydro_pkg->AddParam("cluster_dfloor", dfloor); + hydro_pkg->AddParam("cluster_eceil", eceil); + hydro_pkg->AddParam("cluster_vceil", vceil); + hydro_pkg->AddParam("cluster_vAceil", vAceil); + hydro_pkg->AddParam("cluster_clip_r", clip_r); + + /************************************************************ + * Add derived fields + * NOTE: these must be filled in UserWorkBeforeOutput + ************************************************************/ + + auto m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector({1})); + + // log10 of cell-centered radius + hydro_pkg->AddField("log10_cell_radius", m); + // entropy + hydro_pkg->AddField("entropy", m); + // sonic Mach number v/c_s + hydro_pkg->AddField("mach_sonic", m); + // temperature + hydro_pkg->AddField("temperature", m); + + if (hydro_pkg->Param("enable_cooling") == Cooling::tabular) { + // cooling time + hydro_pkg->AddField("cooling_time", m); + } + + if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { + // alfven Mach number v_A/c_s + hydro_pkg->AddField("mach_alfven", m); + + // plasma beta + hydro_pkg->AddField("plasma_beta", m); + } + + + + + + + + /************************************************************ + * Read Density perturbation + ************************************************************/ + + //const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + //hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + + + /************************************************************ + * Read Velocity perturbation + ************************************************************/ + + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); + if (sigma_v != 0.0) { + // peak of init vel perturb + auto k_peak_v = pin->GetReal("problem/cluster/init_perturb", "k_peak_v"); + auto num_modes_v = + pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_v", 40); + auto sol_weight_v = + pin->GetOrAddReal("problem/cluster/init_perturb", "sol_weight_v", 1.0); + uint32_t rseed_v = pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_v", 1); + // 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 auto t_corr = 1e-10; + + auto k_vec_v = utils::few_modes_ft::MakeRandomModes(num_modes_v, k_peak_v, rseed_v); + + auto few_modes_ft = FewModesFT(pin, hydro_pkg, "cluster_perturb_v", num_modes_v, + k_vec_v, k_peak_v, sol_weight_v, t_corr, rseed_v); + hydro_pkg->AddParam<>("cluster/few_modes_ft_v", few_modes_ft); + + // 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); + } + const auto sigma_b = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_b", 0.0); + if (sigma_b != 0.0) { + PARTHENON_REQUIRE_THROWS(hydro_pkg->Param("fluid") == Fluid::glmmhd, + "Requested initial magnetic field perturbations but not " + "solving the MHD equations.") + // peak of init vel perturb + auto k_peak_b = pin->GetReal("problem/cluster/init_perturb", "k_peak_b"); + auto num_modes_b = + pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_b", 40); + uint32_t rseed_b = pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_b", 2); + // In principle arbitrary because the inital A_hat is 0 and the A_hat_new will contain + // the perturbation (and is normalized in the following to get the desired sigma_b) + const auto t_corr = 1e-10; + // This field should by construction have no compressive modes, so we fix the number. + const auto sol_weight_b = 1.0; + + auto k_vec_b = utils::few_modes_ft::MakeRandomModes(num_modes_b, k_peak_b, rseed_b); + + const bool fill_ghosts = true; // as we fill a vector potential to calc B + auto few_modes_ft = + FewModesFT(pin, hydro_pkg, "cluster_perturb_b", num_modes_b, k_vec_b, k_peak_b, + sol_weight_b, t_corr, rseed_b, fill_ghosts); + hydro_pkg->AddParam<>("cluster/few_modes_ft_b", few_modes_ft); + + // Add field for initial perturation (must not need to be consistent but defining it + // this way is easier for now). Only add if not already done for the velocity. + if (sigma_v == 0.0) { + Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, + std::vector({3})); + hydro_pkg->AddField("tmp_perturb", m); + } + } + +} + +//======================================================================================== +//! \fn void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) +//! \brief Generate problem data for all blocks on rank +// +// Note, this requires that parthenon/mesh/pack_size=-1 during initialization so that +// reductions work +//======================================================================================== + +void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { + // This could be more optimized, but require a refactor of init routines being called. + // However, given that it's just called during initial setup, this should not be a + // performance concern. + for (int b = 0; b < md->NumBlocks(); b++) { + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + auto hydro_pkg = pmb->packages.Get("Hydro"); + + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + + // Initialize the conserved variables + auto &u = pmb->meshblock_data.Get()->Get("cons").data; + + auto &coords = pmb->coords; + + // Get Adiabatic Index + const Real gam = pin->GetReal("hydro", "gamma"); + const Real gm1 = (gam - 1.0); + + /************************************************************ + * Initialize the initial hydro state + ************************************************************/ + const auto &init_uniform_gas = hydro_pkg->Param("init_uniform_gas"); + if (init_uniform_gas) { + const Real rho = hydro_pkg->Param("uniform_gas_rho"); + const Real ux = hydro_pkg->Param("uniform_gas_ux"); + const Real uy = hydro_pkg->Param("uniform_gas_uy"); + const Real uz = hydro_pkg->Param("uniform_gas_uz"); + const Real pres = hydro_pkg->Param("uniform_gas_pres"); + + const Real Mx = rho * ux; + const Real My = rho * uy; + const Real Mz = rho * uz; + const Real E = rho * (0.5 * (ux * ux + uy * uy + uz * uz) + pres / (gm1 * rho)); + + 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) { + u(IDN, k, j, i) = rho; + u(IM1, k, j, i) = Mx; + u(IM2, k, j, i) = My; + u(IM3, k, j, i) = Mz; + u(IEN, k, j, i) = E; + }); + + // end if(init_uniform_gas) + } else { + /************************************************************ + * Initialize a HydrostaticEquilibriumSphere + ************************************************************/ + const auto &he_sphere = + hydro_pkg + ->Param>( + "hydrostatic_equilibirum_sphere"); + + const auto P_rho_profile = he_sphere.generate_P_rho_profile(ib, jb, kb, coords); + + // 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) { + // 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)); + + // Get pressure and density from generated profile + const Real P_r = P_rho_profile.P_from_r(r); + const Real rho_r = P_rho_profile.rho_from_r(r); + + // Fill conserved states, 0 initial velocity + 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; + }); + } + + if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { + /************************************************************ + * Initialize the initial magnetic field state via a vector potential + ************************************************************/ + parthenon::ParArray4D A("A", 3, pmb->cellbounds.ncellsk(IndexDomain::entire), + pmb->cellbounds.ncellsj(IndexDomain::entire), + pmb->cellbounds.ncellsi(IndexDomain::entire)); + + IndexRange a_ib = ib; + a_ib.s -= 1; + a_ib.e += 1; + IndexRange a_jb = jb; + a_jb.s -= 1; + a_jb.e += 1; + IndexRange a_kb = kb; + a_kb.s -= 1; + a_kb.e += 1; + + /************************************************************ + * Initialize an initial magnetic tower + ************************************************************/ + const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); + + magnetic_tower.AddInitialFieldToPotential(pmb.get(), a_kb, a_jb, a_ib, A); + + /************************************************************ + * Add dipole magnetic field to the magnetic potential + ************************************************************/ + const auto &init_dipole_b_field = hydro_pkg->Param("init_dipole_b_field"); + if (init_dipole_b_field) { + const Real mx = hydro_pkg->Param("dipole_b_field_mx"); + const Real my = hydro_pkg->Param("dipole_b_field_my"); + const Real mz = hydro_pkg->Param("dipole_b_field_mz"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential", + parthenon::DevExecSpace(), a_kb.s, a_kb.e, a_jb.s, a_jb.e, a_ib.s, a_ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + // Compute and apply potential + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + + const Real r3 = pow(SQR(x) + SQR(y) + SQR(z), 3. / 2); + + const Real m_cross_r_x = my * z - mz * y; + const Real m_cross_r_y = mz * x - mx * z; + const Real m_cross_r_z = mx * y - mx * y; + + A(0, k, j, i) += m_cross_r_x / (4 * M_PI * r3); + A(1, k, j, i) += m_cross_r_y / (4 * M_PI * r3); + A(2, k, j, i) += m_cross_r_z / (4 * M_PI * r3); + }); + } + + /************************************************************ + * Apply the potential to the conserved variables + ************************************************************/ + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::ApplyMagneticPotential", + 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) { + u(IB1, k, j, i) = + (A(2, k, j + 1, i) - A(2, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0 - + (A(1, k + 1, j, i) - A(1, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0; + u(IB2, k, j, i) = + (A(0, k + 1, j, i) - A(0, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0 - + (A(2, k, j, i + 1) - A(2, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0; + u(IB3, k, j, i) = + (A(1, k, j, i + 1) - A(1, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0 - + (A(0, k, j + 1, i) - A(0, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0; + + u(IEN, k, j, i) += 0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + + SQR(u(IB3, k, j, i))); + }); + + /************************************************************ + * Add uniform magnetic field to the conserved variables + ************************************************************/ + const auto &init_uniform_b_field = hydro_pkg->Param("init_uniform_b_field"); + if (init_uniform_b_field) { + const Real bx = hydro_pkg->Param("uniform_b_field_bx"); + const Real by = hydro_pkg->Param("uniform_b_field_by"); + const Real bz = hydro_pkg->Param("uniform_b_field_bz"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::ApplyUniformBField", + 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) { + const Real bx_i = u(IB1, k, j, i); + const Real by_i = u(IB2, k, j, i); + const Real bz_i = u(IB3, k, j, i); + + u(IB1, k, j, i) += bx; + u(IB2, k, j, i) += by; + u(IB3, k, j, i) += bz; + + // Old magnetic energy is b_i^2, new Magnetic energy should be 0.5*(b_i + + // b)^2, add b_i*b + 0.5b^2 to old energy to accomplish that + u(IEN, k, j, i) += + bx_i * bx + by_i * by + bz_i * bz + 0.5 * (SQR(bx) + SQR(by) + SQR(bz)); + }); + // end if(init_uniform_b_field) + } + + } // END if(hydro_pkg->Param("fluid") == Fluid::glmmhd) + } + + /************************************************************ + * Initial parameters + ************************************************************/ + + auto pmb = md->GetBlockData(0)->GetBlockPointer(); + 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 fluid = hydro_pkg->Param("fluid"); + auto const &cons = md->PackVariables(std::vector{"cons"}); + const auto num_blocks = md->NumBlocks(); + + /************************************************************ + * Set initial density perturbations + ************************************************************/ + + const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + const bool full_box = pin->GetOrAddBoolean("problem/cluster/init_perturb", "full_box", true); + 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); + + hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + + Real passive_scalar = 0.0; // Not useful here + + // Spherical collapse test with an initial overdensity + + if (spherical_collapse == 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 r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); // Computing radius + const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); + 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 an initial + u(IDN, k, j, i) = background_density; + + // For any cell at a radius less then overdensity radius, set the density to foreground_density + if (r < overdensity_radius){ + u(IDN, k, j, i) = foreground_density; + } + + //u(IDN, k, j, i) += foreground_density * std::exp(- 2 * SQR(r) / SQR(overdensity_radius) ) ; + + + }, + passive_scalar); + + } + + /* -------------- Setting up a clumpy atmosphere -------------- + + 1) Extract the values of the density from an input hdf5 file using H5Easy + 2) Initiate the associated density field + 3) Optionnaly, add some overpressure ring to check behavior (overpressure_ring bool) + 4) Optionnaly, add a central overdensity + + */ + + if (init_perturb_rho == true) { + + auto init_perturb_rho_file = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_file"); + auto init_perturb_rho_keys = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_keys"); + + hydro_pkg->AddParam<>("cluster/init_perturb_rho_file", init_perturb_rho_file); + hydro_pkg->AddParam<>("cluster/init_perturb_rho_keys", init_perturb_rho_keys); + + std::cout << "Setting density perturbation"; + + // Read HDF5 file containing the density + std::string filename_rho = "/work/bbd0833/test/rho.h5"; + std::string keys_rho = "data"; + H5Easy::File file(filename_rho, HighFive::File::ReadOnly); + auto rho_init = H5Easy::load, 256>, 256>>(file, keys_rho); + + Real passive_scalar = 0.0; // Useless + + std::cout << "entering initialisation of rho field"; + + 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); + + // Adding the value of the density field + if (full_box){ + + u(IDN, k, j, i) += rho_init[gks + k - 2][gjs + j - 2][gis + i - 2]; + + } + + else { + + const Real z = coords.Xc<3>(k); + + // if (z > -thickness_ism / 2 && z < thickness_ism / 2){ + + u(IDN, k, j, i) += rho_init[gks + k - 2][gjs + j - 2][gis + i - 2] * std::exp(- SQR(z) / SQR(thickness_ism)) ; + + + } + + // Ring of overpressure + // compute radius + if (overpressure_ring){ + + 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)); // Computing radius + const Real width = 0.002; + const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); + + if (r > 0 && r < overdensity_radius){ + + //u(IEN, k, j, i) *= 100; + + // Pushing a ring of gas outward (kinetic boost) + const Real u_x = x / r; + const Real u_y = y / r; + const Real u_z = z / r; + + const Real velocity_blast = pin->GetOrAddReal("problem/cluster/init_perturb", "velocity_blast", 0.0); + + //u(IEN, k, j, i) *= 100; + + u(IM1, k, j, i) = u_x * velocity_blast * u(IDN, k, j, i); // p = (1 Mpc / Gyr) * rho + u(IM2, k, j, i) = u_y * velocity_blast * u(IDN, k, j, i); // p = (1 Mpc / Gyr) * rho + u(IM3, k, j, i) = u_z * velocity_blast * u(IDN, k, j, i); // p = (1 Mpc / Gyr) * rho + u(IEN, k, j, i) += 0.5 * u(IDN, k, j, i) * (pow(velocity_blast,2)); + } + } + + if (spherical_collapse){ + + 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)); // Computing radius + const Real overdensity_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overdensity_radius", 0.01); + + u(IDN, k, j, i) = 3; + u(IDN, k, j, i) *= 50 * std::exp(- 2 * SQR(r) / SQR(overdensity_radius)); + + } + + }, + passive_scalar); + + } + + /************************************************************ + * Set initial velocity perturbations (requires no other velocities for now) + ************************************************************/ + + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); + + if (sigma_v != 0.0) { + auto few_modes_ft = hydro_pkg->Param("cluster/few_modes_ft_v"); + // Init phases on all blocks + for (int b = 0; b < md->NumBlocks(); b++) { + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + few_modes_ft.SetPhases(pmb.get(), pin); + + } + // 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.Generate(md, dt, "tmp_perturb"); + + Real v2_sum = 0.0; // used for normalization + + auto perturb_pack = 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); + auto rho = u(IDN, k, j, i); + // The following restriction could be lifted, but requires refactoring of the + // logic for the normalization/reduction below + 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))) * + coords.CellVolume(k, j, i) / SQR(rho); + }, + v2_sum); + +#ifdef MPI_PARALLEL + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &v2_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 v_norm = std::sqrt(v2_sum / (Lx * Ly * Lz) / (SQR(sigma_v))); + + pmb->par_for( + "Norm 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) { + const auto &u = cons(b); + + u(IM1, k, j, i) /= v_norm; + u(IM2, k, j, i) /= v_norm; + u(IM3, k, j, i) /= v_norm; + + u(IEN, k, j, i) += + 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) / + u(IDN, k, j, i); + }); + } + + /************************************************************ + * Set initial magnetic field perturbations (resets magnetic field field) + ************************************************************/ + const auto sigma_b = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_b", 0.0); + if (sigma_b != 0.0) { + auto few_modes_ft = hydro_pkg->Param("cluster/few_modes_ft_b"); + // Init phases on all blocks + for (int b = 0; b < md->NumBlocks(); b++) { + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + few_modes_ft.SetPhases(pmb.get(), pin); + } + // As for t_corr in few_modes_ft, the choice for dt is + // in principle arbitrary because the inital b_hat is 0 and the b_hat_new will contain + // the perturbation (and is normalized in the following to get the desired sigma_b) + const Real dt = 1.0; + few_modes_ft.Generate(md, dt, "tmp_perturb"); + + Real b2_sum = 0.0; // used for normalization + + auto perturb_pack = md->PackVariables(std::vector{"tmp_perturb"}); + + pmb->par_reduce( + "Init 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, Real &lsum) { + const auto &coords = cons.GetCoords(b); + const auto &u = cons(b); + // The following restriction could be lifted, but requires refactoring of the + // logic for the normalization/reduction below + PARTHENON_REQUIRE( + u(IB1, k, j, i) == 0.0 && u(IB2, k, j, i) == 0.0 && u(IB3, k, j, i) == 0.0, + "Found existing non-zero B when setting magnetic field perturbations."); + u(IB1, k, j, i) = + (perturb_pack(b, 2, k, j + 1, i) - perturb_pack(b, 2, k, j - 1, i)) / + coords.Dxc<2>(j) / 2.0 - + (perturb_pack(b, 1, k + 1, j, i) - perturb_pack(b, 1, k - 1, j, i)) / + coords.Dxc<3>(k) / 2.0; + u(IB2, k, j, i) = + (perturb_pack(b, 0, k + 1, j, i) - perturb_pack(b, 0, k - 1, j, i)) / + coords.Dxc<3>(k) / 2.0 - + (perturb_pack(b, 2, k, j, i + 1) - perturb_pack(b, 2, k, j, i - 1)) / + coords.Dxc<1>(i) / 2.0; + u(IB3, k, j, i) = + (perturb_pack(b, 1, k, j, i + 1) - perturb_pack(b, 1, k, j, i - 1)) / + coords.Dxc<1>(i) / 2.0 - + (perturb_pack(b, 0, k, j + 1, i) - perturb_pack(b, 0, k, j - 1, i)) / + coords.Dxc<2>(j) / 2.0; + + // No need to touch the energy yet as we'll normalize later + lsum += (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i))) * + coords.CellVolume(k, j, i); + }, + b2_sum); + +#ifdef MPI_PARALLEL + 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) { + const auto &u = cons(b); + + u(IB1, k, j, i) /= b_norm; + u(IB2, k, j, i) /= b_norm; + u(IB3, k, j, i) /= b_norm; + + u(IEN, k, j, i) += + 0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i))); + }); + } +} + +void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { + // get hydro + auto pkg = pmb->packages.Get("Hydro"); + const Real gam = pin->GetReal("hydro", "gamma"); + const Real gm1 = (gam - 1.0); + + // get prim vars + auto &data = pmb->meshblock_data.Get(); + auto const &prim = data->Get("prim").data; + + // get derived fields + auto &log10_radius = data->Get("log10_cell_radius").data; + auto &entropy = data->Get("entropy").data; + auto &mach_sonic = data->Get("mach_sonic").data; + auto &temperature = data->Get("temperature").data; + + // for computing temperature from primitives + auto units = pkg->Param("units"); + auto mbar_over_kb = pkg->Param("mbar_over_kb"); + auto mbar = mbar_over_kb * units.k_boltzmann(); + + // fill derived vars (*including ghost cells*) + auto &coords = pmb->coords; + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + + pmb->par_for( + "Cluster::UserWorkBeforeOutput", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real v1 = prim(IV1, k, j, i); + const Real v2 = prim(IV2, k, j, i); + const Real v3 = prim(IV3, k, j, i); + const Real P = prim(IPR, k, j, i); + + // compute radius + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r2 = SQR(x) + SQR(y) + SQR(z); + log10_radius(k, j, i) = 0.5 * std::log10(r2); + + // compute entropy + const Real K = P / std::pow(rho / mbar, gam); + entropy(k, j, i) = K; + + const Real v_mag = std::sqrt(SQR(v1) + SQR(v2) + SQR(v3)); + const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS + const Real M_s = v_mag / c_s; + mach_sonic(k, j, i) = M_s; + + // compute temperature + temperature(k, j, i) = mbar_over_kb * P / rho; + }); + + if (pkg->Param("enable_cooling") == Cooling::tabular) { + auto &cooling_time = data->Get("cooling_time").data; + + // get cooling function + const cooling::TabularCooling &tabular_cooling = + pkg->Param("tabular_cooling"); + const auto cooling_table_obj = tabular_cooling.GetCoolingTableObj(); + + pmb->par_for( + "Cluster::UserWorkBeforeOutput::CoolingTime", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real P = prim(IPR, k, j, i); + + // compute cooling time + const Real eint = P / (rho * gm1); + const Real edot = cooling_table_obj.DeDt(eint, rho); + cooling_time(k, j, i) = (edot != 0) ? -eint / edot : NAN; + }); + } + + if (pkg->Param("fluid") == Fluid::glmmhd) { + auto &plasma_beta = data->Get("plasma_beta").data; + auto &mach_alfven = data->Get("mach_alfven").data; + + pmb->par_for( + "Cluster::UserWorkBeforeOutput::MHD", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real P = prim(IPR, k, j, i); + const Real Bx = prim(IB1, k, j, i); + const Real By = prim(IB2, k, j, i); + const Real Bz = prim(IB3, k, j, i); + const Real B2 = (SQR(Bx) + SQR(By) + SQR(Bz)); + + // compute Alfven mach number + const Real v_A = std::sqrt(B2 / rho); + const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS + mach_alfven(k, j, i) = mach_sonic(k, j, i) * c_s / v_A; + + // compute plasma beta + plasma_beta(k, j, i) = (B2 != 0) ? P / (0.5 * B2) : NAN; + }); + } +} + +} // namespace cluster diff --git a/src/utils/few_modes_ft_lognormal.cpp b/src/utils/few_modes_ft_lognormal.cpp new file mode 100644 index 00000000..d47b30c9 --- /dev/null +++ b/src/utils/few_modes_ft_lognormal.cpp @@ -0,0 +1,431 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//======================================================================================== +//! \file few_modes_ft.cpp +// \brief Helper functions for an inverse (explicit complex to real) FT + +// C++ headers +#include +#include + +// Parthenon headers +#include "basic_types.hpp" +#include "config.hpp" +#include "globals.hpp" +#include "kokkos_abstraction.hpp" +#include "mesh/domain.hpp" +#include "mesh/meshblock_pack.hpp" + +// AthenaPK headers +#include "../main.hpp" +#include "few_modes_ft_lognormal.hpp" +#include "utils/error_checking.hpp" + +namespace utils::few_modes_ft_log { +using Complex = Kokkos::complex; +using parthenon::IndexRange; +using parthenon::Metadata; + +FewModesFTLog::FewModesFTLog(parthenon::ParameterInput *pin, parthenon::StateDescriptor *pkg, + 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), + t_corr_(t_corr), fill_ghosts_(fill_ghosts) { + + 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 + << "If many modes are required in the transform field consider using " + << "the driving mechanism based on full FFTs." << std::endl; + } + // Ensure that all all wavevectors can be represented on the root grid + const auto gnx1 = pin->GetInteger("parthenon/mesh", "nx1"); + const auto gnx2 = pin->GetInteger("parthenon/mesh", "nx2"); + const auto gnx3 = pin->GetInteger("parthenon/mesh", "nx3"); + // 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"); + const auto nx3 = pin->GetInteger("parthenon/meshblock", "nx3"); + const auto ng_tot = fill_ghosts_ ? 2 * parthenon::Globals::nghost : 0; + auto m = Metadata({Metadata::None, Metadata::Derived, Metadata::OneCopy}, + std::vector({2, num_modes, nx1 + ng_tot}), prefix + "_phases_i"); + pkg->AddField(prefix + "_phases_i", m); + m = Metadata({Metadata::None, Metadata::Derived, Metadata::OneCopy}, + std::vector({2, num_modes, nx2 + ng_tot}), prefix + "_phases_j"); + pkg->AddField(prefix + "_phases_j", m); + m = Metadata({Metadata::None, Metadata::Derived, Metadata::OneCopy}, + std::vector({2, num_modes, nx3 + ng_tot}), prefix + "_phases_k"); + pkg->AddField(prefix + "_phases_k", m); + + // Variable (e.g., acceleration field for turbulence driver) in Fourier space using + // complex to real transform. + var_hat_ = ParArray2D(prefix + "_var_hat", 3, num_modes); + var_hat_new_ = ParArray2D(prefix + "_var_hat_new", 3, num_modes); + + PARTHENON_REQUIRE((sol_weight == -1.0) || (sol_weight >= 0.0 && sol_weight <= 1.0), + "sol_weight for projection in few modes fft module needs to be " + "between 0.0 and 1.0 or set to -1.0 (to disable projection).") + sol_weight_ = sol_weight; + + random_num_ = Kokkos::View( + "random_num", 3, num_modes, 2); + random_num_host_ = Kokkos::create_mirror_view(random_num_); + + rng_.seed(rseed); + dist_ = std::uniform_real_distribution<>(-1.0, 1.0); +} + +void FewModesFTLog::SetPhases(MeshBlock *pmb, ParameterInput *pin) { + auto pm = pmb->pmy_mesh; + auto hydro_pkg = pmb->packages.Get("Hydro"); + + // The following restriction could technically be lifted if the turbulence driver is + // directly embedded in the hydro driver rather than a user defined source as well as + // fixing the pack_size=-1 when using the Mesh- (not MeshBlock-)based problem generator. + // The restriction stems from requiring a collective MPI comm to normalize the + // acceleration and magnetic field, respectively. Note, that the restriction does not + // apply here, but for the ProblemGenerator() and Driving() function below. The check is + // just added here for convenience as this function is called during problem + // initializtion. From my (pgrete) point of view, it's currently cleaner to keep things + // 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 " + "to work because of global reductions.") + + auto Lx1 = pm->mesh_size.x1max - pm->mesh_size.x1min; + auto Lx2 = pm->mesh_size.x2max - pm->mesh_size.x2min; + auto Lx3 = pm->mesh_size.x3max - pm->mesh_size.x3min; + + // Adjust (logical) grid size at levels other than the root level. + // This is required for simulation with mesh refinement so that the phases calculated + // 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); // Size of cell + 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 + PARTHENON_REQUIRE_THROWS(((gnx1 == gnx2) && (gnx2 == gnx3)) && + ((Lx1 == Lx2) && (Lx2 == Lx3)), + "FMFT has only been tested with cubic meshes and constant " + "dx/dy/dz. Remove this warning at your own risk.") + + const auto nx1 = pmb->block_size.nx1; + 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; + + // 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); + + auto &base = pmb->meshblock_data.Get(); + auto &phases_i = base->Get(prefix_ + "_phases_i").data; + 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) { + Real gi = static_cast((i + gis - ng) % static_cast(gnx1)); + 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 { + phase = Kokkos::exp(I * w_kx * gi); + } + phases_i(i, m, 0) = phase.real(); + phases_i(i, m, 1) = phase.imag(); + + std::cout << "SetPhases: phases computed \n" ; + } + }); + + pmb->par_for( + "FMFT: calc phases_j", 0, nx2 - 1 + 2 * ng, KOKKOS_LAMBDA(int j) { + Real gj = static_cast((j + gjs - ng) % static_cast(gnx2)); + Real w_ky; + Complex phase; + + for (int m = 0; m < num_modes; m++) { + w_ky = k_vec(1, m) * 2. * M_PI / static_cast(gnx2); + phase = Kokkos::exp(I * w_ky * gj); + phases_j(j, m, 0) = phase.real(); + phases_j(j, m, 1) = phase.imag(); + } + }); + + pmb->par_for( + "FMFT: calc phases_k", 0, nx3 - 1 + 2 * ng, KOKKOS_LAMBDA(int k) { + Real gk = static_cast((k + gks - ng) % static_cast(gnx3)); + Real w_kz; + Complex phase; + + for (int m = 0; m < num_modes; m++) { + w_kz = k_vec(2, m) * 2. * M_PI / static_cast(gnx3); + phase = Kokkos::exp(I * w_kz * gk); + phases_k(k, m, 0) = phase.real(); + phases_k(k, m, 1) = phase.imag(); + } + }); +} + +void FewModesFTLog::Generate(MeshData *md, const Real dt, + const std::string &var_name) { + auto pmb = md->GetBlockData(0)->GetBlockPointer(); + + const auto num_modes = num_modes_; + + Complex I(0.0, 1.0); + auto &random_num = random_num_; + + // get a set of random numbers from the CPU so that they are deterministic + // when run on GPUs + Real v1, v2, v_sqr; + for (int n = 0; n < 3; n++) + for (int m = 0; m < num_modes; m++) { + do { + v1 = dist_(rng_); + v2 = dist_(rng_); + v_sqr = v1 * v1 + v2 * v2; + } while (v_sqr >= 1.0 || v_sqr == 0.0); + + random_num_host_(n, m, 0) = v1; + random_num_host_(n, m, 1) = v2; + } + Kokkos::deep_copy(random_num, random_num_host_); + + // make local ref to capure in lambda + auto &k_vec = k_vec_; + auto &var_hat = var_hat_; + auto &var_hat_new = var_hat_new_; + + const auto kmin = k_min_; + + // generate new power spectrum (injection) + pmb->par_for( + "FMFT: new power spec", 0, 2, 0, num_modes - 1, + KOKKOS_LAMBDA(const int n, const int m) { + Real kmag, tmp, norm, v_sqr; + + Real kx = k_vec(0, m); + Real ky = k_vec(1, m); + Real kz = k_vec(2, m); + + kmag = std::sqrt(kx * kx + ky * ky + kz * kz); + + var_hat_new(n, m) = Complex(0., 0.); + + tmp = std::pow(kmag / kmin, -5./3.); // AMPLITUDE, TO BE REDEFINED TO LOGNORMAL DISTRIBUTION + if (tmp < 0.) tmp = 0.; + v_sqr = SQR(random_num(n, m, 0)) + SQR(random_num(n, m, 1)); + norm = std::sqrt(-2.0 * std::log(v_sqr) / v_sqr); + + var_hat_new(n, m) = + Complex(tmp * norm * random_num(n, m, 0), tmp * norm * random_num(n, m, 1)); + }); + + // enforce symmetry of complex to real transform + pmb->par_for( + "forcing: enforce symmetry", 0, 2, 0, num_modes - 1, + KOKKOS_LAMBDA(const int n, const int m) { + if (k_vec(0, m) == 0.) { + for (int m2 = 0; m2 < m; m2++) { + if (k_vec(1, m) == -k_vec(1, m2) && k_vec(2, m) == -k_vec(2, m2)) + var_hat_new(n, m) = + Complex(var_hat_new(n, m2).real(), -var_hat_new(n, m2).imag()); + } + } + }); + + const auto sol_weight = sol_weight_; + if (sol_weight_ >= 0.0) { + // project + pmb->par_for( + "forcing: projection", 0, num_modes - 1, KOKKOS_LAMBDA(const int m) { + Real kmag; + + Real kx = k_vec(0, m); + Real ky = k_vec(1, m); + Real kz = k_vec(2, m); + + kmag = std::sqrt(kx * kx + ky * ky + kz * kz); + + // setting kmag to 1 as a "continue" doesn't work within the parallel_for + // construct and it doesn't affect anything (there should never be power in the + // k=0 mode) + if (kmag == 0.) kmag = 1.; + + // make it a unit vector + kx /= kmag; + ky /= kmag; + kz /= kmag; + + Complex dot(var_hat_new(0, m).real() * kx + var_hat_new(1, m).real() * ky + + var_hat_new(2, m).real() * kz, + var_hat_new(0, m).imag() * kx + var_hat_new(1, m).imag() * ky + + var_hat_new(2, m).imag() * kz); + + var_hat_new(0, m) = Complex(var_hat_new(0, m).real() * sol_weight + + (1. - 2. * sol_weight) * dot.real() * kx, + var_hat_new(0, m).imag() * sol_weight + + (1. - 2. * sol_weight) * dot.imag() * kx); + var_hat_new(1, m) = Complex(var_hat_new(1, m).real() * sol_weight + + (1. - 2. * sol_weight) * dot.real() * ky, + var_hat_new(1, m).imag() * sol_weight + + (1. - 2. * sol_weight) * dot.imag() * ky); + var_hat_new(2, m) = Complex(var_hat_new(2, m).real() * sol_weight + + (1. - 2. * sol_weight) * dot.real() * kz, + var_hat_new(2, m).imag() * sol_weight + + (1. - 2. * sol_weight) * dot.imag() * kz); + }); + } + + // evolve + const auto c_drift = std::exp(-dt / t_corr_); + const auto c_diff = std::sqrt(1.0 - c_drift * c_drift); + + pmb->par_for( + "FMFT: evolve spec", 0, 2, 0, num_modes - 1, + KOKKOS_LAMBDA(const int n, const int m) { + var_hat(n, m) = + Complex(var_hat(n, m).real() * c_drift + var_hat_new(n, m).real() * c_diff, + var_hat(n, m).imag() * c_drift + var_hat_new(n, m).imag() * c_diff); + }); + + auto domain = fill_ghosts_ ? IndexDomain::entire : IndexDomain::interior; + IndexRange ib = md->GetBlockData(0)->GetBoundsI(domain); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(domain); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(domain); + auto var_pack = md->PackVariables(std::vector{var_name}); + auto phases_i = md->PackVariables(std::vector{prefix_ + "_phases_i"}); + auto phases_j = md->PackVariables(std::vector{prefix_ + "_phases_j"}); + auto phases_k = md->PackVariables(std::vector{prefix_ + "_phases_k"}); + + // implictly assuming cubic box of size L=1 + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "FMFT: Inverse FT", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, 0, 2, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i) { + Complex phase, phase_i, phase_j, phase_k; + var_pack(b, n, k, j, i) = 0.0; + + for (int m = 0; m < num_modes; m++) { + phase_i = + Complex(phases_i(b, 0, i - ib.s, m, 0), phases_i(b, 0, i - ib.s, m, 1)); + phase_j = + Complex(phases_j(b, 0, j - jb.s, m, 0), phases_j(b, 0, j - jb.s, m, 1)); + phase_k = + Complex(phases_k(b, 0, k - kb.s, m, 0), phases_k(b, 0, k - kb.s, m, 1)); + phase = phase_i * phase_j * phase_k; + var_pack(b, n, k, j, i) += 2. * (var_hat(n, m).real() * phase.real() - + var_hat(n, m).imag() * phase.imag()); + } + }); +} + +// Creates a random set of wave vectors from k_min to k_max (limit imposed by cell size) +ParArray2D MakeRandomModesLog(const int num_modes, const Real k_min, const Real k_max, + uint32_t rseed = 31224) { + + auto k_vec = parthenon::ParArray2D("k_vec", 3, num_modes); + auto k_vec_h = Kokkos::create_mirror_view_and_copy(parthenon::HostMemSpace(), k_vec); + + const int k_low = k_min; + const int k_high = k_max; + + std::mt19937 rng; + rng.seed(rseed); + std::uniform_int_distribution<> dist(-k_high, k_high); // Function to be used to define size of the + + int n_mode = 0; + int n_attempt = 0; + constexpr int max_attempts = 1000000; + Real kx1, kx2, kx3, 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); + k_mag = std::sqrt(SQR(kx1) + SQR(kx2) + SQR(kx3)); + + // Expected amplitude of the spectral function. If this is changed, it also needs to + // be changed in the FMFT class (or abstracted). + ampl = std::pow((k_mag / k_low), -5./3.); + + // Check is mode was already picked by chance + mode_exists = false; + for (int n_mode_exsist = 0; n_mode_exsist < n_mode; n_mode_exsist++) { + if (k_vec_h(0, n_mode_exsist) == kx1 && k_vec_h(1, n_mode_exsist) == kx2 && + k_vec_h(2, n_mode_exsist) == kx3) { + mode_exists = true; + } + } + + // kx1 < 0.0 because we use a explicit symmetric Complex to Real transform + if (ampl < 0 || k_mag < k_low || k_mag > k_high || mode_exists || kx1 < 0.0) { + continue; + } + k_vec_h(0, n_mode) = kx1; + k_vec_h(1, n_mode) = kx2; + k_vec_h(2, n_mode) = kx3; + n_mode++; + } + PARTHENON_REQUIRE_THROWS( + n_attempt < max_attempts, + "Cluster init did not succeed in calculating perturbation modes.") + Kokkos::deep_copy(k_vec, k_vec_h); + + return k_vec; +} + + +} // namespace utils::few_modes_ft \ No newline at end of file diff --git a/src/utils/few_modes_ft_lognormal.hpp b/src/utils/few_modes_ft_lognormal.hpp new file mode 100644 index 00000000..a8823947 --- /dev/null +++ b/src/utils/few_modes_ft_lognormal.hpp @@ -0,0 +1,72 @@ + +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//======================================================================================== +//! \file few_modes_ft.hpp +// \brief Helper functions for an inverse (explicit complex to real) FT + +// Parthenon headers +#include "basic_types.hpp" +#include "config.hpp" +#include +#include +#include + +// AthenaPK headers +#include "../main.hpp" +#include "mesh/domain.hpp" + +namespace utils::few_modes_ft_log { +using namespace parthenon::package::prelude; +using parthenon::Real; +using Complex = Kokkos::complex; +using parthenon::IndexRange; +using parthenon::ParArray2D; + +class FewModesFTLog { + private: + int num_modes_; + std::string prefix_; + ParArray2D var_hat_, var_hat_new_; + ParArray2D k_vec_; + Real k_min_; // minimum wave vector + Real k_max_; // maximum wave vector + Kokkos::View random_num_; + Kokkos::View random_num_host_; + std::mt19937 rng_; + std::uniform_real_distribution<> dist_; + Real sol_weight_; // power in solenoidal modes for projection. Set to negative to + // disable projection + Real t_corr_; // correlation time for evolution of Ornstein-Uhlenbeck process + bool fill_ghosts_; // if the inverse transform should also fill ghost zones + + public: + FewModesFTLog(parthenon::ParameterInput *pin, parthenon::StateDescriptor *pkg, + 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 = false); + + ParArray2D GetVarHat() { return var_hat_; } + int GetNumModes() { return num_modes_; } + void SetPhases(MeshBlock *pmb, ParameterInput *pin); + void Generate(MeshData *md, const Real dt, const std::string &var_name); + void RestoreRNG(std::istringstream &iss) { iss >> rng_; } + void RestoreDist(std::istringstream &iss) { iss >> dist_; } + std::string GetRNGState() { + std::ostringstream oss; + oss << rng_; + return oss.str(); + } + std::string GetDistState() { + std::ostringstream oss; + oss << dist_; + return oss.str(); + } +}; + +// Creates a random set of wave vectors with k_mag within k_min and k_max (Nyquist lim.) +ParArray2D MakeRandomModesLog(const int num_modes, const Real k_min, const Real k_max, uint32_t rseed); + +} // namespace utils::few_modes_ft \ No newline at end of file