diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 00000000..73a31d4b --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,23 @@ +# Changelog + +## Current develop (i.e., `main` branch) + +### Added (new features/APIs/variables/...) + +### Changed (changing behavior/API/variables/...) +- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15) + +### Fixed (not changing behavior/API/variables/...) + +### Infrastructure +- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Added `CHANGELOG.md` + +### Removed (removing behavior/API/varaibles/...) + +### Incompatibilities (i.e. breaking changes) +- [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15) + - Updated access to block dimension: `pmb->block_size.nx1` -> `pmb->block_size.nx(X1DIR)` (and similarly x2 and x3) + - Update access to mesh size: `pmesh->mesh_size.x1max` -> `pmesh->mesh_size.xmax(X1DIR)` (and similarly x2, x3, and min) + - Updated Parthenon `GradMinMod` signature for custom prolongation ops + - `GetBlockPointer` returns a raw pointer not a shared one (and updated interfaces to use raw pointers rather than shared ones) + diff --git a/external/parthenon b/external/parthenon index 5b6bb619..88ece809 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 5b6bb61906f7c278f9724ee9f38e79dee8707098 +Subproject commit 88ece809dc56ae16d695db8d9051433bff01feb1 diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 6d0e087c..e5402454 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -775,8 +775,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; @@ -848,7 +848,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/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp index 02d2b15f..4693e755 100644 --- a/src/hydro/prolongation/custom_ops.hpp +++ b/src/hydro/prolongation/custom_ops.hpp @@ -81,34 +81,40 @@ struct ProlongateCellMinModMultiD { Real dx1fm = 0; Real dx1fp = 0; + [[maybe_unused]] Real gx1m = 0, gx1p = 0; Real gx1c = 0; if constexpr (INCLUDE_X1) { Real dx1m, dx1p; GetGridSpacings<1, el>(coords, coarse_coords, cib, ib, i, fi, &dx1m, &dx1p, &dx1fm, &dx1fp); - gx1c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1), - coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p); + gx1c = + GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1), + coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p, gx1m, gx1p); } Real dx2fm = 0; Real dx2fp = 0; + [[maybe_unused]] Real gx2m = 0, gx2p = 0; Real gx2c = 0; if constexpr (INCLUDE_X2) { Real dx2m, dx2p; GetGridSpacings<2, el>(coords, coarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, &dx2fm, &dx2fp); - gx2c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i), - coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p); + gx2c = + GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i), + coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p, gx2m, gx2p); } Real dx3fm = 0; Real dx3fp = 0; + [[maybe_unused]] Real gx3m = 0, gx3p = 0; Real gx3c = 0; if constexpr (INCLUDE_X3) { Real dx3m, dx3p; GetGridSpacings<3, el>(coords, coarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, &dx3fm, &dx3fp); - gx3c = GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i), - coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p); + gx3c = + GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i), + coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p, gx3m, dx3p); } // Max. expected total difference. (dx#fm/p are positive by construction) diff --git a/src/pgen/cloud.cpp b/src/pgen/cloud.cpp index bc9fed25..b496dd2d 100644 --- a/src/pgen/cloud.cpp +++ b/src/pgen/cloud.cpp @@ -213,7 +213,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } void InflowWindX2(std::shared_ptr> &mbd, bool coarse) { - std::shared_ptr pmb = mbd->GetBlockPointer(); + auto pmb = mbd->GetBlockPointer(); auto cons = mbd->PackVariables(std::vector{"cons"}, coarse); // TODO(pgrete) Add par_for_bndry to Parthenon without requiring nb const auto nb = IndexRange{0, 0}; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index a2b10f7e..1e73d8b5 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -25,6 +25,7 @@ #include // c_str() // Parthenon headers +#include "Kokkos_MathematicalFunctions.hpp" #include "kokkos_abstraction.hpp" #include "mesh/domain.hpp" #include "mesh/mesh.hpp" @@ -347,6 +348,11 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddField("mach_sonic", m); // temperature hydro_pkg->AddField("temperature", m); + // radial velocity + hydro_pkg->AddField("v_r", m); + + // spherical theta + hydro_pkg->AddField("theta_sph", m); if (hydro_pkg->Param("enable_cooling") == Cooling::tabular) { // cooling time @@ -359,6 +365,9 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd // plasma beta hydro_pkg->AddField("plasma_beta", m); + + // plasma beta + hydro_pkg->AddField("B_mag", m); } /************************************************************ @@ -563,7 +572,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { ************************************************************/ const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); - magnetic_tower.AddInitialFieldToPotential(pmb.get(), a_kb, a_jb, a_ib, A); + magnetic_tower.AddInitialFieldToPotential(pmb, a_kb, a_jb, a_ib, A); /************************************************************ * Add dipole magnetic field to the magnetic potential @@ -666,7 +675,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { // Init phases on all blocks for (int b = 0; b < md->NumBlocks(); b++) { auto pmb = md->GetBlockData(b)->GetBlockPointer(); - few_modes_ft.SetPhases(pmb.get(), pin); + few_modes_ft.SetPhases(pmb, pin); } // As for t_corr in few_modes_ft, the choice for dt is // in principle arbitrary because the inital v_hat is 0 and the v_hat_new will contain @@ -705,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( @@ -734,7 +743,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { // Init phases on all blocks for (int b = 0; b < md->NumBlocks(); b++) { auto pmb = md->GetBlockData(b)->GetBlockPointer(); - few_modes_ft.SetPhases(pmb.get(), pin); + few_modes_ft.SetPhases(pmb, pin); } // As for t_corr in few_modes_ft, the choice for dt is // in principle arbitrary because the inital b_hat is 0 and the b_hat_new will contain @@ -783,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( @@ -818,6 +827,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto &entropy = data->Get("entropy").data; auto &mach_sonic = data->Get("mach_sonic").data; auto &temperature = data->Get("temperature").data; + auto &v_r = data->Get("v_r").data; + auto &theta_sph = data->Get("theta_sph").data; // for computing temperature from primitives auto units = pkg->Param("units"); @@ -844,8 +855,12 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { const Real x = coords.Xc<1>(i); const Real y = coords.Xc<2>(j); const Real z = coords.Xc<3>(k); - const Real r2 = SQR(x) + SQR(y) + SQR(z); - log10_radius(k, j, i) = 0.5 * std::log10(r2); + const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); + log10_radius(k, j, i) = std::log10(r); + + v_r(k, j, i) = ((v1 * x) + (v2 * y) + (v3 * z)) / r; + + theta_sph(k, j, i) = std::acos(z / r); // compute entropy const Real K = P / std::pow(rho / mbar, gam); @@ -884,6 +899,7 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { if (pkg->Param("fluid") == Fluid::glmmhd) { auto &plasma_beta = data->Get("plasma_beta").data; auto &mach_alfven = data->Get("mach_alfven").data; + auto &b_mag = data->Get("B_mag").data; pmb->par_for( "Cluster::UserWorkBeforeOutput::MHD", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -896,6 +912,8 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { const Real Bz = prim(IB3, k, j, i); const Real B2 = (SQR(Bx) + SQR(By) + SQR(Bz)); + b_mag(k, j, i) = Kokkos::sqrt(B2); + // compute Alfven mach number const Real v_A = std::sqrt(B2 / rho); const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS 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..45eeb60e 100644 --- a/src/pgen/field_loop.cpp +++ b/src/pgen/field_loop.cpp @@ -135,13 +135,16 @@ 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..2c9ff9a7 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,12 @@ 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_;