Skip to content

Commit

Permalink
Fixing cluster.cpp
Browse files Browse the repository at this point in the history
  • Loading branch information
HPC user committed Oct 11, 2023
1 parent 94a302a commit f188f33
Showing 1 changed file with 0 additions and 125 deletions.
125 changes: 0 additions & 125 deletions src/pgen/cluster.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,131 +64,6 @@ using namespace parthenon::package::prelude;
using utils::few_modes_ft::FewModesFT;
using utils::few_modes_ft_log::FewModesFTLog;

<<<<<<< HEAD
template <class EOS>
void ApplyClusterClips(MeshData<Real> *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<Real>("cluster_dfloor");
const auto &eceil = hydro_pkg->Param<Real>("cluster_eceil");
const auto &vceil = hydro_pkg->Param<Real>("cluster_vceil");
const auto &vAceil = hydro_pkg->Param<Real>("cluster_vAceil");
const auto &clip_r = hydro_pkg->Param<Real>("cluster_clip_r");

if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits<Real>::infinity() ||
vceil < std::numeric_limits<Real>::infinity() ||
vAceil < std::numeric_limits<Real>::infinity())) {
// Grab some necessary variables
const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
const auto &cons_pack = md->PackVariables(std::vector<std::string>{"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<int>("nhydro");
const auto nscalars = hydro_pkg->Param<int>("nscalars");

const Real clip_r2 = SQR(clip_r);
const Real vceil2 = SQR(vceil);
const Real vAceil2 = SQR(vAceil);
const Real gm1 = (hydro_pkg->Param<Real>("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); // 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<Real>::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<Real>::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<Real>::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<Real> *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>("fluid");
if (fluid == Fluid::euler) {
ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param<AdiabaticHydroEOS>("eos"));
} else if (fluid == Fluid::glmmhd) {
ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos"));
} else {
PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS");
}
}

//void ClusterSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
// const Real beta_dt) {

void ClusterUnsplitSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
const Real beta_dt) {

Expand Down

0 comments on commit f188f33

Please sign in to comment.