Skip to content

Commit

Permalink
Updated Parthenon to PR930
Browse files Browse the repository at this point in the history
  • Loading branch information
Forrest Wolfgang Glines committed Aug 29, 2023
1 parent 160d10a commit 6853b2d
Show file tree
Hide file tree
Showing 9 changed files with 78 additions and 66 deletions.
2 changes: 1 addition & 1 deletion external/parthenon
Submodule parthenon updated from 5b6bb6 to 60d547
6 changes: 3 additions & 3 deletions src/hydro/hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -794,8 +794,8 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
int il, iu, jl, ju, kl, ku;
jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e;
// TODO(pgrete): are these looop limits are likely too large for 2nd order
if (pmb->block_size.nx2 > 1) {
if (pmb->block_size.nx3 == 1) // 2D
if (pmb->block_size.nx(X2DIR) > 1) {
if (pmb->block_size.nx(X3DIR)== 1) // 2D
jl = jb.s - 1, ju = jb.e + 1, kl = kb.s, ku = kb.e;
else // 3D
jl = jb.s - 1, ju = jb.e + 1, kl = kb.s - 1, ku = kb.e + 1;
Expand Down Expand Up @@ -867,7 +867,7 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
parthenon::ScratchPad2D<Real>::shmem_size(num_scratch_vars, nx1) * 3;
// set the loop limits
il = ib.s - 1, iu = ib.e + 1, kl = kb.s, ku = kb.e;
if (pmb->block_size.nx3 == 1) // 2D
if (pmb->block_size.nx(X3DIR) == 1) // 2D
kl = kb.s, ku = kb.e;
else // 3D
kl = kb.s - 1, ku = kb.e + 1;
Expand Down
64 changes: 36 additions & 28 deletions src/pgen/cluster.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,9 @@ using namespace parthenon::driver::prelude;
using namespace parthenon::package::prelude;
using utils::few_modes_ft::FewModesFT;


void ClusterUnsplitSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
const Real beta_dt) {
const Real beta_dt) {
auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");

const bool &gravity_srcterm = hydro_pkg->Param<bool>("gravity_srcterm");
Expand All @@ -78,11 +79,17 @@ void ClusterUnsplitSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
const auto &magnetic_tower = hydro_pkg->Param<MagneticTower>("magnetic_tower");
magnetic_tower.FixedFieldSrcTerm(md, beta_dt, tm);

const auto &snia_feedback = hydro_pkg->Param<SNIAFeedback>("snia_feedback");
snia_feedback.FeedbackSrcTerm(md, beta_dt, tm);
//const auto &snia_feedback = hydro_pkg->Param<SNIAFeedback>("snia_feedback");
//snia_feedback.FeedbackSrcTerm(md, beta_dt, tm);

//ApplyClusterClips(md, tm, beta_dt);

const auto &stellar_feedback = hydro_pkg->Param<StellarFeedback>("stellar_feedback");
stellar_feedback.FeedbackSrcTerm(md, beta_dt, tm);

};
void ClusterSplitSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
const Real dt) {
const Real dt) {
auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");

const auto &stellar_feedback = hydro_pkg->Param<StellarFeedback>("stellar_feedback");
Expand Down Expand Up @@ -288,49 +295,50 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
made and a history output is not, then the mass/energy between the last
history output and the restart dump is lost
*/
std::string reduction_strs[] = {"stellar_mass", "added_dfloor_mass",
"removed_eceil_energy", "removed_vceil_energy",
"added_vAceil_mass"};
std::string reduction_strs[] = {"stellar_mass","added_dfloor_mass", "removed_eceil_energy",
"removed_vceil_energy", "added_vAceil_mass"};

// Add a param for each reduction, then add it as a summation reduction for
// history outputs
//Add a param for each reduction, then add it as a summation reduction for
//history outputs
auto hst_vars = hydro_pkg->Param<parthenon::HstVar_list>(parthenon::hist_param_key);

for (auto reduction_str : reduction_strs) {
hydro_pkg->AddParam(reduction_str, 0.0, true);
for( auto reduction_str : reduction_strs ) {
hydro_pkg->AddParam(reduction_str, 0.0, true);
hst_vars.emplace_back(parthenon::HistoryOutputVar(
parthenon::UserHistoryOperation::sum,
[reduction_str](MeshData<Real> *md) {
auto pmb = md->GetBlockData(0)->GetBlockPointer();
auto hydro_pkg = pmb->packages.Get("Hydro");
const Real reduction = hydro_pkg->Param<Real>(reduction_str);
// Reset the running count for this reduction between history outputs
hydro_pkg->UpdateParam(reduction_str, 0.0);
//Reset the running count for this reduction between history outputs
hydro_pkg->UpdateParam(reduction_str,0.0);
return reduction;
},
reduction_str));
}

// Add history reduction for total cold gas using stellar mass threshold
//Add history reduction for total cold gas using stellar mass threshold
const Real cold_thresh =
pin->GetOrAddReal("problem/cluster/reductions", "cold_temp_thresh", 0.0);
if (cold_thresh > 0) {
if( cold_thresh > 0){
hydro_pkg->AddParam("reduction_cold_threshold", cold_thresh);
hst_vars.emplace_back(parthenon::HistoryOutputVar(
parthenon::UserHistoryOperation::sum, LocalReduceColdGas, "cold_mass"));
parthenon::UserHistoryOperation::sum, LocalReduceColdGas,
"cold_mass"));
}
const Real agn_tracer_thresh =
pin->GetOrAddReal("problem/cluster/reductions", "agn_tracer_thresh", -1.0);
if (agn_tracer_thresh >= 0) {
auto mbar_over_kb = hydro_pkg->Param<Real>("mbar_over_kb");
PARTHENON_REQUIRE(
pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false),
"AGN Tracer must be enabled to reduce AGN tracer extent");
if( agn_tracer_thresh >= 0){
auto mbar_over_kb = hydro_pkg->Param<Real>("mbar_over_kb");
PARTHENON_REQUIRE(pin->GetOrAddBoolean("problem/cluster/agn_feedback", "enable_tracer", false),
"AGN Tracer must be enabled to reduce AGN tracer extent");
hydro_pkg->AddParam("reduction_agn_tracer_threshold", agn_tracer_thresh);
hst_vars.emplace_back(parthenon::HistoryOutputVar(
parthenon::UserHistoryOperation::max, LocalReduceAGNExtent, "agn_extent"));
parthenon::UserHistoryOperation::max, LocalReduceAGNExtent,
"agn_extent"));
}


