diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 73d6dcc2..ac5ed0f6 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -288,7 +288,19 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { namespace cr = radmoment_cons; namespace ir = radmoment_internal; namespace impl = internal_variables; + + using parthenon::MakePackDescriptor; auto *pmesh = rc->GetMeshPointer(); + auto &resolved_pkgs = pmesh->resolved_packages; + const int ndim = pmesh->ndim; + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto v = desc.GetPack(); + const int nblocks = v.GetNBlocks(); + IndexRange ib = rc->GetBoundsI(domain); IndexRange jb = rc->GetBoundsJ(domain); IndexRange kb = rc->GetBoundsK(domain); @@ -306,39 +318,6 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { if (!enable_floors) return TaskStatus::complete; - const std::vector vars( - {p::density::name(), c::density::name(), p::velocity::name(), c::momentum::name(), - p::energy::name(), c::energy::name(), p::bfield::name(), p::ye::name(), - c::ye::name(), p::pressure::name(), p::temperature::name(), p::gamma1::name(), - pr::J::name(), pr::H::name(), cr::E::name(), cr::F::name(), - impl::cell_signal_speed::name(), impl::fail::name(), ir::tilPi::name()}); - - PackIndexMap imap; - auto v = rc->PackVariables(vars, imap); - - const int prho = imap[p::density::name()].first; - const int crho = imap[c::density::name()].first; - const int pvel_lo = imap[p::velocity::name()].first; - const int pvel_hi = imap[p::velocity::name()].second; - const int cmom_lo = imap[c::momentum::name()].first; - const int cmom_hi = imap[c::momentum::name()].second; - const int peng = imap[p::energy::name()].first; - const int ceng = imap[c::energy::name()].first; - const int prs = imap[p::pressure::name()].first; - const int tmp = imap[p::temperature::name()].first; - const int gm1 = imap[p::gamma1::name()].first; - const int slo = imap[impl::cell_signal_speed::name()].first; - const int shi = imap[impl::cell_signal_speed::name()].second; - const int pb_lo = imap[p::bfield::name()].first; - const int pb_hi = imap[p::bfield::name()].second; - int pye = imap[p::ye::name()].second; // negative if not present - int cye = imap[c::ye::name()].second; - auto idx_J = imap.GetFlatIdx(pr::J::name(), false); - auto idx_H = imap.GetFlatIdx(pr::H::name(), false); - auto idx_E = imap.GetFlatIdx(cr::E::name(), false); - auto idx_F = imap.GetFlatIdx(cr::F::name(), false); - auto iTilPi = imap.GetFlatIdx(ir::tilPi::name(), false); - const int num_species = enable_rad_floors ? rad_pkg->Param("num_species") : 0; auto eos = eos_pkg->Param("d.EOS"); @@ -354,16 +333,16 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { c2p_floor_scale_fac, c2p_fail_on_floors, c2p_fail_on_ceilings); - Coordinates_t coords = rc->GetParentPointer()->coords; - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "ApplyFloors", DevExecSpace(), 0, v.GetDim(5) - 1, kb.s, kb.e, + DEFAULT_LOOP_PATTERN, "ApplyFloors", DevExecSpace(), 0, nblocks - 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) { + auto &coords = v.GetCoordinates(b); double eos_lambda[2]; // used for stellarcollapse eos and // other EOS's that require root // finding. - eos_lambda[1] = std::log10(v(b, tmp, k, j, i)); // use last temp as initial guess + eos_lambda[1] = + std::log10(v(b, p::temperature(), k, j, i)); // use last temp as initial guess double rho_floor, sie_floor; bounds.GetFloors(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), @@ -400,10 +379,11 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { if (enable_mhd_floors) { Real Bsq = 0.0; Real Bdotv = 0.0; - const Real vp[3] = {v(b, pvel_lo, k, j, i), v(b, pvel_lo + 1, k, j, i), - v(b, pvel_lo + 2, k, j, i)}; - const Real bp[3] = {v(b, pb_lo, k, j, i), v(pb_lo + 1, k, j, i), - v(pb_lo + 2, k, j, i)}; + const Real vp[3] = {v(b, p::velocity(0), k, j, i), + v(b, p::velocity(1), k, j, i), + v(b, p::velocity(2), k, j, i)}; + const Real bp[3] = {v(b, p::bfield(0), k, j, i), v(p::bfield(1), k, j, i), + v(p::bfield(2), k, j, i)}; const Real W = phoebus::GetLorentzFactor(vp, gcov); const Real iW = 1.0 / W; SPACELOOP2(ii, jj) { @@ -417,11 +397,12 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { u_floor_max = std::max(u_floor_max, bsq / bsqou_max); rho_floor_max = std::max( - rho_floor_max, std::max(v(b, peng, k, j, i), u_floor_max) / e_max); + rho_floor_max, + std::max(v(b, p::energy(), k, j, i), u_floor_max) / e_max); } - Real drho = rho_floor_max - v(b, prho, k, j, i); - Real du = u_floor_max - v(b, peng, k, j, i); + Real drho = rho_floor_max - v(b, p::density(), k, j, i); + Real du = u_floor_max - v(b, p::energy(), k, j, i); if (drho > 0. || du > 0.) { floor_applied = true; drho = std::max(drho, du / sie_floor); @@ -430,8 +411,8 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { Real dcrho, dS[3], dBcons[3], dtau, dyecons; Real bp[3] = {0}; - if (pb_hi > 0) { - SPACELOOP(ii) { bp[ii] = v(b, pb_lo + ii, k, j, i); } + if (v.Contains(b,p::bfield())) { + SPACELOOP(ii) { bp[ii] = v(b, p::bfield(ii), k, j, i); } } if (floor_applied) { @@ -448,90 +429,93 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { dyecons); // Update cons vars (not B field) - v(b, crho, k, j, i) += dcrho; - SPACELOOP(ii) { v(b, cmom_lo + ii, k, j, i) += dS[ii]; } - v(b, ceng, k, j, i) += dtau; - if (pye > 0) { - v(b, cye, k, j, i) += dyecons; + v(b, c::density(), k, j, i) += dcrho; + SPACELOOP(ii) { v(b, c::momentum(ii), k, j, i) += dS[ii]; } + v(b, c::energy(), k, j, i) += dtau; + if (v.Contains(b, p::ye())) { + v(b, c::ye(), k, j, i) += dyecons; } // fluid c2p auto status = invert(geom, eos, coords, k, j, i); if (status == con2prim_robust::ConToPrimStatus::failure) { // If fluid c2p fails, set to floors - v(b, prho, k, j, i) = drho; - SPACELOOP(ii) { v(b, pvel_lo + ii, k, j, i) = vp_normalobs[ii]; } - v(b, peng, k, j, i) = du; - if (pye > 0) { - v(b, pye, k, j, i) = ye_prim_default; + v(b, p::density(), k, j, i) = drho; + SPACELOOP(ii) { v(b, p::velocity(ii), k, j, i) = vp_normalobs[ii]; } + v(b, p::energy(), k, j, i) = du; + if (v.Contains(b, p::ye())) { + v(b, p::ye(), k, j, i) = ye_prim_default; } // Update auxiliary primitives - if (pye > 0) { - eos_lambda[0] = v(b, pye, k, j, i); + if (v.Contains(b, p::ye())) { + eos_lambda[0] = v(b, p::ye(), k, j, i); } - v(b, tmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy( - v(b, prho, k, j, i), v(b, peng, k, j, i) / v(b, prho, k, j, i), - eos_lambda); + v(b, p::temperature(), k, j, i) = eos.TemperatureFromDensityInternalEnergy( + v(b, p::density(), k, j, i), + v(b, p::energy(), k, j, i) / v(b, p::density(), k, j, i), eos_lambda); // Update cons vars (not B field) prim2con::p2c(drho, vp_normalobs, bp, du, ye_prim_default, dprs, dgm1, gcov, - gammacon, betacon, alpha, sdetgam, v(b, crho, k, j, i), dS, - dBcons, v(b, ceng, k, j, i), dyecons); - SPACELOOP(ii) { v(b, cmom_lo + ii, k, j, i) = dS[ii]; } - if (pye > 0) { - v(b, cye, k, j, i) = dyecons; + gammacon, betacon, alpha, sdetgam, v(b, c::density(), k, j, i), + dS, dBcons, v(b, c::energy(), k, j, i), dyecons); + SPACELOOP(ii) { v(b, c::momentum(ii), k, j, i) = dS[ii]; } + if (v.Contains(b, p::ye())) { + v(b, c::ye(), k, j, i) = dyecons; } } } // Fluid ceilings? - Real vpcon[3] = {v(b, pvel_lo, k, j, i), v(b, pvel_lo + 1, k, j, i), - v(b, pvel_lo + 2, k, j, i)}; + Real vpcon[3] = {v(b, p::velocity(0), k, j, i), v(b, p::velocity(1), k, j, i), + v(b, p::velocity(2), k, j, i)}; const Real W = phoebus::GetLorentzFactor(vpcon, gcov); - if (W > gamma_max || v(b, peng, k, j, i) / v(b, prho, k, j, i) > e_max) { + if (W > gamma_max || + v(b, p::energy(), k, j, i) / v(b, p::density(), k, j, i) > e_max) { ceiling_applied = true; } if (ceiling_applied) { const Real rescale = std::sqrt(gamma_max * gamma_max - 1.) / (W * W - 1.); SPACELOOP(ii) { vpcon[ii] *= rescale; } - SPACELOOP(ii) { v(b, pvel_lo + ii, k, j, i) = vpcon[ii]; } + SPACELOOP(ii) { v(b, p::velocity(ii), k, j, i) = vpcon[ii]; } Real ye = 0.; - if (pye > 0) { - ye = v(b, pye, k, j, i); + if (v.Contains(b, p::ye())) { + ye = v(b, p::ye(), k, j, i); } - prim2con::p2c(v(b, prho, k, j, i), vpcon, bp, v(b, peng, k, j, i), ye, - v(b, prs, k, j, i), v(b, gm1, k, j, i), gcov, gammacon, betacon, - alpha, sdetgam, v(b, crho, k, j, i), dS, dBcons, - v(b, ceng, k, j, i), dyecons); - SPACELOOP(ii) { v(b, cmom_lo + ii, k, j, i) = dS[ii]; } - if (pye > 0) { - v(b, cye, k, j, i) = dyecons; + prim2con::p2c(v(b, p::density(), k, j, i), vpcon, bp, v(b, p::energy(), k, j, i), + ye, v(b, p::pressure(), k, j, i), v(b, p::gamma1(), k, j, i), + gcov, gammacon, betacon, alpha, sdetgam, + v(b, c::density(), k, j, i), dS, dBcons, v(b, c::energy(), k, j, i), + dyecons); + SPACELOOP(ii) { v(b, c::momentum(ii), k, j, i) = dS[ii]; } + if (v.Contains(b, p::ye())) { + v(b, c::ye(), k, j, i) = dyecons; } } if (floor_applied || ceiling_applied) { // Update derived prims - if (pye > 0) eos_lambda[0] = v(b, pye, k, j, i); - v(b, tmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy( - v(b, prho, k, j, i), ratio(v(b, peng, k, j, i), v(b, prho, k, j, i)), - eos_lambda); - v(b, prs, k, j, i) = eos.PressureFromDensityTemperature( - v(b, prho, k, j, i), v(b, tmp, k, j, i), eos_lambda); - v(b, gm1, k, j, i) = - ratio(eos.BulkModulusFromDensityTemperature(v(b, prho, k, j, i), - v(b, tmp, k, j, i), eos_lambda), - v(b, prs, k, j, i)); + if (v.Contains(b, p::ye())) eos_lambda[0] = v(b, p::ye(), k, j, i); + v(b, p::temperature(), k, j, i) = eos.TemperatureFromDensityInternalEnergy( + v(b, p::density(), k, j, i), + ratio(v(b, p::energy(), k, j, i), v(b, p::density(), k, j, i)), eos_lambda); + v(b, p::pressure(), k, j, i) = eos.PressureFromDensityTemperature( + v(b, p::density(), k, j, i), v(b, p::temperature, k, j, i), eos_lambda); + v(b, p::gamma1(), k, j, i) = + ratio(eos.BulkModulusFromDensityTemperature(v(b, p::density(), k, j, i), + v(b, p::temperature(), k, j, i), + eos_lambda), + v(b, p::pressure(), k, j, i)); } if (rad_active) { - Vec con_vp{v(b, pvel_lo, k, j, i), v(b, pvel_lo + 1, k, j, i), - v(b, pvel_lo + 2, k, j, i)}; + Vec con_vp{v(b, p::velocity(0), k, j, i), v(b, p::velocity(1), k, j, i), + v(b, p::velocity(2), k, j, i)}; const Real W = phoebus::GetLorentzFactor(con_vp.data, gcov); - Vec con_v{v(b, pvel_lo, k, j, i) / W, v(b, pvel_lo + 1, k, j, i) / W, - v(b, pvel_lo + 2, k, j, i) / W}; + Vec con_v{v(b, p::velocity(0), k, j, i) / W, v(b, p::velocity(1), k, j, i) / W, + v(b, p::velocity(2), k, j, i) / W}; typename CLOSURE::LocalGeometryType g(geom, CellLocation::Cent, b, k, j, i); Real E; Real J; @@ -546,13 +530,13 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { // rad c2p for (int ispec = 0; ispec < num_species; ++ispec) { - E = v(b, idx_E(ispec), k, j, i) / sdetgam; - SPACELOOP(ii) { cov_F(ii) = v(b, idx_F(ispec, ii), k, j, i) / sdetgam; } + E = v(b, cr::E(ispec), k, j, i) / sdetgam; + SPACELOOP(ii) { cov_F(ii) = v(b, cr::F(ispec, ii), k, j, i) / sdetgam; } // We need the real conTilPi - if (iTilPi.IsValid()) { + if (v.Contains(b, ir::tilPi())) { SPACELOOP2(ii, jj) { - con_TilPi(ii, jj) = v(b, iTilPi(ispec, ii, jj), k, j, i); + con_TilPi(ii, jj) = v(b, ir::tilPi(ispec, ii, jj), k, j, i); } } else { // TODO(BRR) don't be lazy and actually retrieve these @@ -562,8 +546,8 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { } c.Con2Prim(E, cov_F, con_TilPi, &J, &cov_H); - v(b, idx_J(ispec), k, j, i) = J; - SPACELOOP(ii) { v(b, idx_H(ispec, ii), k, j, i) = cov_H(ii) / J; } + v(b, pr::J(ispec), k, j, i) = J; + SPACELOOP(ii) { v(b, pr::H(ispec, ii), k, j, i) = cov_H(ii) / J; } } } @@ -572,7 +556,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { CLOSURE c_iso(con_v_normalobs, &g); Tens2 con_tilPi_iso{0}; for (int ispec = 0; ispec < num_species; ++ispec) { - Real dJ = J_floor - v(b, idx_J(ispec), k, j, i); + Real dJ = J_floor - v(b, pr::J(ispec), k, j, i); if (dJ > 0.) { constexpr bool update_cons_vars = true; // false; @@ -588,16 +572,16 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { PARTHENON_DEBUG_REQUIRE(!std::isnan(cov_H(2)), "bad"); c_iso.Prim2Con(J, cov_H, con_tilPi_iso, &E, &cov_F); - E += v(b, idx_E(ispec), k, j, i) / sdetgam; - SPACELOOP(ii) { cov_F(ii) += v(b, idx_F(ispec, ii), k, j, i) / sdetgam; } + E += v(b, cr::E(ispec), k, j, i) / sdetgam; + SPACELOOP(ii) { cov_F(ii) += v(b, cr::F(ispec, ii), k, j, i) / sdetgam; } - v(b, idx_E(ispec), k, j, i) = E * sdetgam; - SPACELOOP(ii) { v(b, idx_F(ispec, ii), k, j, i) = cov_F(ii) * sdetgam; } + v(b, cr::E(ispec), k, j, i) = E * sdetgam; + SPACELOOP(ii) { v(b, cr::F(ispec, ii), k, j, i) = cov_F(ii) * sdetgam; } // We need the real conTilPi - if (iTilPi.IsValid()) { + if (v.Contains(b, ir::tilPi)) { SPACELOOP2(ii, jj) { - con_TilPi(ii, jj) = v(b, iTilPi(ispec, ii, jj), k, j, i); + con_TilPi(ii, jj) = v(b, ir::tilPi(ispec, ii, jj), k, j, i); } } else { // TODO(BRR) don't be lazy and actually retrieve these @@ -608,18 +592,18 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { c.Con2Prim(E, cov_F, con_TilPi, &J, &cov_H); - v(b, idx_J(ispec), k, j, i) = J; - SPACELOOP(ii) { v(b, idx_H(ispec, ii), k, j, i) = cov_H(ii) / J; } + v(b, pr::J(ispec), k, j, i) = J; + SPACELOOP(ii) { v(b, pr::H(ispec, ii), k, j, i) = cov_H(ii) / J; } } else { - v(b, idx_J(ispec), k, j, i) += dJ; + v(b, pr::J(ispec), k, j, i) += dJ; - J = v(b, idx_J(ispec), k, j, i); - SPACELOOP(ii) { cov_H(ii) = v(b, idx_H(ispec, ii), k, j, i) * J; } + J = v(b, pr::J(ispec), k, j, i); + SPACELOOP(ii) { cov_H(ii) = v(b, pr::H(ispec, ii), k, j, i) * J; } // We need the real conTilPi - if (iTilPi.IsValid()) { + if (v.Contains(b, ir::tilPi())) { SPACELOOP2(ii, jj) { - con_TilPi(ii, jj) = v(b, iTilPi(ispec, ii, jj), k, j, i); + con_TilPi(ii, jj) = v(b, ir::tilPi(ispec, ii, jj), k, j, i); } } else { c.GetConTilPiFromPrim(J, cov_H, &con_TilPi); @@ -631,29 +615,29 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { PARTHENON_DEBUG_REQUIRE(!std::isnan(cov_H(2)), "bad"); c.Prim2Con(J, cov_H, con_TilPi, &E, &cov_F); - v(b, idx_E(ispec), k, j, i) = E * sdetgam; - SPACELOOP(ii) { v(b, idx_F(ispec, ii), k, j, i) = cov_F(ii) * sdetgam; } + v(b, cr::E(ispec), k, j, i) = E * sdetgam; + SPACELOOP(ii) { v(b, cr::F(ispec, ii), k, j, i) = cov_F(ii) * sdetgam; } } } - Vec cov_H{v(b, idx_H(ispec, 0), k, j, i), v(b, idx_H(ispec, 1), k, j, i), - v(b, idx_H(ispec, 2), k, j, i)}; + Vec cov_H{v(b, pr::H(ispec, 0), k, j, i), v(b, pr::H(ispec, 1), k, j, i), + v(b, pr::H(ispec, 2), k, j, i)}; const Real xi = std::sqrt(g.contractCov3Vectors(cov_H, cov_H) - std::pow(g.contractConCov3Vectors(con_v, cov_H), 2)); if (xi > xi_max) { - J = v(b, idx_J(ispec), k, j, i); + J = v(b, pr::J(ispec), k, j, i); SPACELOOP(ii) { - cov_H(ii) = (xi_max / xi) * v(b, idx_H(ispec, ii), k, j, i) * J; - v(b, idx_H(ispec, ii), k, j, i) = cov_H(ii) / J; + cov_H(ii) = (xi_max / xi) * v(b, pr::H(ispec, ii), k, j, i) * J; + v(b, pr::H(ispec, ii), k, j, i) = cov_H(ii) / J; } // We need the real conTilPi - if (iTilPi.IsValid()) { + if (v.Contains(b, ir::tilPi())) { SPACELOOP2(ii, jj) { - con_TilPi(ii, jj) = v(b, iTilPi(ispec, ii, jj), k, j, i); + con_TilPi(ii, jj) = v(b, ir::tilPi(ispec, ii, jj), k, j, i); } } else { c.GetConTilPiFromPrim(J, cov_H, &con_TilPi); @@ -663,8 +647,8 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { PARTHENON_DEBUG_REQUIRE(!std::isnan(cov_H(1)), "bad"); PARTHENON_DEBUG_REQUIRE(!std::isnan(cov_H(2)), "bad"); c.Prim2Con(J, cov_H, con_TilPi, &E, &cov_F); - v(b, idx_E(ispec), k, j, i) = E * sdetgam; - SPACELOOP(ii) { v(b, idx_F(ispec, ii), k, j, i) = cov_F(ii) * sdetgam; } + v(b, cr::E(ispec), k, j, i) = E * sdetgam; + SPACELOOP(ii) { v(b, cr::F(ispec, ii), k, j, i) = cov_F(ii) * sdetgam; } } } } @@ -703,7 +687,7 @@ TaskStatus ApplyFloors(T *rc) { return TaskStatus::fail; } -template TaskStatus ApplyFloors>(MeshBlockData *rc); +template TaskStatus ApplyFloors>(MeshData *rc); TaskStatus FixFluxes(MeshBlockData *rc) { using parthenon::BoundaryFace; diff --git a/src/fixup/fixup_c2p.cpp b/src/fixup/fixup_c2p.cpp index eb9bf5c2..0b8e988d 100644 --- a/src/fixup/fixup_c2p.cpp +++ b/src/fixup/fixup_c2p.cpp @@ -49,71 +49,27 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { namespace ir = radmoment_internal; namespace pr = radmoment_prim; namespace cr = radmoment_cons; + using parthenon::MakePackDescriptor; Mesh *pmesh = rc->GetMeshPointer(); IndexRange ib = rc->GetBoundsI(IndexDomain::interior); IndexRange jb = rc->GetBoundsJ(IndexDomain::interior); IndexRange kb = rc->GetBoundsK(IndexDomain::interior); + auto &resolved_pkgs = pmesh->resolved_packages; + const int ndim = pmesh->ndim; + static auto desc = + MakePackDescriptor( + resolved_pkgs.get()); + auto v = desc.GetPack(); + const int nblocks = v.GetNBlocks(); StateDescriptor *fix_pkg = pmesh->packages.Get("fixup").get(); StateDescriptor *fluid_pkg = pmesh->packages.Get("fluid").get(); StateDescriptor *rad_pkg = pmesh->packages.Get("radiation").get(); StateDescriptor *eos_pkg = pmesh->packages.Get("eos").get(); - const std::vector vars({p::density::name(), - c::density::name(), - p::velocity::name(), - c::momentum::name(), - p::energy::name(), - c::energy::name(), - p::bfield::name(), - p::ye::name(), - c::ye::name(), - p::pressure::name(), - p::temperature::name(), - p::gamma1::name(), - impl::cell_signal_speed::name(), - impl::fail::name(), - ir::c2pfail::name(), - ir::tilPi::name(), - pr::J::name(), - pr::H::name(), - cr::E::name(), - cr::F::name(), - ir::xi::name(), - ir::phi::name()}); - - PackIndexMap imap; - auto v = rc->PackVariables(vars, imap); - - const int prho = imap[p::density::name()].first; - const int crho = imap[c::density::name()].first; - const int pvel_lo = imap[p::velocity::name()].first; - const int pvel_hi = imap[p::velocity::name()].second; - const int cmom_lo = imap[c::momentum::name()].first; - const int cmom_hi = imap[c::momentum::name()].second; - const int peng = imap[p::energy::name()].first; - const int ceng = imap[c::energy::name()].first; - const int prs = imap[p::pressure::name()].first; - const int tmp = imap[p::temperature::name()].first; - const int gm1 = imap[p::gamma1::name()].first; - const int slo = imap[impl::cell_signal_speed::name()].first; - const int shi = imap[impl::cell_signal_speed::name()].second; - const int pb_lo = imap[p::bfield::name()].first; - const int pb_hi = imap[p::bfield::name()].second; - int pye = imap[p::ye::name()].second; // negative if not present - int cye = imap[c::ye::name()].second; - - int ifail = imap[impl::fail::name()].first; - int irfail = imap[ir::c2pfail::name()].first; - - auto idx_E = imap.GetFlatIdx(cr::E::name(), false); - auto idx_F = imap.GetFlatIdx(cr::F::name(), false); - auto idx_J = imap.GetFlatIdx(pr::J::name(), false); - auto idx_H = imap.GetFlatIdx(pr::H::name(), false); - auto iTilPi = imap.GetFlatIdx(ir::tilPi::name(), false); - auto iXi = imap.GetFlatIdx(ir::xi::name(), false); - auto iPhi = imap.GetFlatIdx(ir::phi::name(), false); - const bool rad_active = rad_pkg->Param("active"); const int num_species = rad_active ? rad_pkg->Param("num_species") : 0; @@ -126,9 +82,9 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { int nfail_total; parthenon::par_reduce( parthenon::loop_pattern_mdrange_tag, "ConToPrim::Solve fixup failures", - DevExecSpace(), 0, v.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DevExecSpace(), 0, nblocks - 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, int &nf) { - if (v(b, ifail, k, j, i) == con2prim_robust::FailFlags::fail) { + if (v(b, impl::fail(), k, j, i) == con2prim_robust::FailFlags::fail) { nf++; } }, @@ -140,9 +96,9 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { nfail_total = 0; parthenon::par_reduce( parthenon::loop_pattern_mdrange_tag, "Rad ConToPrim::Solve fixup failures", - DevExecSpace(), 0, v.GetDim(5) - 1, kbi.s, kbi.e, jbi.s, jbi.e, ibi.s, ibi.e, + DevExecSpace(), 0, nblocks - 1, kbi.s, kbi.e, jbi.s, jbi.e, ibi.s, ibi.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, int &nf) { - if (v(b, ifail, k, j, i) == con2prim_robust::FailFlags::fail) { + if (v(b, impl::fail(), k, j, i) == con2prim_robust::FailFlags::fail) { nf++; } }, @@ -150,14 +106,10 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { printf("total interior nfail: %i\n", nfail_total); } - const int ndim = pmesh->ndim; - auto eos = eos_pkg->Param("d.EOS"); auto geom = Geometry::GetCoordinateSystem(rc); auto bounds = fix_pkg->Param("bounds"); - Coordinates_t coords = rc->GetParentPointer()->coords; - auto fluid_c2p_failure_strategy = fix_pkg->Param("fluid_c2p_failure_strategy"); const auto c2p_failure_force_fixup_both = @@ -167,22 +119,23 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { "As currently implemented this is a race condition!"); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "ConToPrim::Solve fixup", DevExecSpace(), 0, v.GetDim(5) - 1, + DEFAULT_LOOP_PATTERN, "ConToPrim::Solve fixup", DevExecSpace(), 0, nblocks - 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) { + auto &coords = v.GetCoordinates(b); Real eos_lambda[2]; // use last temp as initial guess eos_lambda[0] = 0.5; - eos_lambda[1] = std::log10(v(b, tmp, k, j, i)); + eos_lambda[1] = std::log10(v(b, p::temperature(), k, j, i)); Real gamma_max, e_max; bounds.GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i), coords.Xc<3>(k, j, i), gamma_max, e_max); if (c2p_failure_force_fixup_both && rad_active) { - if (v(b, ifail, k, j, i) == con2prim_robust::FailFlags::fail || - v(b, irfail, k, j, i) == radiation::FailFlags::fail) { - v(b, ifail, k, j, i) = con2prim_robust::FailFlags::fail; - v(b, irfail, k, j, i) = radiation::FailFlags::fail; + if (v(b, impl::fail(), k, j, i) == con2prim_robust::FailFlags::fail || + v(b, ir::c2pfail(), k, j, i) == radiation::FailFlags::fail) { + v(b, impl::fail(), k, j, i) = con2prim_robust::FailFlags::fail; + v(b, ir::c2pfail(), k, j, i) = radiation::FailFlags::fail; } } @@ -193,9 +146,9 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { auto fail = [&](const int k, const int j, const int i) { if (c2p_failure_force_fixup_both) { - return v(b, ifail, k, j, i) * v(b, irfail, k, j, i); + return v(b, impl::fail(), k, j, i) * v(b, ir::c2pfail(), k, j, i); } else { - return v(b, ifail, k, j, i); + return v(b, impl::fail(), k, j, i); } }; auto fixup = [&](const int iv, const Real inv_mask_sum) { @@ -213,11 +166,15 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { }; // When using IndexDomain::entire, can't stencil from e.g. i = 0 because 0 - 1 = // -1 is a segfault - if (v(b, ifail, k, j, i) == con2prim_robust::FailFlags::fail) { + if (v(b, impl::fail(), k, j, i) == con2prim_robust::FailFlags::fail) { Real num_valid = 0; - num_valid = v(b, ifail, k, j, i - 1) + v(b, ifail, k, j, i + 1); - if (ndim > 1) num_valid += v(b, ifail, k, j - 1, i) + v(b, ifail, k, j + 1, i); - if (ndim == 3) num_valid += v(b, ifail, k - 1, j, i) + v(b, ifail, k + 1, j, i); + num_valid = v(b, impl::fail(), k, j, i - 1) + v(b, impl::fail(), k, j, i + 1); + if (ndim > 1) + num_valid += + v(b, impl::fail(), k, j - 1, i) + v(b, impl::fail(), k, j + 1, i); + if (ndim == 3) + num_valid += + v(b, impl::fail(), k - 1, j, i) + v(b, impl::fail(), k + 1, j, i); // if (num_valid > 0.5 && // fluid_c2p_failure_strategy == FAILURE_STRATEGY::interpolate && i > ib.s && @@ -225,32 +182,32 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { if (num_valid > 0.5 && fluid_c2p_failure_strategy == FAILURE_STRATEGY::interpolate) { const Real norm = 1.0 / num_valid; - v(b, prho, k, j, i) = fixup(prho, norm); - for (int pv = pvel_lo; pv <= pvel_hi; pv++) { - v(b, pv, k, j, i) = fixup(pv, norm); + v(b, p::density(), k, j, i) = fixup(p::density(), norm); + for (int pv = 0; pv <= 2; pv++) { + v(b, p::velocity(pv), k, j, i) = fixup(p::velocity(pv), norm); } - v(b, peng, k, j, i) = fixup(peng, norm); + v(b, p::energy(), k, j, i) = fixup(p::energy(), norm); - if (pye > 0) v(b, pye, k, j, i) = fixup(pye, norm); + if (v.Contains(b, p::ye())) v(b, p::ye(), k, j, i) = fixup(p::ye(), norm); - v(b, prho, k, j, i) = - std::max(v(b, prho, k, j, i), 100. * robust::SMALL()); - v(b, peng, k, j, i) = - std::max(v(b, peng, k, j, i), 100. * robust::SMALL()); + v(b, p::density(), k, j, i) = + std::max(v(b, p::density(), k, j, i), 100. * robust::SMALL()); + v(b, p::energy(), k, j, i) = + std::max(v(b, p::energy(), k, j, i), 100. * robust::SMALL()); } else { // No valid neighbors; set fluid mass/energy to near-zero and set primitive // velocities to zero - v(b, prho, k, j, i) = 100. * robust::SMALL(); - v(b, peng, k, j, i) = 100. * robust::SMALL(); + v(b, p::density(), k, j, i) = 100. * robust::SMALL(); + v(b, p::energy(), k, j, i) = 100. * robust::SMALL(); // Safe value for ye - if (pye > 0) { - v(b, pye, k, j, i) = 0.5; + if (v.Contains(b, p::ye())) { + v(b, p::ye(), k, j, i) = 0.5; } // Zero primitive velocities - SPACELOOP(ii) { v(b, pvel_lo + ii, k, j, i) = 0.; } + SPACELOOP(ii) { v(b, p::velocity(ii), k, j, i) = 0.; } } const Real sdetgam = geom.DetGamma(CellLocation::Cent, k, j, i); @@ -263,74 +220,76 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { geom.MetricInverse(CellLocation::Cent, k, j, i, gcon); // Clamp velocity now (for rad inversion) - Real vpcon[3] = {v(b, pvel_lo, k, j, i), v(b, pvel_lo + 1, k, j, i), - v(b, pvel_lo + 2, k, j, i)}; + Real vpcon[3] = {v(b, p::velocity(0), k, j, i), v(b, p::velocity(1), k, j, i), + v(b, p::velocity(2), k, j, i)}; Real W = phoebus::GetLorentzFactor(vpcon, gcov); if (W > gamma_max) { const Real rescale = std::sqrt(gamma_max * gamma_max - 1.) / (W * W - 1.); SPACELOOP(ii) { vpcon[ii] *= rescale; } - SPACELOOP(ii) { v(b, pvel_lo + ii, k, j, i) = vpcon[ii]; } + SPACELOOP(ii) { v(b, p::velocity(ii), k, j, i) = vpcon[ii]; } W = gamma_max; } // Update dependent primitives - if (pye > 0) eos_lambda[0] = v(b, pye, k, j, i); - v(b, tmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy( - v(b, prho, k, j, i), ratio(v(b, peng, k, j, i), v(b, prho, k, j, i)), - eos_lambda); - v(b, prs, k, j, i) = eos.PressureFromDensityTemperature( - v(b, prho, k, j, i), v(b, tmp, k, j, i), eos_lambda); - v(b, gm1, k, j, i) = - ratio(eos.BulkModulusFromDensityTemperature(v(b, prho, k, j, i), - v(b, tmp, k, j, i), eos_lambda), - v(b, prs, k, j, i)); + if (v.Contains(b, p::ye())) eos_lambda[0] = v(b, p::ye(), k, j, i); + v(b, p::temperature(), k, j, i) = eos.TemperatureFromDensityInternalEnergy( + v(b, p::density(), k, j, i), + ratio(v(b, p::energy(), k, j, i), v(b, p::density(), k, j, i)), eos_lambda); + v(b, p::pressure(), k, j, i) = eos.PressureFromDensityTemperature( + v(b, p::density(), k, j, i), v(b, p::temperature(), k, j, i), eos_lambda); + v(b, p::gamma1(), k, j, i) = + ratio(eos.BulkModulusFromDensityTemperature(v(b, p::density(), k, j, i), + v(b, p::temperature(), k, j, i), + eos_lambda), + v(b, p::pressure(), k, j, i)); // Update conserved variables Real S[3]; Real bcons[3]; Real bp[3] = {0.0, 0.0, 0.0}; - if (pb_hi > 0) { - bp[0] = v(b, pb_lo, k, j, i); - bp[1] = v(b, pb_lo + 1, k, j, i); - bp[2] = v(b, pb_hi, k, j, i); + if (v.Contains(b, p::bfield(2))) { + bp[0] = v(b, p::bfield(0), k, j, i); + bp[1] = v(b, p::bfield(1), k, j, i); + bp[2] = v(b, p::bfield(2), k, j, i); } Real ye_cons; Real ye_prim = 0.5; - if (pye > 0) { - ye_prim = v(b, pye, k, j, i); + if (v.Contains(b, p::ye())) { + ye_prim = v(b, p::ye(), k, j, i); } Real sig[3]; - prim2con::p2c(v(b, prho, k, j, i), vpcon, bp, v(b, peng, k, j, i), ye_prim, - v(b, prs, k, j, i), v(b, gm1, k, j, i), gcov, gcon, beta, alpha, - sdetgam, v(b, crho, k, j, i), S, bcons, v(b, ceng, k, j, i), + prim2con::p2c(v(b, p::density(), k, j, i), vpcon, bp, + v(b, p::energy(), k, j, i), ye_prim, v(b, p::pressure(), k, j, i), + v(b, p::gamma1(), k, j, i), gcov, gcon, beta, alpha, sdetgam, + v(b, c::density(), k, j, i), S, bcons, v(b, c::energy(), k, j, i), ye_cons, sig); - v(b, cmom_lo, k, j, i) = S[0]; - v(b, cmom_lo + 1, k, j, i) = S[1]; - v(b, cmom_hi, k, j, i) = S[2]; - if (pye > 0) v(b, cye, k, j, i) = ye_cons; - for (int m = slo; m <= shi; m++) { - v(b, m, k, j, i) = sig[m - slo]; + v(b, c::momentum(0), k, j, i) = S[0]; + v(b, c::momentum(1), k, j, i) = S[1]; + v(b, c::momentum(2), k, j, i) = S[2]; + if (v.Contains(b, p::ye())) v(b, c::ye(), k, j, i) = ye_cons; + for (int m = 0; m <= 2; m++) { + v(b, impl::cell_signal_speed(m), k, j, i) = sig[m]; } - if (irfail >= 0) { + if (v.Contains(b, ir::c2pfail())) { // If rad c2p failed, we'll fix that up subsequently - if (v(b, irfail, k, j, i) == radiation::FailFlags::success) { + if (v(b, ir::c2pfail(), k, j, i) == radiation::FailFlags::success) { for (int ispec = 0; ispec < num_species; ispec++) { typename CLOSURE::LocalGeometryType g(geom, CellLocation::Cent, b, k, j, i); Vec con_v{vpcon[0] / W, vpcon[1] / W, vpcon[2] / W}; CLOSURE c(con_v, &g); - Real E = v(b, idx_E(ispec), k, j, i) / sdetgam; + Real E = v(b, cr::E(ispec), k, j, i) / sdetgam; Vec cov_F; - SPACELOOP(ii) { cov_F(ii) = v(b, idx_F(ispec, ii), k, j, i) / sdetgam; } + SPACELOOP(ii) { cov_F(ii) = v(b, cr::F(ispec, ii), k, j, i) / sdetgam; } Tens2 con_tilPi; Real J; Vec cov_H; - if (iTilPi.IsValid()) { + if (v.Contains(b, ir::tilPi())) { SPACELOOP2(ii, jj) { - con_tilPi(ii, jj) = v(b, iTilPi(ispec, ii, jj), k, j, i); + con_tilPi(ii, jj) = v(b, ir::tilPi(ispec, ii, jj), k, j, i); } } else { Real xi = 0.; @@ -341,8 +300,8 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) { c.Con2Prim(E, cov_F, con_tilPi, &J, &cov_H); - v(b, idx_J(ispec), k, j, i) = J; - SPACELOOP(ii) { v(b, idx_H(ispec, ii), k, j, i) = cov_H(ii) / J; } + v(b, pr::J(ispec), k, j, i) = J; + SPACELOOP(ii) { v(b, pr::H(ispec, ii), k, j, i) = cov_H(ii) / J; } // Floors and ceilings will be applied subsequently by ApplyFloors task } @@ -380,7 +339,5 @@ TaskStatus ConservedToPrimitiveFixup(T *rc) { return TaskStatus::fail; } -template TaskStatus -ConservedToPrimitiveFixup>(MeshBlockData *rc); - +template TaskStatus ConservedToPrimitiveFixup>(MeshData *rc); } // namespace fixup diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp index 82aa5826..b61fc69c 100644 --- a/src/fluid/fluid.cpp +++ b/src/fluid/fluid.cpp @@ -65,9 +65,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { std::string c2p_method = pin->GetOrAddString("fluid", "c2p_method", "robust"); params.Add("c2p_method", c2p_method); if (c2p_method == "robust") { - params.Add("c2p_func", ConservedToPrimitiveRobust>); + params.Add("c2p_func", ConservedToPrimitiveRobust>); } else if (c2p_method == "classic") { - params.Add("c2p_func", ConservedToPrimitiveClassic>); + params.Add("c2p_func", ConservedToPrimitiveClassic>); } else { PARTHENON_THROW("Invalid c2p_method."); } @@ -330,7 +330,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add(parthenon::hist_param_key, hst_vars); // Fill Derived and Estimate Timestep - physics->FillDerivedBlock = ConservedToPrimitive>; + physics->FillDerivedMesh = ConservedToPrimitive>; physics->EstimateTimestepBlock = EstimateTimestepBlock; return physics; @@ -451,6 +451,13 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa const IndexRange &kb) { // TODO(JMM): This one will not work with meshblock packs auto *pmb = rc->GetParentPointer(); + using parthenon::MakePackDescriptor; + Mesh *pmesh = rc->GetMeshPointer(); + auto &resolved_pkgs = pmesh->resolved_packages; + const int ndim = pmesh->ndim; + static auto desc = MakePackDescriptor<>(resolved_pkgs.get()); + auto v = desc.GetPack(rc); + const int nblocks = v.GetNBlocks(); StateDescriptor *fix_pkg = pmb->packages.Get("fixup").get(); auto bounds = fix_pkg->Param("bounds"); @@ -468,15 +475,15 @@ TaskStatus ConservedToPrimitiveRobust(T *rc, const IndexRange &ib, const IndexRa StateDescriptor *eos_pkg = pmb->packages.Get("eos").get(); auto eos = eos_pkg->Param("d.EOS"); auto geom = Geometry::GetCoordinateSystem(rc); - auto coords = pmb->coords; // TODO(JCD): move the setting of this into the solver so we can call this on MeshData auto fail = rc->Get(internal_variables::fail::name()).data; parthenon::par_for( - DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, invert.NumBlocks() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, nblocks - 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) { + auto coords = v.GetCoordinates(b); auto status = invert(geom, eos, coords, k, j, i); fail(k, j, i) = (status == con2prim_robust::ConToPrimStatus::success ? con2prim_robust::FailFlags::success @@ -491,6 +498,12 @@ TaskStatus ConservedToPrimitiveClassic(T *rc, const IndexRange &ib, const IndexR const IndexRange &kb) { using namespace con2prim; auto *pmesh = rc->GetMeshPointer(); + using parthenon::MakePackDescriptor; + auto &resolved_pkgs = pmesh->resolved_packages; + const int ndim = pmesh->ndim; + static auto desc = MakePackDescriptor<>(resolved_pkgs.get()); + auto v = desc.GetPack(rc); + const int nblocks = v.GetNBlocks(); StateDescriptor *fix_pkg = pmesh->packages.Get("fixup").get(); auto bounds = fix_pkg->Param("bounds"); @@ -510,22 +523,22 @@ TaskStatus ConservedToPrimitiveClassic(T *rc, const IndexRange &ib, const IndexR // TODO(JCD): revisit this. don't think it's required anymore. in fact the // original performance thing was related to the loop being a reduce parthenon::par_for( - DEFAULT_LOOP_PATTERN, "ConToPrim::Setup", DevExecSpace(), 0, invert.NumBlocks() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DEFAULT_LOOP_PATTERN, "ConToPrim::Setup", DevExecSpace(), 0, nblocks - 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) { invert.Setup(geom, k, j, i); }); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, invert.NumBlocks() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DEFAULT_LOOP_PATTERN, "ConToPrim::Solve", DevExecSpace(), 0, nblocks - 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) { auto status = invert(eos, k, j, i); fail(k, j, i) = (status == ConToPrimStatus::success ? FailFlags::success : FailFlags::fail); }); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "ConToPrim::Finalize", DevExecSpace(), 0, - invert.NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DEFAULT_LOOP_PATTERN, "ConToPrim::Finalize", DevExecSpace(), 0, nblocks - 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) { invert.Finalize(eos, geom, k, j, i); }); diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 38ceb8f3..7f5fccb9 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -111,19 +111,11 @@ void PhoebusDriver::PostInitializationCommunication() { const int num_partitions = pmesh->DefaultNumPartitions(); - TaskRegion &async_region_1 = tc.AddRegion(blocks.size()); - for (int ib = 0; ib < blocks.size(); ib++) { - auto pmb = blocks[ib].get(); - auto &tl = async_region_1[ib]; - auto &sc = pmb->meshblock_data.Get(); - auto apply_floors = - tl.AddTask(none, fixup::ApplyFloors>, sc.get()); - } - TaskRegion &sync_region_1 = tc.AddRegion(num_partitions); for (int ib = 0; ib < num_partitions; ib++) { auto &md = pmesh->mesh_data.GetOrAdd("base", ib); auto &tl = sync_region_1[ib]; + auto apply_floors = tl.AddTask(none, fixup::ApplyFloors>, sc.get()); const auto any = parthenon::BoundaryType::any; const auto local = parthenon::BoundaryType::local; @@ -523,6 +515,19 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { } // Fix up flux failures + TaskRegion &sync_region_new = tc.AddRegion(num_partitions); + for (int ip = 0; ip < num_partitions; ip++) { + auto &tl = sync_region_new[ip]; + auto &sc1 = pmesh->mesh_data.GetOrAdd(stage_name[stage], ip); + // fill in derived fields + auto fill_derived = + tl.AddTask(none, parthenon::Update::FillDerived>, sc1.get()); + + auto fixup = tl.AddTask(fill_derived, + fixup::ConservedToPrimitiveFixup>, sc1.get()); + auto floors = tl.AddTask(radfixup, fixup::ApplyFloors>, sc1.get()); + } + TaskRegion &async_region_2 = tc.AddRegion(num_independent_task_lists); for (int ib = 0; ib < num_independent_task_lists; ib++) { auto pmb = blocks[ib].get(); @@ -535,19 +540,9 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { auto &sc1 = pmb->meshblock_data.Get(stage_name[stage]); using MDT = std::remove_pointer::type; - // fill in derived fields - auto fill_derived = - tl.AddTask(none, parthenon::Update::FillDerived>, sc1.get()); - - auto fixup = tl.AddTask( - fill_derived, fixup::ConservedToPrimitiveFixup>, sc1.get()); - auto radfixup = tl.AddTask( fixup, fixup::RadConservedToPrimitiveFixup>, sc1.get()); - auto floors = - tl.AddTask(radfixup, fixup::ApplyFloors>, sc1.get()); - TaskID gas_rad_int(0); if (rad_mocmc_active) { auto impl_update = @@ -562,12 +557,11 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { beta * dt, fluid_active); gas_rad_int = gas_rad_int | impl_update; } - if (rad_moments_active) { // Only apply floors because MomentFluidSource already ensured that a sensible state // was returned auto floors = - tl.AddTask(gas_rad_int, fixup::ApplyFloors>, sc1.get()); + tl.AddTask(gas_rad_int, fixup::ApplyFloors>, sc1.get()); } } diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp index d510af9e..76ad8d4e 100644 --- a/src/phoebus_utils/variables.hpp +++ b/src/phoebus_utils/variables.hpp @@ -17,6 +17,8 @@ #include #include +using parthenon::variable_names::ANYDIM; + #define VARIABLE(ns, varname) \ struct varname : public parthenon::variable_names::base_t { \ template \ @@ -41,6 +43,13 @@ static std::string name() { return #varstring; } \ } +#define TENSOR_VARIABLE(ns, varname, ...) \ + struct varname : public parthenon::variable_names::base_t { \ + template \ + KOKKOS_INLINE_FUNCTION varname(Ts &&...args) : parthenon::variable_names::base_t(std::forward(args)...) {} \ + static std::string name() { return #ns "." #varname; } \ + } + namespace fluid_prim { VARIABLE(p, density); VARIABLE(p, velocity); @@ -62,12 +71,12 @@ VARIABLE(c, ye); namespace radmoment_prim { VARIABLE(r.p, J); -VARIABLE(r.p, H); + TENSOR_VARIABLE(r.p, H, ANYDIM, 3); } // namespace radmoment_prim namespace radmoment_cons { VARIABLE(r.c, E); -VARIABLE(r.c, F); + TENSOR_VARIABLE(r.c, F, ANYDIM, 3); } // namespace radmoment_cons namespace radmoment_internal { @@ -81,7 +90,7 @@ VARIABLE(r.i, dJ); VARIABLE(r.i, kappaJ); VARIABLE(r.i, kappaH); VARIABLE(r.i, JBB); -VARIABLE(r.i, tilPi); +TENSOR_VARIABLE(r.i, tilPi, ANYDIM, 3); VARIABLE(r.i, kappaH_mean); VARIABLE(r.i, c2pfail); VARIABLE(r.i, srcfail);