Address XDMF/Visit Issues

CHANGELOG.md

@@ -3,7 +3,8 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
-- [[PR1020]](https://github.com/parthenon-hpc-lab/parthenon/pull/1020) Add bi- and trilinear interpolation routines
+- [[PR 1079]](https://github.com/parthenon-hpc-lab/parthenon/pull/1079) Address XDMF/Visit Issues
+- [[PR 1020]](https://github.com/parthenon-hpc-lab/parthenon/pull/1020) Add bi- and trilinear interpolation routines
 - [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
 - [[PR 1037]](https://github.com/parthenon-hpc-lab/parthenon/pull/1037) Add SwarmPacks
 - [[PR 1068]](https://github.com/parthenon-hpc-lab/parthenon/pull/1068) Add ability to dump sparse pack contents as a string
@@ -19,6 +20,7 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1079]](https://github.com/parthenon-hpc-lab/parthenon/pull/1079) Address XDMF/Visit Issues
 - [[PR 1071]](https://github.com/parthenon-hpc-lab/parthenon/pull/1070) Fix bug in static mesh refinement related to redefinition of Mesh::root_level
 - [[PR 1073]](https://github.com/parthenon-hpc-lab/parthenon/pull/1073) Fix bug in AMR and sparse restarts
 - [[PR 1070]](https://github.com/parthenon-hpc-lab/parthenon/pull/1070) Correctly exclude flux vars from searches by default

src/outputs/output_parameters.hpp

@@ -50,12 +50,13 @@ struct OutputParameters {
   bool sparse_seed_nans;
   int hdf5_compression_level;
   bool write_xdmf;
+  bool write_swarm_xdmf;
   // TODO(felker): some of the parameters in this class are not initialized in constructor
   OutputParameters()
       : block_number(0), next_time(0.0), dt(-1.0), file_number(0),
         include_ghost_zones(false), cartesian_vector(false),
         single_precision_output(false), sparse_seed_nans(false),
-        hdf5_compression_level(5), write_xdmf(false) {}
+        hdf5_compression_level(5), write_xdmf(false), write_swarm_xdmf(false) {}
 };
 
 } // namespace parthenon

src/outputs/outputs.cpp

@@ -280,6 +280,9 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) {
         }
 #ifdef ENABLE_HDF5
         op.write_xdmf = pin->GetOrAddBoolean(op.block_name, "write_xdmf", true);
+        op.write_swarm_xdmf =
+            (restart) ? false
+                      : pin->GetOrAddBoolean(op.block_name, "write_swarm_xdmf", false);
         pnew_type = new PHDF5Output(op, restart);
 #else
         msg << "### FATAL ERROR in Outputs constructor" << std::endl

src/outputs/parthenon_hdf5.cpp

@@ -1,6 +1,6 @@
 //========================================================================================
 // Parthenon performance portable AMR framework
-// Copyright(C) 2020-2023 The Parthenon collaboration
+// Copyright(C) 2020-2024 The Parthenon collaboration
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
 // (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
@@ -34,6 +34,7 @@
 
 #include "driver/driver.hpp"
 #include "interface/metadata.hpp"
+#include "interface/swarm_default_names.hpp"
 #include "mesh/mesh.hpp"
 #include "mesh/meshblock.hpp"
 #include "outputs/output_utils.hpp"
@@ -485,6 +486,17 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
       local_count[rank] = swinfo.count_on_rank;
       global_count[rank] = swinfo.global_count;
     };
+    auto SetCountsParticlePositions = [&](const SwarmInfo &swinfo) {
+      for (int i = 0; i < 6; ++i) {
+        local_offset[i] = 0; // reset everything
+        local_count[i] = 0;
+        global_count[i] = 0;
+      }
+      local_offset[0] = swinfo.global_offset;
+      local_count[0] = swinfo.count_on_rank;
+      global_count[0] = swinfo.global_count;
+      local_count[1] = global_count[1] = 3;
+    };
     auto &int_vars = std::get<SwarmInfo::MapToVarVec<int>>(swinfo.vars);
     for (auto &[vname, swmvarvec] : int_vars) {
       const auto &vinfo = swinfo.var_info.at(vname);
@@ -493,14 +505,40 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
       HDF5WriteND(g_var, vname, host_data.data(), vinfo.tensor_rank + 1, local_offset,
                   local_count, global_count, pl_xfer, H5P_DEFAULT);
     }
+    std::vector<Real> pos_tmp; // tmp vector to (potentially) hold particle positions
     auto &rvars = std::get<SwarmInfo::MapToVarVec<Real>>(swinfo.vars);
     for (auto &[vname, swmvarvec] : rvars) {
       const auto &vinfo = swinfo.var_info.at(vname);
       auto host_data = swinfo.FillHostBuffer(vname, swmvarvec);
       SetCounts(swinfo, vinfo);
       HDF5WriteND(g_var, vname, host_data.data(), vinfo.tensor_rank + 1, local_offset,
                   local_count, global_count, pl_xfer, H5P_DEFAULT);
+      if (output_params.write_swarm_xdmf &&
+          (vname == swarm_position::x::name() || vname == swarm_position::y::name() ||
+           vname == swarm_position::z::name())) {
+        pos_tmp.insert(pos_tmp.end(), host_data.begin(), host_data.end());
+      }
     }
+    if (output_params.write_swarm_xdmf) {
+      // TODO(@pdmullen): Here and above, we have worked with temp vectors pos_tmp and
+      // swarm_positions so that we can take the existing swarm position data structures
+      // and recast them into a format that XDMF/VisIt prefers.  Future efforts may (1)
+      // eliminate the extra geometry dump and/or (2) directly write the auxillary
+      // positions via low-level HDF writes, such that the vector manipulation below is
+      // unnecessary.
+      const int npart = pos_tmp.size() / 3;
+      std::vector<Real> swarm_positions(pos_tmp.size());
+      int spcnt = 0;
+      for (int i = 0; i < npart; ++i) {
+        swarm_positions[spcnt++] = pos_tmp[i];
+        swarm_positions[spcnt++] = pos_tmp[npart + i];
+        swarm_positions[spcnt++] = pos_tmp[2 * npart + i];
+      }
+      SetCountsParticlePositions(swinfo);
+      HDF5WriteND(g_var, "swarm_positions", swarm_positions.data(), 2, local_offset,
+                  local_count, global_count, pl_xfer, H5P_DEFAULT);
+    }
+
     // If swarm does not contain an "id" object, generate a sequential
     // one for vis.
     if (swinfo.var_info.count("id") == 0) {
@@ -515,10 +553,11 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
   }
   Kokkos::Profiling::popRegion(); // write particle data
 
-  if (output_params.write_xdmf) {
+  if (output_params.write_xdmf || output_params.write_swarm_xdmf) {
     Kokkos::Profiling::pushRegion("genXDMF");
     // generate XDMF companion file
-    XDMF::genXDMF(filename, pm, tm, theDomain, nx1, nx2, nx3, all_vars_info, swarm_info);
+    XDMF::genXDMF(filename, pm, tm, theDomain, nx1, nx2, nx3, all_vars_info, swarm_info,
+                  output_params.write_xdmf, output_params.write_swarm_xdmf);
     Kokkos::Profiling::popRegion(); // genXDMF
   }