custom coords

CHANGELOG.md

@@ -3,7 +3,8 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
-
+- [[PR1039]](https://github.com/parthenon-hpc-lab/parthenon/pull/1039) Add ability to output custom coordinate positions for Visit/Paraview
+- [[PR1019](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Enable output for non-cell-centered variables
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR 1024]](https://github.com/parthenon-hpc-lab/parthenon/pull/1024) Add .outN. to history output filenames
@@ -21,6 +22,7 @@
 
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR1019](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Remove support for file formats < 3
 
 
 ## Release 24.03
@@ -39,7 +41,6 @@ Date: 2024-03-21
 - [[PR 996]](https://github.com/parthenon-hpc-lab/parthenon/pull/996) Remove dynamic allocations from swarm particle creation
 
 ### Changed (changing behavior/API/variables/...)
-- [[PR1019](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Enable output for non-cell-centered variables
 - [[PR 973]](https://github.com/parthenon-hpc-lab/parthenon/pull/973) Multigrid performance upgrades
 
 ### Fixed (not changing behavior/API/variables/...)
@@ -58,7 +59,6 @@ Date: 2024-03-21
 - [[PR 982]](https://github.com/parthenon-hpc-lab/parthenon/pull/982) add some gut check testing for parthenon-VIBE
 
 ### Incompatibilities (i.e. breaking changes)
-- [[PR1019](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Remove support for file formats < 3
 - [[PR 987]](https://github.com/parthenon-hpc-lab/parthenon/pull/987) Change the API for what was IterativeTasks
 - [[PR 974]](https://github.com/parthenon-hpc-lab/parthenon/pull/974) Change GetParentPointer to always return T*
 - [[PR 996]](https://github.com/parthenon-hpc-lab/parthenon/pull/996) Remove dynamic allocations from swarm particle creation

CMakeLists.txt

@@ -62,9 +62,9 @@ include(cmake/Format.cmake)
 include(cmake/Lint.cmake)
 
 # regression test reference data
-set(REGRESSION_GOLD_STANDARD_VER 21 CACHE STRING "Version of gold standard to download and use")
+set(REGRESSION_GOLD_STANDARD_VER 22 CACHE STRING "Version of gold standard to download and use")
 set(REGRESSION_GOLD_STANDARD_HASH
-  "SHA512=e16b14272915b4607965e5900961402f6da96dc13da8ea3c3d213d61f82d3a1dded08c40a9ab644aa3409d93a045bba360a90a43dc289b24f525878f9ba50890"
+  "SHA512=c64f34b6841569c74c3918aa33d5c1fe6795b6ba124e25be7184c48f589ebe634062955e97969e68d252e59d8d0fce5d65e32b4252e2c2537f9f49ebfdf37ee0"
   CACHE STRING "Hash of default gold standard file to download")
 option(REGRESSION_GOLD_STANDARD_SYNC "Automatically sync gold standard files." ON)
 

doc/sphinx/src/interface/metadata.rst

@@ -79,6 +79,9 @@ These flags specify how a variable interacts with I/O. Enable them to
 enable output properties.
 
 -  ``Metadata::Restart`` implies a variable is required in restart files
+- ``Metadata::CoordinatesVec`` implies a variable should be used to
+   describe coordinate positions of nodes. Must be a node-centered
+   variable with 3 components.
 
 Tensor properties and boundaries
 --------------------------------

doc/sphinx/src/outputs.rst

@@ -408,6 +408,65 @@ cases, the ``.xdmf`` files should be opened. In ParaView, select the
    for ParaView and VisIt. However, our python tooling does support
    all mesh locations.
 
+Tying non-standard coordinates to visualization tools
+------------------------------------------------------
+
+By default, Parthenon outputs the positions of faces on each block in
+``X1``, ``X2`` and ``X3`` and assumes these correspond to the ``x``,
+``y``, and ``z`` components of the nodes on the mesh. However, some
+applications may apply cooridnate transformations or use moving
+meshes. In these cases, the above strategy will not provide intuitive
+plots.
+
+For these applications we provide a special ``Metadata`` flag. If you
+mark a node-centered 3-vector variable with the the flag
+``Metadata::CoordinatesVec``, and fill it with the ``x``, ``y``, and
+``z`` values of your node positions, Parthenon will specify these
+values should be used by visualization software such as Visit or
+Paraview.
+
+For example, in your package ``Initialize`` function, you might
+declare something like:
+
+.. code:: cpp
+
+   pkg->AddField("locations",
+     Metadata({Metadata::Node, Metadata::CoordinatesVec, Metadata::Derived, Metadata::OneCopy},
+     std::vector<int>{3}));
+
+and then (trivially) if you set
+
+.. code:: cpp
+
+   pman.app_input->InitMeshBlockUserData = SetGeometryBlock;
+
+for
+
+.. code:: cpp
+
+   void SetGeometryBlock(MeshBlock *pmb, ParameterInput *pin) {
+     /* boiler plate to build a pack object */
+     parthenon::par_for(DEFAULT_LOOP_PATTERN, "positions", DevExecSpace(), 0, pack.GetNBlocks() - 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) {
+          const auto &coords = pack.GetCoordinates(b);
+          pack(b, 0, k, j, i) = coords.X<X1DIR, parthenon::TopologicalElement::NN>(k, j, i);
+          pack(b, 1, k, j, i) = coords.X<X2DIR, parthenon::TopologicalElement::NN>(k, j, i);
+          pack(b, 2, k, j, i) = coords.X<X3DIR, parthenon::TopologicalElement::NN>(k, j, i);
+     });
+     return;
+   }
+
+then the code will set the nodal values to their trivial coordinate
+values and these will be used for visualization. In a more non-trivial
+example, ``SetGeometryBlock`` might apply a coordinate
+transformation. Or actually evolve ``"locations"``.
+
+.. warning::
+
+   Non-standard coordinates are not supported in XDMF for 1D meshes
+   and Parthenon will revert to the traditional output in 1D.
+
 Preparing outputs for ``yt``
 ----------------------------
 

example/sparse_advection/sparse_advection_package.cpp

@@ -72,8 +72,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   // add sparse field
   {
     Metadata m({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes,
-                Metadata::FillGhost, Metadata::Sparse},
-               std::vector<int>({1}));
+                Metadata::FillGhost, Metadata::Sparse});
     SparsePool pool("sparse", m);
 
     for (int sid = 0; sid < NUM_FIELDS; ++sid) {

src/interface/metadata.cpp

@@ -184,6 +184,88 @@ std::ostream &operator<<(std::ostream &os, const parthenon::Metadata &m) {
   return os;
 }
 
+// Return true if the flags constraints are satisfied, false otherwise. If throw_on_fail
+// is true, throw a descriptive exception when invalid
+bool Metadata::IsValid(bool throw_on_fail) const {
+  bool valid = true;
+
+  // Topology
+  if (CountSet({None, Node, Edge, Face, Cell}) != 1) {
+    valid = false;
+    if (throw_on_fail) {
+      PARTHENON_THROW("Exactly one topology flag must be set");
+    }
+  }
+
+  // Role
+  if (CountSet({Private, Provides, Requires, Overridable}) != 1) {
+    valid = false;
+    if (throw_on_fail) {
+      PARTHENON_THROW("Exactly one role flag must be set");
+    }
+  }
+
+  // Shape
+  if (CountSet({Vector, Tensor}) > 1) {
+    valid = false;
+    if (throw_on_fail) {
+      PARTHENON_THROW("At most one shape flag can be set");
+    }
+  }
+
+  // Coordinates
+  if (IsSet(CoordinatesVec)) {
+    if (Where() != Node) {
+      valid = false;
+      if (throw_on_fail) {
+        PARTHENON_THROW("Coordinate field must be node-centered");
+      }
+    }
+    if (shape_.size() != 1) {
+      valid = false;
+      if (throw_on_fail) {
+        PARTHENON_THROW("Coordinate field must be tensor rank 1");
+      }
+    }
+    if (shape_[0] != 3) {
+      valid = false;
+      if (throw_on_fail) {
+        PARTHENON_THROW(
+            "Coordinate field must be 3-vector. (Does not need Vector metadata flag).");
+      }
+    }
+  }
+
+  // Datatype
+  if (CountSet({Boolean, Integer, Real}) != 1) {
+    valid = false;
+    if (throw_on_fail) {
+      PARTHENON_THROW("Exactly one data type flag must be set");
+    }
+  }
+
+  // Independent
+  if (CountSet({Independent, Derived}) != 1) {
+    valid = false;
+    if (throw_on_fail) {
+      PARTHENON_THROW("Either the Independent or Derived flag must be set");
+    }
+  }
+
+  // Prolongation/restriction
+  if (IsRefined()) {
+    if (refinement_funcs_.label().size() == 0) {
+      valid = false;
+      if (throw_on_fail) {
+        PARTHENON_THROW(
+            "Registered for refinment but no prolongation/restriction ops found");
+      }
+    }
+  }
+
+  return valid;
+}
+
 std::vector<MetadataFlag> Metadata::Flags() const {
   std::vector<MetadataFlag> set_flags;
   const auto &flags = metadata_state.AllFlags();

src/interface/metadata.hpp

@@ -115,7 +115,11 @@
   /** the variable participate in GMG calculations */                                    \
   PARTHENON_INTERNAL_FOR_FLAG(GMGRestrict)                                               \
   /** the variable must always be allocated for new blocks **/                           \
-  PARTHENON_INTERNAL_FOR_FLAG(ForceAllocOnNewBlocks)
+  PARTHENON_INTERNAL_FOR_FLAG(ForceAllocOnNewBlocks)                                     \
+  /************************************************/                                     \
+  /** Vars specifying coordinates for visualization purposes **/                         \
+  /** You can specify a single 3D var **/                                                \
+  PARTHENON_INTERNAL_FOR_FLAG(CoordinatesVec)
 namespace parthenon {
 
 namespace internal {
@@ -383,62 +387,7 @@ class Metadata {
 
   // Return true if the flags constraints are satisfied, false otherwise. If throw_on_fail
   // is true, throw a descriptive exception when invalid
-  bool IsValid(bool throw_on_fail = false) const {
-    bool valid = true;
-
-    // Topology
-    if (CountSet({None, Node, Edge, Face, Cell}) != 1) {
-      valid = false;
-      if (throw_on_fail) {
-        PARTHENON_THROW("Exactly one topology flag must be set");
-      }
-    }
-
-    // Role
-    if (CountSet({Private, Provides, Requires, Overridable}) != 1) {
-      valid = false;
-      if (throw_on_fail) {
-        PARTHENON_THROW("Exactly one role flag must be set");
-      }
-    }
-
-    // Shape
-    if (CountSet({Vector, Tensor}) > 1) {
-      valid = false;
-      if (throw_on_fail) {
-        PARTHENON_THROW("At most one shape flag can be set");
-      }
-    }
-
-    // Datatype
-    if (CountSet({Boolean, Integer, Real}) != 1) {
-      valid = false;
-      if (throw_on_fail) {
-        PARTHENON_THROW("Exactly one data type flag must be set");
-      }
-    }
-
-    // Independent
-    if (CountSet({Independent, Derived}) != 1) {
-      valid = false;
-      if (throw_on_fail) {
-        PARTHENON_THROW("Either the Independent or Derived flag must be set");
-      }
-    }
-
-    // Prolongation/restriction
-    if (IsRefined()) {
-      if (refinement_funcs_.label().size() == 0) {
-        valid = false;
-        if (throw_on_fail) {
-          PARTHENON_THROW(
-              "Registered for refinment but no prolongation/restriction ops found");
-        }
-      }
-    }
-
-    return valid;
-  }
+  bool IsValid(bool throw_on_fail = false) const;
 
   /*--------------------------------------------------------*/
   // Getters for attributes
@@ -500,6 +449,9 @@ class Metadata {
             IsSet(GMGProlongate) || IsSet(GMGRestrict));
   }
 
+  // Returns true if this variable is a coords var
+  bool IsCoordinateField() const { return IsSet(CoordinatesVec); }
+
   const std::vector<int> &Shape() const { return shape_; }
 
   /*--------------------------------------------------------*/

src/outputs/output_utils.hpp

@@ -55,6 +55,7 @@ struct VarInfo {
   MetadataFlag where;
   bool is_sparse;
   bool is_vector;
+  bool is_coordinate_field;
   IndexShape cellbounds;
   std::vector<std::string> component_labels;
   // list of topological elements in variable... e.g., Face1, Face2, etc
@@ -132,7 +133,8 @@ struct VarInfo {
         tensor_rank(metadata.Shape().size()), where(metadata.Where()),
         topological_elements(topological_elements), is_sparse(is_sparse),
         is_vector(is_vector), cellbounds(cellbounds), rnx_(nx_.rbegin(), nx_.rend()),
-        ntop_elems(topological_elements.size()), element_matters(ntop_elems > 1) {
+        ntop_elems(topological_elements.size()), element_matters(ntop_elems > 1),
+        is_coordinate_field(metadata.IsCoordinateField()) {
     if (num_components <= 0) {
       std::stringstream msg;
       msg << "### ERROR: Got variable " << label << " with " << num_components

src/outputs/parthenon_hdf5.cpp

@@ -241,15 +241,27 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
   // simulation, but not all variables may be allocated on all blocks
 
   auto get_vars = [=](const std::shared_ptr<MeshBlock> pmb) {
-    auto &var_vec = pmb->meshblock_data.Get()->GetVariableVector();
+    const VariableVector<Real> &var_vec = pmb->meshblock_data.Get()->GetVariableVector();
+    VariableVector<Real> coords_vars =
+        GetAnyVariables(var_vec, {parthenon::Metadata::CoordinatesVec});
+    VariableVector<Real> out;
     if (restart_) {
       // get all vars with flag Independent OR restart
-      return GetAnyVariables(
+      out = GetAnyVariables(
           var_vec, {parthenon::Metadata::Independent, parthenon::Metadata::Restart});
     } else {
-      return GetAnyVariables(var_vec, output_params.variables);
+      out = GetAnyVariables(var_vec, output_params.variables);
     }
+    auto coords_loc = std::find_if(out.begin(), out.end(), [](const auto &v) {
+      return v->metadata().IsCoordinateField();
+    });
+    // if we need to add the coords var to the list of output variables, do so.
+    if ((coords_vars.size() == 1) && (coords_loc == out.end())) {
+      out.push_back(coords_vars[0]); // there can be only one
+    }
+    return out;
   };
+
   // get list of all vars, just use the first block since the list is
   // the same for all blocks
   auto all_vars_info = VarInfo::GetAll(get_vars(pm->block_list.front()), cellbounds);