WIP: Update submodules, API changes, tracer swarm packs, formatting CI

.github/workflows/formatting.yml

@@ -18,6 +18,11 @@ jobs:
           run: export DEBIAN_FRONTEND=noninteractive
         - name: install dependencies
           run: |
-            sudo apt-get install -y --force-yes -qq git clang-format-12
+            sudo apt-get install -y --force-yes -qq git
+        - name: Install pip, clang-format
+          uses: BSFishy/pip-action@v1
+          with:
+            packages: |
+               clang-format==12.0.1
         - name: check formatting
-          run: find . -regex '.*\.\(cpp\|hpp\)' | xargs clang-format-12 -style=file -i && git diff --exit-code --ignore-submodules
+          run: find . -regex '.*\.\(cpp\|hpp\)' | xargs clang-format -style=file -i && git diff --exit-code --ignore-submodules

CMakeLists.txt

@@ -166,6 +166,7 @@ set(SINGULARITY_USE_FORTRAN OFF CACHE BOOL "" FORCE)
 set(SINGULARITY_FORCE_SUBMODULE_MODE ON CACHE BOOL "" FORCE)
 set(SINGULARITY_USE_KOKKOSKERNELS OFF CACHE BOOL "" FORCE)
 set(SINGULARITY_BUILD_CLOSURE OFF CACHE BOOL "" FORCE)
