Skip to content

Commit

Permalink
Cleaning
Browse files Browse the repository at this point in the history
  • Loading branch information
HPC user committed Nov 30, 2023
1 parent e4ed473 commit 06c5835
Show file tree
Hide file tree
Showing 21 changed files with 6,483 additions and 373 deletions.
5 changes: 0 additions & 5 deletions inputs/cluster/hse.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ nlim = -1 # cycle limit
tlim = 1e-3 # time limit
integrator = vl2 # time integration algorithm


<parthenon/mesh>
refinement = static
nghost = 2
Expand All @@ -47,7 +46,6 @@ x3max = 0.1 # maximum value of X3
ix3_bc = outflow # inner-X3 boundary flag
ox3_bc = outflow # outer-X3 boundary flag


<parthenon/static_refinement0>
x1min = -0.0125
x1max = 0.0125
Expand All @@ -57,7 +55,6 @@ x3min = -0.0125
x3max = 0.0125
level = 2


<parthenon/meshblock>
nx1 = 32 # Number of zones in X1-direction
nx2 = 32 # Number of zones in X2-direction
Expand Down Expand Up @@ -106,7 +103,6 @@ g_smoothing_radius = 1e-6
#Include gravity as a source term
gravity_srcterm = true


<problem/cluster/entropy_profile>
#Entropy profile parameters
k_0 = 8.851337676479303e-121
Expand All @@ -115,7 +111,6 @@ r_k = 0.1
alpha_k = 1.1

<problem/cluster/hydrostatic_equilibrium>