hydro_pkg->UpdateParam(parthenon::hist_param_key, hst_vars);

/************************************************************
Expand Down Expand Up @@ -706,9 +714,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
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;
const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR);
const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR);
const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR);
auto v_norm = std::sqrt(v2_sum / (Lx * Ly * Lz) / (SQR(sigma_v)));

pmb->par_for(
Expand Down Expand Up @@ -784,9 +792,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
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;
const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR);
const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR);
const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR);
auto b_norm = std::sqrt(b2_sum / (Lx * Ly * Lz) / (SQR(sigma_b)));

pmb->par_for(
Expand Down
6 changes: 4 additions & 2 deletions src/pgen/cpaw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@

namespace cpaw {
using namespace parthenon::driver::prelude;
using namespace parthenon::package::prelude;

// Parameters which define initial solution -- made global so that they can be shared
// with functions A1,2,3 which compute vector potentials
Real den, pres, gm1, b_par, b_perp, v_perp, v_par;
Expand Down Expand Up @@ -213,8 +215,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
}

// write errors
std::fprintf(pfile, "%d %d", mesh->mesh_size.nx1, mesh->mesh_size.nx2);
std::fprintf(pfile, " %d %d %e", mesh->mesh_size.nx3, tm.ncycle, rms_err);
std::fprintf(pfile, "%d %d", mesh->mesh_size.nx(X1DIR), mesh->mesh_size.nx(X2DIR));
std::fprintf(pfile, " %d %d %e", mesh->mesh_size.nx(X3DIR), tm.ncycle, rms_err);
std::fprintf(pfile, " %e %e %e %e", err[IDN], err[IM1], err[IM2], err[IM3]);
std::fprintf(pfile, " %e", err[IEN]);
std::fprintf(pfile, " %e %e %e", err[IB1], err[IB2], err[IB3]);
Expand Down
6 changes: 3 additions & 3 deletions src/pgen/field_loop.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,13 +135,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
int iprob = pin->GetInteger("problem/field_loop", "iprob");
Real ang_2, cos_a2(0.0), sin_a2(0.0), lambda(0.0);

Real x1size = pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min;
Real x2size = pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min;
Real x1size = pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR);
Real x2size = pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR);

const bool two_d = pmb->pmy_mesh->ndim < 3;
// for 2D sim set x3size to zero so that v_z is 0 below
Real x3size =
two_d ? 0 : pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min;
two_d ? 0 : pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR);

// For (iprob=4) -- rotated cylinder in 3D -- set up rotation angle and wavelength
if (iprob == 4) {
Expand Down
11 changes: 6 additions & 5 deletions src/pgen/linear_wave.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@

namespace linear_wave {
using namespace parthenon::driver::prelude;
using namespace parthenon::package::prelude;

// TODO(pgrete) temp fix to address removal in Parthenon. Update when merging with MHD
constexpr int NWAVE = 5;
Expand Down Expand Up @@ -280,9 +281,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
if (parthenon::Globals::my_rank == 0) {
// normalize errors by number of cells
const auto mesh_size = mesh->mesh_size;
const auto vol = (mesh_size.x1max - mesh_size.x1min) *
(mesh_size.x2max - mesh_size.x2min) *
(mesh_size.x3max - mesh_size.x3min);
const auto vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) *
(mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) *
(mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR));
for (int i = 0; i < (NHYDRO + NFIELD); ++i)
l1_err[i] = l1_err[i] / vol;
// compute rms error
Expand Down Expand Up @@ -320,8 +321,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
}

// write errors
std::fprintf(pfile, "%d %d", mesh_size.nx1, mesh_size.nx2);
std::fprintf(pfile, " %d %d", mesh_size.nx3, tm.ncycle);
std::fprintf(pfile, "%d %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR));
std::fprintf(pfile, " %d %d", mesh_size.nx(X2DIR), tm.ncycle);
std::fprintf(pfile, " %e %e", rms_err, l1_err[IDN]);
std::fprintf(pfile, " %e %e %e", l1_err[IM1], l1_err[IM2], l1_err[IM3]);
std::fprintf(pfile, " %e", l1_err[IEN]);
Expand Down
11 changes: 6 additions & 5 deletions src/pgen/linear_wave_mhd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@

namespace linear_wave_mhd {
using namespace parthenon::driver::prelude;
using namespace parthenon::package::prelude;

constexpr int NMHDWAVE = 7;
// Parameters which define initial solution -- made global so that they can be shared
Expand Down Expand Up @@ -300,9 +301,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
if (parthenon::Globals::my_rank == 0) {
// normalize errors by number of cells
const auto mesh_size = mesh->mesh_size;
const auto vol = (mesh_size.x1max - mesh_size.x1min) *
(mesh_size.x2max - mesh_size.x2min) *
(mesh_size.x3max - mesh_size.x3min);
const auto vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) *
(mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) *
(mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR));
for (int i = 0; i < (NGLMMHD); ++i)
l1_err[i] = l1_err[i] / vol;
// compute rms error
Expand Down Expand Up @@ -342,8 +343,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
}

// write errors
std::fprintf(pfile, "%d %d", mesh_size.nx1, mesh_size.nx2);
std::fprintf(pfile, " %d %d", mesh_size.nx3, tm.ncycle);
std::fprintf(pfile, "%d %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR));
std::fprintf(pfile, " %d %d", mesh_size.nx(X3DIR), tm.ncycle);
std::fprintf(pfile, " %e %e", rms_err, l1_err[IDN]);
std::fprintf(pfile, " %e %e %e", l1_err[IM1], l1_err[IM2], l1_err[IM3]);
std::fprintf(pfile, " %e", l1_err[IEN]);
Expand Down
14 changes: 7 additions & 7 deletions src/pgen/turbulence.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,10 +217,10 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
const auto gm1 = pin->GetReal("hydro", "gamma") - 1.0;
const auto p0 = pin->GetReal("problem/turbulence", "p0");
const auto rho0 = pin->GetReal("problem/turbulence", "rho0");
const auto x3min = pmesh->mesh_size.x3min;
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;
const auto x3min = pmesh->mesh_size.xmin(X3DIR);
const auto Lx = pmesh->mesh_size.xmax(X1DIR) - pmesh->mesh_size.xmin(X1DIR);
const auto Ly = pmesh->mesh_size.xmax(X2DIR) - pmesh->mesh_size.xmin(X2DIR);
const auto Lz = pmesh->mesh_size.xmax(X3DIR) - pmesh->mesh_size.xmin(X3DIR);
const auto kz = 2.0 * M_PI / Lz;

// already pack data here to get easy access to coords in kernels
Expand Down Expand Up @@ -402,9 +402,9 @@ void Perturb(MeshData<Real> *md, const Real dt) {
MPI_SUM, MPI_COMM_WORLD));
#endif // MPI_PARALLEL

const auto Lx = pmb->pmy_mesh->mesh_size.x1max - pmb->pmy_mesh->mesh_size.x1min;
const auto Ly = pmb->pmy_mesh->mesh_size.x2max - pmb->pmy_mesh->mesh_size.x2min;
const auto Lz = pmb->pmy_mesh->mesh_size.x3max - pmb->pmy_mesh->mesh_size.x3min;
const auto Lx = pmb->pmy_mesh->mesh_size.xmax(X1DIR) - pmb->pmy_mesh->mesh_size.xmin(X1DIR);
const auto Ly = pmb->pmy_mesh->mesh_size.xmax(X2DIR) - pmb->pmy_mesh->mesh_size.xmin(X2DIR);
const auto Lz = pmb->pmy_mesh->mesh_size.xmax(X3DIR) - pmb->pmy_mesh->mesh_size.xmin(X3DIR);
const auto accel_rms = hydro_pkg->Param<Real>("turbulence/accel_rms");
auto norm = accel_rms / std::sqrt(sums[0] / (Lx * Ly * Lz));

Expand Down
24 changes: 12 additions & 12 deletions src/utils/few_modes_ft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,18 +106,18 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) {
"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;
const auto Lx1 = pm->mesh_size.xmax(X1DIR) - pm->mesh_size.xmin(X1DIR);
const auto Lx2 = pm->mesh_size.xmax(X2DIR) - pm->mesh_size.xmin(X2DIR);
const auto Lx3 = pm->mesh_size.xmax(X3DIR) - pm->mesh_size.xmin(X3DIR);

// 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);
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);
auto gnx1 = pm->mesh_size.nx(X1DIR) * std::pow(2, pmb->loc.level() - root_level);
auto gnx2 = pm->mesh_size.nx(X2DIR) * std::pow(2, pmb->loc.level() - root_level);
auto gnx3 = pm->mesh_size.nx(X3DIR) * std::pow(2, pmb->loc.level() - root_level);

// Restriction should also be easily fixed, just need to double check transforms and
// volume weighting everywhere
Expand All @@ -126,13 +126,13 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) {
"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 nx1 = pmb->block_size.nx(X1DIR);
const auto nx2 = pmb->block_size.nx(X2DIR);
const auto nx3 = pmb->block_size.nx(X3DIR);

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;
const auto gis = pmb->loc.lx1() * pmb->block_size.nx(X1DIR);
const auto gjs = pmb->loc.lx2() * pmb->block_size.nx(X2DIR);
const auto gks = pmb->loc.lx3() * pmb->block_size.nx(X3DIR);

// make local ref to capure in lambda
const auto num_modes = num_modes_;
Expand Down

0 comments on commit 6853b2d

Please sign in to comment.