+set(SINGULARITY_USE_TRUE_LOG_GRIDDING ON CACHE BOOL "" FORCE)
 
 # Kill cmake's package registry because it can interfere
 set(CMAKE_FIND_PACKAGE_NO_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
@@ -184,6 +185,9 @@ include(cmake/Geometry.cmake)
 # Fluid
 include(cmake/Fluid.cmake)
 
+# Transport
+include(cmake/Transport.cmake)
+
 # Phoebus src
 message("\nConfiguring src")
 add_subdirectory(src)

cmake/Transport.cmake

@@ -0,0 +1,13 @@
+# © 2021-2025. Triad National Security, LLC. All rights reserved.  This
+# program was produced under U.S. Government contract 89233218CNA000001
+# for Los Alamos National Laboratory (LANL), which is operated by Triad
+# National Security, LLC for the U.S.  Department of Energy/National
+# Nuclear Security Administration. All rights in the program are
+# reserved by Triad National Security, LLC, and the U.S. Department of
+# Energy/National Nuclear Security Administration. The Government is
+# granted for itself and others acting on its behalf a nonexclusive,
+# paid-up, irrevocable worldwide license in this material to reproduce,
+# prepare derivative works, distribute copies to the public, perform
+# publicly and display publicly, and to permit others to do so.
+
+set(MOCMC_NUM_SPECIES 3 CACHE INTEGER "Number of radiation species for MOCMC.")

src/compile_constants.hpp.in

@@ -16,6 +16,8 @@
 
 #define NCONS_MAX (@MAX_NUMBER_CONSERVED_VARS@)
 
+#define MOCMC_NUM_SPECIES (@MOCMC_NUM_SPECIES@)
+
 #define PHOEBUS_GEOMETRY Geometry::@PHOEBUS_GEOMETRY@
 #define GEOMETRY_MESH @GEOMETRY_MESH@
 #define GEOMETRY_MESH_BLOCK @GEOMETRY_MESH_BLOCK@

src/fixup/fixup.cpp

@@ -15,7 +15,6 @@
 
 #include "fixup.hpp"
 
-#include <bvals/bvals_interfaces.hpp>
 #include <defs.hpp>
 
 #include "analysis/history.hpp"
@@ -731,6 +730,8 @@ TaskStatus FixFluxes(MeshBlockData<Real> *rc) {
   const std::string ox1_bc = fluid->Param<std::string>("ox1_bc");
   const std::string ix2_bc = fluid->Param<std::string>("ix2_bc");
   const std::string ox2_bc = fluid->Param<std::string>("ox2_bc");
+  const std::string ix3_bc = fluid->Param<std::string>("ix3_bc");
+  const std::string ox3_bc = fluid->Param<std::string>("ox3_bc");
 
   auto rad = pmb->packages.Get("radiation");
   int num_species = 0;
@@ -993,7 +994,7 @@ TaskStatus FixFluxes(MeshBlockData<Real> *rc) {
   if (ndim == 2) return TaskStatus::complete;
 
   // x3-direction
-  if (pmb->boundary_flag[BoundaryFace::inner_x3] == BoundaryFlag::outflow) {
+  if (ix3_bc == "outflow") {
     PackIndexMap imap;
     auto v =
         rc->PackVariablesAndFluxes(std::vector<std::string>({c::density::name()}),
@@ -1004,7 +1005,7 @@ TaskStatus FixFluxes(MeshBlockData<Real> *rc) {
         ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) {
           v.flux(X3DIR, crho, k, j, i) = std::min(v.flux(X3DIR, crho, k, j, i), 0.0);
         });
-  } else if (pmb->boundary_flag[BoundaryFace::inner_x3] == BoundaryFlag::reflect) {
+  } else if (ix3_bc == "reflect") {
     PackIndexMap imap;
     auto v = rc->PackVariablesAndFluxes(
         std::vector<std::string>({c::density::name(), c::energy::name(), cr::E::name()}),
@@ -1025,7 +1026,7 @@ TaskStatus FixFluxes(MeshBlockData<Real> *rc) {
           }
         });
   }
-  if (pmb->boundary_flag[BoundaryFace::outer_x3] == BoundaryFlag::outflow) {
+  if (ox3_bc == "outflow") {
     PackIndexMap imap;
     auto v =
         rc->PackVariablesAndFluxes(std::vector<std::string>({c::density::name()}),
@@ -1036,7 +1037,7 @@ TaskStatus FixFluxes(MeshBlockData<Real> *rc) {
         jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) {
           v.flux(X3DIR, crho, k, j, i) = std::max(v.flux(X3DIR, crho, k, j, i), 0.0);
         });
-  } else if (pmb->boundary_flag[BoundaryFace::outer_x3] == BoundaryFlag::reflect) {
+  } else if (ox3_bc == "reflect") {
     PackIndexMap imap;
     auto v = rc->PackVariablesAndFluxes(
         std::vector<std::string>({c::density::name(), c::energy::name(), cr::E::name()}),

src/fixup/fixup_c2p.cpp

@@ -15,7 +15,6 @@
 
 #include "fixup.hpp"
 
-#include <bvals/bvals_interfaces.hpp>
 #include <defs.hpp>
 
 #include "fluid/con2prim_robust.hpp"

src/fixup/fixup_netfield.cpp

@@ -15,7 +15,6 @@
 
 #include "fixup.hpp"
 
-#include <bvals/bvals_interfaces.hpp>
 #include <defs.hpp>
 #include <globals.hpp>
 

src/fixup/fixup_particles.cpp

@@ -38,11 +38,12 @@ TaskStatus PurgeParticles(MeshBlockData<Real> *rc, const std::string swarmName)
   /* only do this when FMKS is used. */
   if constexpr (std::is_same<PHOEBUS_GEOMETRY, Geometry::FMKS>::value) {
     auto *pmb = rc->GetParentPointer();
-    auto &swarm = pmb->swarm_data.Get()->Get(swarmName);
+    // auto &swarm = rc->swarm_data.Get()->Get(swarmName);
+    auto &swarm = rc->GetSwarmData()->Get(swarmName);
 
-    auto &x = swarm->Get<Real>("x").Get();
-    auto &y = swarm->Get<Real>("y").Get();
-    auto &z = swarm->Get<Real>("z").Get();
+    auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+    auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+    auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
 
     auto swarm_d = swarm->GetDeviceContext();
 

src/fixup/fixup_radc2p.cpp

@@ -15,7 +15,6 @@
 
 #include "fixup.hpp"
 
-#include <bvals/bvals_interfaces.hpp>
 #include <defs.hpp>
 
 #include "fluid/con2prim_robust.hpp"

src/fixup/fixup_src.cpp

@@ -15,7 +15,6 @@
 
 #include "fixup.hpp"
 
-#include <bvals/bvals_interfaces.hpp>
 #include <defs.hpp>
 
 #include "fluid/con2prim_robust.hpp"

src/fluid/fluid.cpp

@@ -183,6 +183,10 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   params.Add("ix2_bc", ix2_bc);
   const std::string ox2_bc = pin->GetString("phoebus", "ox2_bc");
   params.Add("ox2_bc", ox2_bc);
+  const std::string ix3_bc = pin->GetString("phoebus", "ix3_bc");
+  params.Add("ix3_bc", ix3_bc);
+  const std::string ox3_bc = pin->GetString("phoebus", "ox3_bc");
+  params.Add("ox3_bc", ox3_bc);
 
   int ndim = 1;
   if (pin->GetInteger("parthenon/mesh", "nx3") > 1)
@@ -768,6 +772,7 @@ TaskStatus CalculateFluxes(MeshBlockData<Real> *rc) {
 }
 
 TaskStatus FluxCT(MeshBlockData<Real> *rc) {
+  using parthenon::MakePackDescriptor;
   Mesh *pmesh = rc->GetMeshPointer();
   auto &fluid = pmesh->packages.Get("fluid");
   if (!fluid->Param<bool>("mhd") || !fluid->Param<bool>("active") ||
@@ -780,49 +785,54 @@ TaskStatus FluxCT(MeshBlockData<Real> *rc) {
   IndexRange jb = rc->GetBoundsJ(IndexDomain::interior);
   IndexRange kb = rc->GetBoundsK(IndexDomain::interior);
 
-  auto f1 = rc->Get(fluid_cons::bfield::name()).flux[X1DIR];
-  auto f2 = rc->Get(fluid_cons::bfield::name()).flux[X2DIR];
-  auto f3 = rc->Get(fluid_cons::bfield::name()).flux[X3DIR];
+  auto &resolved_pkgs = pmesh->resolved_packages;
+  static auto desc = MakePackDescriptor<fluid_cons::bfield>(
+      resolved_pkgs.get(), {}, {parthenon::PDOpt::WithFluxes});
+  auto v = desc.GetPack(rc);
   auto emf = rc->Get(internal_variables::emf::name()).data;
 
   if (ndim == 2) {
     parthenon::par_for(
         DEFAULT_LOOP_PATTERN, "FluxCT::EMF::2D", DevExecSpace(), kb.s, kb.e, jb.s,
         jb.e + 1, ib.s, ib.e + 1, KOKKOS_LAMBDA(const int k, const int j, const int i) {
-          emf(k, j, i) = 0.25 * (f1(1, k, j, i) + f1(1, k, j - 1, i) - f2(0, k, j, i) -
-                                 f2(0, k, j, i - 1));
+          emf(k, j, i) =
+              0.25 * (v.flux(0, X1DIR, 1, k, j, i) + v.flux(0, X1DIR, 1, k, j - 1, i) -
+                      v.flux(0, X2DIR, 0, k, j, i) - v.flux(0, X2DIR, 0, k, j, i - 1));
         });
     parthenon::par_for(
         DEFAULT_LOOP_PATTERN, "FluxCT::Flux::2D", DevExecSpace(), kb.s, kb.e, jb.s,
         jb.e + 1, ib.s, ib.e + 1, KOKKOS_LAMBDA(const int k, const int j, const int i) {
-          f1(0, k, j, i) = 0.0;
-          f1(1, k, j, i) = 0.5 * (emf(k, j, i) + emf(k, j + 1, i));
-          f2(0, k, j, i) = -0.5 * (emf(k, j, i) + emf(k, j, i + 1));
-          f2(1, k, j, i) = 0.0;
+          v.flux(0, X1DIR, 0, k, j, i) = 0.0;
+          v.flux(0, X1DIR, 1, k, j, i) = 0.5 * (emf(k, j, i) + emf(k, j + 1, i));
+          v.flux(0, X2DIR, 0, k, j, i) = -0.5 * (emf(k, j, i) + emf(k, j, i + 1));
+          v.flux(0, X2DIR, 1, k, j, i) = 0.0;
         });
   } else if (ndim == 3) {
     parthenon::par_for(
         DEFAULT_LOOP_PATTERN, "FluxCT::EMF::3D", DevExecSpace(), kb.s, kb.e + 1, jb.s,
         jb.e + 1, ib.s, ib.e + 1, KOKKOS_LAMBDA(const int k, const int j, const int i) {
-          emf(0, k, j, i) = 0.25 * (f2(2, k, j, i) + f2(2, k - 1, j, i) - f3(1, k, j, i) -
-                                    f3(1, k, j - 1, i));
-          emf(1, k, j, i) = -0.25 * (f1(2, k, j, i) + f1(2, k - 1, j, i) -
-                                     f3(0, k, j, i) - f3(0, k, j, i - 1));
-          emf(2, k, j, i) = 0.25 * (f1(1, k, j, i) + f1(1, k, j - 1, i) - f2(0, k, j, i) -
-                                    f2(0, k, j, i - 1));
+          emf(0, k, j, i) =
+              0.25 * (v.flux(0, X2DIR, 2, k, j, i) + v.flux(0, X2DIR, 2, k - 1, j, i) -
+                      v.flux(0, X3DIR, 1, k, j, i) - v.flux(0, X3DIR, 1, k, j - 1, i));
+          emf(1, k, j, i) =
+              -0.25 * (v.flux(0, X1DIR, 2, k, j, i) + v.flux(0, X1DIR, 2, k - 1, j, i) -
+                       v.flux(0, X3DIR, 0, k, j, i) - v.flux(0, X3DIR, 0, k, j, i - 1));
+          emf(2, k, j, i) =
+              0.25 * (v.flux(0, X1DIR, 1, k, j, i) + v.flux(0, X1DIR, 1, k, j - 1, i) -
+                      v.flux(0, X2DIR, 0, k, j, i) - v.flux(0, X2DIR, 0, k, j, i - 1));
         });
     parthenon::par_for(
         DEFAULT_LOOP_PATTERN, "FluxCT::Flux::3D", DevExecSpace(), kb.s, kb.e + 1, jb.s,
         jb.e + 1, ib.s, ib.e + 1, KOKKOS_LAMBDA(const int k, const int j, const int i) {
-          f1(0, k, j, i) = 0.0;
-          f1(1, k, j, i) = 0.5 * (emf(2, k, j, i) + emf(2, k, j + 1, i));
-          f1(2, k, j, i) = -0.5 * (emf(1, k, j, i) + emf(1, k + 1, j, i));
-          f2(0, k, j, i) = -0.5 * (emf(2, k, j, i) + emf(2, k, j, i + 1));
-          f2(1, k, j, i) = 0.0;
-          f2(2, k, j, i) = 0.5 * (emf(0, k, j, i) + emf(0, k + 1, j, i));
-          f3(0, k, j, i) = 0.5 * (emf(1, k, j, i) + emf(1, k, j, i + 1));
-          f3(1, k, j, i) = -0.5 * (emf(0, k, j, i) + emf(0, k, j + 1, i));
-          f3(2, k, j, i) = 0.0;
+          v.flux(0, X1DIR, 0, k, j, i) = 0.0;
+          v.flux(0, X1DIR, 1, k, j, i) = 0.5 * (emf(2, k, j, i) + emf(2, k, j + 1, i));
+          v.flux(0, X1DIR, 2, k, j, i) = -0.5 * (emf(1, k, j, i) + emf(1, k + 1, j, i));
+          v.flux(0, X2DIR, 0, k, j, i) = -0.5 * (emf(2, k, j, i) + emf(2, k, j, i + 1));
+          v.flux(0, X2DIR, 1, k, j, i) = 0.0;
+          v.flux(0, X2DIR, 2, k, j, i) = 0.5 * (emf(0, k, j, i) + emf(0, k + 1, j, i));
+          v.flux(0, X3DIR, 0, k, j, i) = 0.5 * (emf(1, k, j, i) + emf(1, k, j, i + 1));
+          v.flux(0, X3DIR, 1, k, j, i) = -0.5 * (emf(0, k, j, i) + emf(0, k, j + 1, i));
+          v.flux(0, X3DIR, 2, k, j, i) = 0.0;
         });
   }
 

src/fluid/riemann.hpp

@@ -70,8 +70,16 @@ class FluxState {
   }
   static void FluxVars(const std::string &var) { flux_vars.push_back(var); }
 
-  static std::vector<std::string> ReconVars() { return recon_vars; }
-  static std::vector<std::string> FluxVars() { return flux_vars; }
+  static std::vector<std::string> ReconVars() {
+    for (auto &v : recon_vars) {
+    }
+    return recon_vars;
+  }
+  static std::vector<std::string> FluxVars() {
+    for (auto &v : flux_vars) {
+    }
+    return flux_vars;
+  }
 
   KOKKOS_FORCEINLINE_FUNCTION
   int NumConserved() const { return ncons; }
@@ -229,12 +237,12 @@ class FluxState {
         pye(imap[fluid_prim::ye::name()].second),
         prs(imap[fluid_prim::pressure::name()].first),
         gm1(imap[fluid_prim::gamma1::name()].first),
-        crho(imap[fluid_cons::density::name()].first),
-        cmom_lo(imap[fluid_cons::momentum::name()].first),
-        ceng(imap[fluid_cons::energy::name()].first),
-        cb_lo(imap[fluid_cons::bfield::name()].first),
-        cb_hi(imap[fluid_cons::bfield::name()].second),
-        cye(imap[fluid_cons::ye::name()].first),
+        crho(imap[fluid_prim::density::name()].first),
+        cmom_lo(imap[fluid_prim::velocity::name()].first),
+        ceng(imap[fluid_prim::energy::name()].first),
+        cb_lo(imap[fluid_prim::bfield::name()].first),
+        cb_hi(imap[fluid_prim::bfield::name()].second),
+        cye(imap[fluid_prim::ye::name()].first),
         ncons(5 + (pb_hi - pb_lo + 1) + (cye > 0)) {
     PARTHENON_REQUIRE_THROWS(
         ncons <= NCONS_MAX,

src/geometry/geometry.cpp

@@ -54,8 +54,10 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   Metadata gcoord_cell =
       Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, cell_shape);
   // TODO(JMM): Make this actual node-centered data when available
+  // TODO: Warning: g.n.coord output broken in Parthenon currently
   std::vector<int> node_shape = {dims.nx1 + 1, dims.nx2 + 1, dims.nx3 + 1, 4};
-  Metadata gcoord_node = Metadata({Metadata::Derived, Metadata::OneCopy}, node_shape);
+  Metadata gcoord_node =
+      Metadata({Metadata::Derived, Metadata::OneCopy, Metadata::None}, node_shape);
   geometry->AddField(geometric_variables::cell_coords::name(), gcoord_cell);
   geometry->AddField(geometric_variables::node_coords::name(), gcoord_node);
 

src/microphysics/opac_phoebus/opac_phoebus.cpp

@@ -144,8 +144,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   }
 
   {
-    auto opacity_host =
-        params.Get<singularity::neutrinos::Opacity>("h.opacity_baseunits");
+    auto opacity_host = params.Get<singularity::neutrinos::Opacity>("h.opacity");
     const Real YeMin = pin->GetOrAddReal("mean_opacity", "yemin", 0.1);
     const Real YeMax = pin->GetOrAddReal("mean_opacity", "yemax", 0.5);
     const int NYe = pin->GetOrAddInteger("mean_opacity", "nye", 10);
@@ -159,10 +158,10 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       const Real lNuMin = std::log10(pin->GetOrAddReal("mean_opacity", "numin", 0.1));
       const Real lNuMax = std::log10(pin->GetOrAddReal("mean_opacity", "numax", 10.));
       const int NNu = pin->GetOrAddInteger("mean_opacity", "nnu", 100);
-      MeanOpacity mean_opac_host =
-          MeanOpacityScaleFree(opacity_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT,
-                               YeMin, YeMax, NYe, lNuMin, lNuMax, NNu);
-      auto mean_opac_device = mean_opac_host.GetOnDevice();
+      auto mean_opac_host =
+          MeanOpacityBase(opacity_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin,
+                          YeMax, NYe, lNuMin, lNuMax, NNu);
+      MeanOpacity mean_opac_device = mean_opac_host.GetOnDevice();
       params.Add("h.mean_opacity", mean_opac_host);
       params.Add("d.mean_opacity", mean_opac_device);
     } else {
@@ -176,13 +175,13 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       const Real lNuMin = std::log10(pin->GetOrAddReal("mean_opacity", "numin", 1.e10));
       const Real lNuMax = std::log10(pin->GetOrAddReal("mean_opacity", "numax", 1.e24));
       const int NNu = pin->GetOrAddInteger("mean_opacity", "nnu", 100);
-      MeanOpacityCGS cgs_mean_opacity =
-          MeanOpacityCGS(opacity_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin,
-                         YeMax, NYe, lNuMin, lNuMax, NNu);
-      MeanOpacity mean_opac_host =
-          MeanNonCGSUnits<MeanOpacityCGS>(std::forward<MeanOpacityCGS>(cgs_mean_opacity),
-                                          time_unit, mass_unit, length_unit, temp_unit);
-      auto mean_opac_device = mean_opac_host.GetOnDevice();
+      auto cgs_mean_opacity =
+          MeanOpacityBase(opacity_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin,
+                          YeMax, NYe, lNuMin, lNuMax, NNu);
+      auto mean_opac_host = MeanNonCGSUnits<MeanOpacityBase>(
+          std::forward<MeanOpacityBase>(cgs_mean_opacity), time_unit, mass_unit,
+          length_unit, temp_unit);
+      MeanOpacity mean_opac_device = mean_opac_host.GetOnDevice();
       params.Add("h.mean_opacity", mean_opac_host);
       params.Add("d.mean_opacity", mean_opac_device);
     }