#Fix density at radius to close system of equations
r_fix = 2.0
rho_fix = 0.01477557589278723
Expand Down
53 changes: 27 additions & 26 deletions src/eos/adiabatic_glmmhd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,50 +70,51 @@ class AdiabaticGLMMHDEOS : public EquationOfState {

auto velocity_ceiling_ = GetVelocityCeiling();
auto e_ceiling_ = GetInternalECeiling();

Real &u_d = cons(IDN, k, j, i);
Real &u_m1 = cons(IM1, k, j, i);
Real &u_m2 = cons(IM2, k, j, i);
Real &u_m3 = cons(IM3, k, j, i);
Real &u_e = cons(IEN, k, j, i);
Real &u_b1 = cons(IB1, k, j, i);
Real &u_b2 = cons(IB2, k, j, i);
Real &u_b3 = cons(IB3, k, j, i);
Real &u_d = cons(IDN, k, j, i);
Real &u_m1 = cons(IM1, k, j, i);
Real &u_m2 = cons(IM2, k, j, i);
Real &u_m3 = cons(IM3, k, j, i);
Real &u_e = cons(IEN, k, j, i);
Real &u_b1 = cons(IB1, k, j, i);
Real &u_b2 = cons(IB2, k, j, i);
Real &u_b3 = cons(IB3, k, j, i);
Real &u_psi = cons(IPS, k, j, i);

Real &w_d = prim(IDN, k, j, i);
Real &w_vx = prim(IV1, k, j, i);
Real &w_vy = prim(IV2, k, j, i);
Real &w_vz = prim(IV3, k, j, i);
Real &w_p = prim(IPR, k, j, i);
Real &w_Bx = prim(IB1, k, j, i);
Real &w_By = prim(IB2, k, j, i);
Real &w_Bz = prim(IB3, k, j, i);
Real &w_d = prim(IDN, k, j, i);
Real &w_vx = prim(IV1, k, j, i);
Real &w_vy = prim(IV2, k, j, i);
Real &w_vz = prim(IV3, k, j, i);
Real &w_p = prim(IPR, k, j, i);
Real &w_Bx = prim(IB1, k, j, i);
Real &w_By = prim(IB2, k, j, i);
Real &w_Bz = prim(IB3, k, j, i);
Real &w_psi = prim(IPS, k, j, i);

// 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 reasonble density floor.");
"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_;
w_d = u_d;

Real di = 1.0 / u_d;
w_vx = u_m1 * di;
w_vy = u_m2 * di;
w_vz = u_m3 * di;

w_Bx = u_b1;
w_By = u_b2;
w_Bz = u_b3;
w_psi = u_psi;

Real e_k = 0.5 * di * (SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
Real e_B = 0.5 * (SQR(u_b1) + SQR(u_b2) + SQR(u_b3));
w_p = gm1 * (u_e - e_k - e_B);

// apply velocity ceiling. By default ceiling is std::numeric_limits<Real>::infinity()
const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz);
if (w_v2 > SQR(velocity_ceiling_)) {
Expand All @@ -130,12 +131,12 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
u_e -= e_k - e_k_new;
e_k = e_k_new;
}

// 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 reasonble pressure or temperature floor.");
"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_)) {
Expand Down
14 changes: 7 additions & 7 deletions src/eos/adiabatic_hydro.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,18 +61,18 @@ class AdiabaticHydroEOS : public EquationOfState {
auto velocity_ceiling_ = GetVelocityCeiling();
auto e_ceiling_ = GetInternalECeiling();

Real &u_d = cons(IDN, k, j, i);
Real &u_d = cons(IDN, k, j, i);
Real &u_m1 = cons(IM1, k, j, i);
Real &u_m2 = cons(IM2, k, j, i);
Real &u_m3 = cons(IM3, k, j, i);
Real &u_e = cons(IEN, k, j, i);

Real &w_d = prim(IDN, k, j, i);
Real &u_e = cons(IEN, k, j, i);
Real &w_d = prim(IDN, k, j, i);
Real &w_vx = prim(IV1, k, j, i);
Real &w_vy = prim(IV2, k, j, i);
Real &w_vz = prim(IV3, k, j, i);
Real &w_p = prim(IPR, k, j, i);
Real &w_p = prim(IPR, k, j, i);

// 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,
Expand All @@ -86,7 +86,7 @@ class AdiabaticHydroEOS : public EquationOfState {
w_vx = u_m1 * di;
w_vy = u_m2 * di;
w_vz = u_m3 * di;

Real e_k = 0.5 * di * (SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
w_p = gm1 * (u_e - e_k);

Expand Down
20 changes: 10 additions & 10 deletions src/hydro/srcterms/tabular_cooling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,8 @@ TabularCooling::TabularCooling(ParameterInput *pin,
} else if (integrator_str == "rk45") {
std::cout << "Cooling integrator is rk45\n";
integrator_ = CoolIntegrator::rk45;
} else if (integrator_str == "townsend\n") {
std::cout << "Cooling integrator is Townsend";
} else if (integrator_str == "townsend") {
std::cout << "Cooling integrator is Townsend\n";
integrator_ = CoolIntegrator::townsend;
} else {
integrator_ = CoolIntegrator::undefined;
Expand Down Expand Up @@ -226,7 +226,7 @@ TabularCooling::TabularCooling(ParameterInput *pin,
if (integrator_ == CoolIntegrator::townsend) {
lambdas_ = ParArray1D<Real>("lambdas_", n_temp_);
temps_ = ParArray1D<Real>("temps_", n_temp_);

// Read log_lambdas in host_lambdas, changing to code units along the way
auto host_lambdas = Kokkos::create_mirror_view(lambdas_);
auto host_temps = Kokkos::create_mirror_view(temps_);
Expand All @@ -238,7 +238,7 @@ TabularCooling::TabularCooling(ParameterInput *pin,
// Copy host_lambdas into device memory
Kokkos::deep_copy(lambdas_, host_lambdas);
Kokkos::deep_copy(temps_, host_temps);

// Coeffs are for intervals, i.e., only n_temp_ - 1 entries
const auto n_bins = n_temp_ - 1;
townsend_Y_k_ = ParArray1D<Real>("townsend_Y_k_", n_bins);
Expand Down Expand Up @@ -504,7 +504,7 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,

// Grab member variables for compiler
const auto dt = dt_; // HACK capturing parameters still broken with Cuda 11.6 ...

const auto units = hydro_pkg->Param<Units>("units");
const auto gm1 = (hydro_pkg->Param<Real>("AdiabaticIndex") - 1.0);
const auto mbar_gm1_over_kb = hydro_pkg->Param<Real>("mbar_over_kb") * gm1;
Expand All @@ -515,10 +515,10 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
const auto temps = temps_;
const auto alpha_k = townsend_alpha_k_;
const auto Y_k = townsend_Y_k_;

const auto internal_e_floor = T_floor_ / mbar_gm1_over_kb;
const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table

// 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"});
Expand All @@ -527,13 +527,13 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire);
IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire);
IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire);

const auto nbins = alpha_k.extent_int(0);

// Get reference values
const auto temp_final = std::pow(10.0, log_temp_final_);
const auto lambda_final = lambda_final_;

par_for(
DEFAULT_LOOP_PATTERN, "TabularCooling::TownsendSrcTerm", DevExecSpace(), 0,
cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
Expand Down
1 change: 1 addition & 0 deletions src/hydro/srcterms/tabular_cooling.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,7 @@ class CoolingTableObj {
const Real lambda = pow(10., log_lambda);
const Real gamma = gamma_units_;
const Real de_dt = -lambda * x_H_over_m_h2_ * rho + gamma*pow(x_H_over_m_h2_,0.5);
//const Real de_dt = -lambda * x_H_over_m_h2_ * rho;

//std::cout << "cool = " << -lambda * x_H_over_m_h2_ * rho << "\t heat = " << gamma*pow(x_H_over_m_h2_,0.5) << "\t e = " << e << "\t rho = " << rho << "\n" ;

Expand Down
Loading

0 comments on commit 06c5835

Please sign in to comment.