diff --git a/external/parthenon b/external/parthenon index 5b6bb619..60d54786 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 5b6bb61906f7c278f9724ee9f38e79dee8707098 +Subproject commit 60d5478693e26fa9b4836beb60ba7aeedfc608a4 diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index d22dddef..846adb9d 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -794,8 +794,8 @@ TaskStatus CalculateFluxes(std::shared_ptr> &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; @@ -867,7 +867,7 @@ TaskStatus CalculateFluxes(std::shared_ptr> &md) { parthenon::ScratchPad2D::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; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 96a781ba..de489177 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -59,8 +59,9 @@ using namespace parthenon::driver::prelude; using namespace parthenon::package::prelude; using utils::few_modes_ft::FewModesFT; + void ClusterUnsplitSrcTerm(MeshData *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("gravity_srcterm"); @@ -78,11 +79,17 @@ void ClusterUnsplitSrcTerm(MeshData *md, const parthenon::SimTime &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); + //const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); + //snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + + //ApplyClusterClips(md, tm, beta_dt); + + const auto &stellar_feedback = hydro_pkg->Param("stellar_feedback"); + stellar_feedback.FeedbackSrcTerm(md, beta_dt, tm); + }; void ClusterSplitSrcTerm(MeshData *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("stellar_feedback"); @@ -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::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 *md) { auto pmb = md->GetBlockData(0)->GetBlockPointer(); auto hydro_pkg = pmb->packages.Get("Hydro"); const Real reduction = hydro_pkg->Param(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("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("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); /************************************************************ @@ -706,9 +714,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *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( @@ -784,9 +792,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *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( diff --git a/src/pgen/cpaw.cpp b/src/pgen/cpaw.cpp index c269c378..0dcefacb 100644 --- a/src/pgen/cpaw.cpp +++ b/src/pgen/cpaw.cpp @@ -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; @@ -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]); diff --git a/src/pgen/field_loop.cpp b/src/pgen/field_loop.cpp index 66940d56..0ccdb69a 100644 --- a/src/pgen/field_loop.cpp +++ b/src/pgen/field_loop.cpp @@ -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) { diff --git a/src/pgen/linear_wave.cpp b/src/pgen/linear_wave.cpp index bc504b06..5480875e 100644 --- a/src/pgen/linear_wave.cpp +++ b/src/pgen/linear_wave.cpp @@ -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; @@ -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 @@ -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]); diff --git a/src/pgen/linear_wave_mhd.cpp b/src/pgen/linear_wave_mhd.cpp index 7f1b549b..67affbe7 100644 --- a/src/pgen/linear_wave_mhd.cpp +++ b/src/pgen/linear_wave_mhd.cpp @@ -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 @@ -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 @@ -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]); diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 58882f52..9d0e96be 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -217,10 +217,10 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *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 @@ -402,9 +402,9 @@ void Perturb(MeshData *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("turbulence/accel_rms"); auto norm = accel_rms / std::sqrt(sums[0] / (Lx * Ly * Lz)); diff --git a/src/utils/few_modes_ft.cpp b/src/utils/few_modes_ft.cpp index f477b211..e42edde2 100644 --- a/src/utils/few_modes_ft.cpp +++ b/src/utils/few_modes_ft.cpp @@ -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 @@ -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_;