Skip to content

Commit

Permalink
Formating
Browse files Browse the repository at this point in the history
  • Loading branch information
HPC user committed Dec 13, 2023
1 parent 1529075 commit b956414
Show file tree
Hide file tree
Showing 24 changed files with 115 additions and 12,068 deletions.
1 change: 0 additions & 1 deletion src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@ add_executable(
refinement/gradient.cpp
refinement/other.cpp
utils/few_modes_ft.cpp
utils/few_modes_ft_lognormal.cpp
)

add_subdirectory(pgen)
Expand Down
7 changes: 2 additions & 5 deletions src/eos/adiabatic_glmmhd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,6 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData<Real> *md) const {
const auto nhydro = pkg->Param<int>("nhydro");
const auto nscalars = pkg->Param<int>("nscalars");


auto this_on_device = (*this);

parthenon::par_for(
Expand All @@ -51,15 +50,13 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData<Real> *md) const {
const auto &cons = cons_pack(b);
auto &prim = prim_pack(b);
// auto &nu = entropy_pack(b);



// Getting the global indexing

auto pmb = md->GetBlockData(b)->GetBlockPointer();
auto pm = pmb->pmy_mesh;
auto hydro_pkg = pmb->packages.Get("Hydro");


return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
});
}
63 changes: 32 additions & 31 deletions src/eos/adiabatic_glmmhd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,10 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
gamma_{gamma} {}

void ConservedToPrimitive(MeshData<Real> *md) const override;

KOKKOS_INLINE_FUNCTION
Real GetGamma() const { return gamma_; }

//----------------------------------------------------------------------------------------
// \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO])
// \brief returns adiabatic sound speed given vector of primitive variables
Expand All @@ -53,7 +53,7 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
return std::sqrt(0.5 * (qsq + std::sqrt(tmp * tmp + 4.0 * asq * ct2)) / d);
}
//

//----------------------------------------------------------------------------------------
// \!fn Real EquationOfState::ConsToPrim(View4D cons, View4D prim, const int& k, const
// int& j, const int& i) \brief Fills an array of primitives given an array of
Expand All @@ -70,51 +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,
// PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0,
// "Got negative density. Consider enabling first-order flux "
// "correction or setting a reasonable density floor.");

// apply density floor, without changing momentum or energy
u_d = (u_d > density_floor_) ? u_d : density_floor_;
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 @@ -131,12 +131,13 @@ 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,
// PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
// "Got negative pressure. Consider enabling first-order flux "
// "correction or setting a reasonable pressure or temperature floor.");
// "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
13 changes: 6 additions & 7 deletions src/eos/adiabatic_hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,27 +37,26 @@ void AdiabaticHydroEOS::ConservedToPrimitive(MeshData<Real> *md) const {
auto ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire);
auto jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire);
auto kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire);

auto pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
const auto nhydro = pkg->Param<int>("nhydro");
const auto nscalars = pkg->Param<int>("nscalars");

auto this_on_device = (*this);

parthenon::par_for(
DEFAULT_LOOP_PATTERN, "ConservedToPrimitive", parthenon::DevExecSpace(), 0,
cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {

// Getting the global indexing

auto pmb = md->GetBlockData(b)->GetBlockPointer();
auto pm = pmb->pmy_mesh;
auto hydro_pkg = pmb->packages.Get("Hydro");

const auto &cons = cons_pack(b);
auto &prim = prim_pack(b);

return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
});
}
26 changes: 13 additions & 13 deletions src/eos/adiabatic_hydro.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,27 +52,27 @@ class AdiabaticHydroEOS : public EquationOfState {
KOKKOS_INLINE_FUNCTION void ConsToPrim(View4D cons, View4D prim, const int &nhydro,
const int &nscalars, const int &k, const int &j,
const int &i) const {

Real gm1 = GetGamma() - 1.0;
auto density_floor_ = GetDensityFloor();
auto pressure_floor_ = GetPressureFloor();
auto e_floor_ = GetInternalEFloor();

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,10 +86,10 @@ 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);

// 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 @@ -106,10 +106,10 @@ class AdiabaticHydroEOS : 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.

// The first argument check whether one of the conditions is true
// By default, floors are deactivated, ie. pressure floor and e_floor
// are -1. The code will eventually crash when w_p turns < 0.
Expand Down
18 changes: 11 additions & 7 deletions src/hydro/srcterms/gravitational_field.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,18 +36,22 @@ void GravitationalFieldSrcTerm(parthenon::MeshData<parthenon::Real> *md,
IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);

/*
std::cout << "Values of gravitationnal field" << std::endl;
std::cout
std::cout
<< gravitationalField.g_from_r(1e-4) << ", " << gravitationalField.g_from_r(3.16e-4)
<< ", " << gravitationalField.g_from_r(1e-3) << ", " << gravitationalField.g_from_r(3.16e-3)
<< ", " << gravitationalField.g_from_r(1e-2) << ", " << gravitationalField.g_from_r(3.16e-2)
<< ", " << gravitationalField.g_from_r(1e-1) << ", " << gravitationalField.g_from_r(3.16e-1)
<< ", " << gravitationalField.g_from_r(1.0) << ", " << gravitationalField.g_from_r(3.16)
<< ", " << gravitationalField.g_from_r(1e-3) << ", " <<
gravitationalField.g_from_r(3.16e-3)
<< ", " << gravitationalField.g_from_r(1e-2) << ", " <<
gravitationalField.g_from_r(3.16e-2)
<< ", " << gravitationalField.g_from_r(1e-1) << ", " <<
gravitationalField.g_from_r(3.16e-1)
<< ", " << gravitationalField.g_from_r(1.0) << ", " <<
gravitationalField.g_from_r(3.16)
<< ", " << gravitationalField.g_from_r(10.0) << std::endl;
*/

parthenon::par_for(
DEFAULT_LOOP_PATTERN, "GravitationalFieldSrcTerm", parthenon::DevExecSpace(), 0,
cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
Expand Down
Loading

0 comments on commit b956414

Please sign in to comment.