Converted to the new STK simple_fields workflow (#1237)

STK is migrating to a new strategy for registering and managing
Fields, where sizing information is purely specified at run-time
instead of the previous technique of specifying it in a confusing
blend of both compile-time and run-time information.  The compile-time
specification was just a suggestion, as it could be overridden
(possibly inconsistently) at run-time to support variable-length
Fields.  This made it unclear what the true size of a Field was
and where it should be specified.

As an example, registering a vector field on the entire mesh
previously looked like this:

  using VectorField = stk::mesh::Field<double, stk::mesh::Cartesian3d>;
  VectorField & field = meta.declare_field<VectorField>(stk::topology::NODE_RANK, "velocity");
  stk::mesh::put_field_on_mesh(field, meta.universal_part(), 3, nullptr);

and now, it looks like this:

  using VectorField = stk::mesh::Field<double>;
  VectorField & field = meta.declare_field<double>(stk::topology::NODE_RANK, "velocity");
  stk::mesh::put_field_on_mesh(field, meta.universal_part(), 3, nullptr);

  stk::io::set_field_output_type(field, stk::io::FieldOutputType::VECTOR_3D); // Optional

The only template parameter for a Field is now the datatype parameter.
Sizing information now exclusively comes from put_field_on_mesh() calls.
The optional set_field_output_type() function call registers with the
IO sub-system how a multi-component Field should be subscripted in
Exodus files.  If this call is left off, you will get the default
[_1, _2, _3] subscripting.  With the above call, you will instead get
[_x, _y, _z] subscripting.

The MetaData::use_simple_fields() flag is set everywhere possible in
the code to prevent accidental regressions before the old behavior
is formally deprecated and removed.  This will yield a run-time error
if the old-style extra template parameters are used anywhere.  These
calls to use_simple_fields() can be removed in the future once the
STK Mesh back-end has removed support for the old behavior.

This wasn't a completely straightforward conversion due to nalu-wind
making heavy use of various algorithm selections based on the
templated Field type.  The ScalarFieldType, VectorFieldType,
TensorFieldType, and GenericFieldType types are now all identical,
so different techniques had to be used to switch behaviors.

docs/source/theory/codeAbstractions.rst

@@ -259,12 +259,11 @@ elemental, nodal, face and edge of desired size.
       const int numStates = realm_.number_of_states();
 
       // declare it
-      density_ 
-        =  &(metaData_.declare_field<ScalarFieldType>(stk::topology::NODE_RANK, 
-                                                     "density", numStates));
+      density_  = &(metaData_.declare_field<double>(stk::topology::NODE_RANK,
+                                                    "density", numStates));
 
       // put it on this part
-      stk::mesh::put_field(*density_, *part);
+      stk::mesh::put_field_on_mesh(*density_, *part, nullptr);
     }
 
 :math:`\bullet` edge field registration:
@@ -277,9 +276,11 @@ elemental, nodal, face and edge of desired size.
     {
       const int nDim = metaData_.spatial_dimension();
       edgeAreaVec_ 
-        = &(metaData_.declare_field<VectorFieldType>(stk::topology::EDGE_RANK, 
-                                                    "edge_area_vector"));
-      stk::mesh::put_field(*edgeAreaVec_, *part, nDim);
+        = &(metaData_.declare_field<double>(stk::topology::EDGE_RANK,
+                                            "edge_area_vector"));
+      stk::mesh::put_field_on_mesh(*edgeAreaVec_, *part, nDim, nullptr);
+      stk::io::set_field_output_type(*edgeAreaVec_,
+                                     stk::io::FieldOutputType::VECTOR_3D);
     }
 
 :math:`\bullet` side/face field registration:
@@ -297,9 +298,10 @@ elemental, nodal, face and edge of desired size.
         stk::topology::rank_t sideRank 
           = static_cast<stk::topology::rank_t>(metaData_.side_rank());
         GenericFieldType *wallFrictionVelocityBip 
-          =  &(metaData_.declare_field<GenericFieldType>
+          =  &(metaData_.declare_field<double>
               (sideRank, "wall_friction_velocity_bip"));
-        stk::mesh::put_field(*wallFrictionVelocityBip, *part, numIp);
+        stk::mesh::put_field_on_mesh(*wallFrictionVelocityBip, *part, numIp,
+                                     nullptr);
       }
     }
 
@@ -317,23 +319,21 @@ addition to the field template type,
 .. code-block:: c++
 
     VectorFieldType *velocityRTM 
-      = metaData_.get_field<VectorFieldType>(stk::topology::NODE_RANK, 
-                                            "velocity");
+      = metaData_.get_field<double>(stk::topology::NODE_RANK, "velocity");
     ScalarFieldType *density 
-      = metaData_.get_field<ScalarFieldType>(stk::topology::NODE_RANK, 
-                                            "density");}
+      = metaData_.get_field<double>(stk::topology::NODE_RANK, "density");}
 
     VectorFieldType *edgeAreaVec 
-      = metaData_.get_field<VectorFieldType>(stk::topology::EDGE_RANK, 
-                                            "edge_area_vector");
+      = metaData_.get_field<double>(stk::topology::EDGE_RANK,
+                                    "edge_area_vector");
 
-    GenericFieldType  *massFlowRate
-      = metaData_.get_field<GenericFieldType>(stk::topology::ELEMENT_RANK, 
-                                             "mass_flow_rate_scs");
+    GenericFieldType *massFlowRate
+      = metaData_.get_field<double>(stk::topology::ELEMENT_RANK,
+                                    "mass_flow_rate_scs");
 
     GenericFieldType *wallFrictionVelocityBip_ 
-      = metaData_.get_field<GenericFieldType>(metaData_.side_rank(), 
-                                             "wall_friction_velocity_bip");
+      = metaData_.get_field<double>(metaData_.side_rank(),
+                                    "wall_friction_velocity_bip");
 
 State
 ~~~~~
@@ -346,8 +346,7 @@ extracted,
 .. code-block:: c++
 
     ScalarFieldType *density 
-      = metaData_.get_field<ScalarFieldType>(stk::topology::NODE_RANK, 
-                                            "density");
+      = metaData_.get_field<double>(stk::topology::NODE_RANK, "density");
     densityNm1_ = &(density->field_of_state(stk::mesh::StateNM1));
     densityN_ = &(density->field_of_state(stk::mesh::StateN));
     densityNp1_ = &(density->field_of_state(stk::mesh::StateNP1));

include/FieldDefinitions.h

@@ -15,30 +15,33 @@
 namespace sierra {
 namespace nalu {
 
-template <typename T>
+enum class FieldLayout { SCALAR, VECTOR, TENSOR, ARRAY };
+
+template <typename T, FieldLayout Layout = FieldLayout::SCALAR>
 struct FieldDefinition
 {
-  using FieldType = T;
+  using DataType = T;
   const stk::topology::rank_t rank;
   const int num_states{1};
   const int num_components{1};
+  const FieldLayout layout{Layout};
 };
 
-using FieldDefScalar = FieldDefinition<ScalarFieldType>;
-using FieldDefVector = FieldDefinition<VectorFieldType>;
-using FieldDefTensor = FieldDefinition<TensorFieldType>;
-using FieldDefGeneric = FieldDefinition<GenericFieldType>;
-using FieldDefGenericInt = FieldDefinition<GenericIntFieldType>;
-using FieldDefTpetraId = FieldDefinition<TpetIDFieldType>;
-using FieldDefLocalId = FieldDefinition<LocalIdFieldType>;
-using FieldDefGlobalId = FieldDefinition<GlobalIdFieldType>;
-using FieldDefHypreId = FieldDefinition<HypreIDFieldType>;
-using FieldDefScalarInt = FieldDefinition<ScalarIntFieldType>;
+using FieldDefScalar = FieldDefinition<double>;
+using FieldDefVector = FieldDefinition<double, FieldLayout::VECTOR>;
+using FieldDefTensor = FieldDefinition<double, FieldLayout::TENSOR>;
+using FieldDefGeneric = FieldDefinition<double, FieldLayout::ARRAY>;
+using FieldDefGenericInt = FieldDefinition<int, FieldLayout::ARRAY>;
+using FieldDefTpetraId = FieldDefinition<TpetIdType>;
+using FieldDefLocalId = FieldDefinition<LocalId>;
+using FieldDefGlobalId = FieldDefinition<stk::mesh::EntityId>;
+using FieldDefHypreId = FieldDefinition<HypreIntType>;
+using FieldDefScalarInt = FieldDefinition<int>;
 
 // Type redundancy can occur between HypreId and ScalarInt
 // which will break std::variant
 using FieldDefTypes = std::conditional<
-  std::is_same_v<ScalarIntFieldType, HypreIDFieldType>,
+  std::is_same_v<int, HypreIntType>,
   std::variant<
     FieldDefScalar,
     FieldDefVector,
@@ -61,29 +64,22 @@ using FieldDefTypes = std::conditional<
     FieldDefScalarInt,
     FieldDefHypreId>>::type;
 
+// Trouble!
 using FieldPointerTypes = std::conditional<
-  std::is_same_v<ScalarIntFieldType, HypreIDFieldType>,
+  std::is_same_v<int, HypreIntType>,
   std::variant<
-    ScalarFieldType*,
-    VectorFieldType*,
-    TensorFieldType*,
-    GenericFieldType*,
-    GenericIntFieldType*,
-    TpetIDFieldType*,
-    LocalIdFieldType*,
-    GlobalIdFieldType*,
-    ScalarIntFieldType*>,
+    stk::mesh::Field<double>*,
+    stk::mesh::Field<int>*,
+    stk::mesh::Field<LocalId>*,
+    stk::mesh::Field<stk::mesh::EntityId>*,
+    stk::mesh::Field<TpetIdType>*>,
   std::variant<
-    ScalarFieldType*,
-    VectorFieldType*,
-    TensorFieldType*,
-    GenericFieldType*,
-    GenericIntFieldType*,
-    TpetIDFieldType*,
-    LocalIdFieldType*,
-    GlobalIdFieldType*,
-    ScalarIntFieldType*,
-    HypreIDFieldType*>>::type;
+    stk::mesh::Field<double>*,
+    stk::mesh::Field<int>*,
+    stk::mesh::Field<LocalId>*,
+    stk::mesh::Field<stk::mesh::EntityId>*,
+    stk::mesh::Field<TpetIdType>*,
+    stk::mesh::Field<HypreIntType>*>>::type;
 
 } // namespace nalu
 } // namespace sierra

include/FieldManager.h

@@ -57,7 +57,7 @@ class FieldManager
   /// ScalarFieldType* and so on.
   ///
   template <typename T>
-  T* register_field(
+  stk::mesh::Field<T>* register_field(
     const std::string& name,
     const stk::mesh::PartVector& parts,
     const void* init_val = nullptr,
@@ -86,28 +86,26 @@ class FieldManager
     stk::mesh::FieldState state = stk::mesh::FieldState::StateNone)
   {
     register_field(name, parts, numStates, numComponents, init_val);
-    return get_field_ptr<GenericFieldType>(name, state);
+    return get_field_ptr<GenericFieldType::value_type>(name, state);
   }
 
   // Return a field by the given name and of type T, the template parameter.
   // This function will throw if the named field is not of the type
   // specified by the template parameter: ScalarFieldType, VectorFieldType,
   // ScalarIntFieldType, GlobalIdFieldType,....
   template <typename T>
-  T* get_field_ptr(
+  stk::mesh::Field<T>* get_field_ptr(
     const std::string& name,
     stk::mesh::FieldState state = stk::mesh::FieldState::StateNone) const
   {
     FieldDefTypes fieldDef =
       FieldRegistry::query(numDimensions_, numStates_, name);
     FieldPointerTypes pointerSet = std::visit(
       [&](auto def) -> FieldPointerTypes {
-        return &meta_
-                  .get_field<typename decltype(def)::FieldType>(def.rank, name)
-                  ->field_of_state(state);
+        return &meta_.get_field<T>(def.rank, name)->field_of_state(state);
       },
       fieldDef);
-    return std::get<T*>(pointerSet);
+    return std::get<stk::mesh::Field<T>*>(pointerSet);
   }
 
   /// Register a field with the option to override default parameters that
@@ -139,9 +137,7 @@ class FieldManager
       FieldRegistry::query(numDimensions_, numStates_, name);
     const stk::mesh::FieldBase& stkField = std::visit(
       [&](auto def) -> stk::mesh::FieldBase& {
-        return meta_
-          .get_field<typename decltype(def)::FieldType>(def.rank, name)
-          ->field_of_state(state);
+        return meta_.get_field<T>(def.rank, name)->field_of_state(state);
       },
       fieldDef);
     stk::mesh::NgpField<T>& tmp = stk::mesh::get_updated_ngp_field<T>(stkField);
@@ -154,16 +150,16 @@ class FieldManager
     std::string name,
     stk::mesh::FieldState state = stk::mesh::FieldState::StateNone) const
   {
-    return MakeSmartField<tags::DEVICE, ACCESS>().template operator()<T>(
+    return MakeSmartField<tags::DEVICE, ACCESS>().template operator()(
       get_ngp_field_ptr<T>(name, state));
   }
 
   template <typename T, typename ACCESS>
-  SmartField<T, tags::LEGACY, ACCESS> get_legacy_smart_field(
+  SmartField<stk::mesh::Field<T>, tags::LEGACY, ACCESS> get_legacy_smart_field(
     std::string name,
     stk::mesh::FieldState state = stk::mesh::FieldState::StateNone) const
   {
-    return MakeSmartField<tags::LEGACY, ACCESS>().template operator()<T>(
+    return MakeSmartField<tags::LEGACY, ACCESS>().template operator()(
       get_field_ptr<T>(name, state));
   }
 };

include/FieldTypeDef.h

@@ -11,7 +11,6 @@
 #define FieldTypeDef_h
 
 #include <stk_mesh/base/FieldBase.hpp>
-#include <stk_mesh/base/CoordinateSystems.hpp>
 #include <stk_mesh/base/Ngp.hpp>
 #include <stk_mesh/base/NgpField.hpp>
 
@@ -27,43 +26,38 @@
 namespace sierra {
 namespace nalu {
 
-// define scalar field typedef
-typedef stk::mesh::Field<double> ScalarFieldType;
-typedef stk::mesh::Field<stk::mesh::EntityId> GlobalIdFieldType;
-typedef stk::mesh::Field<int> ScalarIntFieldType;
-typedef stk::mesh::NgpField<double> NGPDoubleFieldType;
-typedef stk::mesh::NgpField<stk::mesh::EntityId> NGPGlobalIdFieldType;
-typedef stk::mesh::NgpField<int> NGPScalarIntFieldType;
+using ScalarFieldType = stk::mesh::Field<double>;
+using GlobalIdFieldType = stk::mesh::Field<stk::mesh::EntityId>;
+using ScalarIntFieldType = stk::mesh::Field<int>;
+using NGPDoubleFieldType = stk::mesh::NgpField<double>;
+using NGPGlobalIdFieldType = stk::mesh::NgpField<stk::mesh::EntityId>;
+using NGPScalarIntFieldType = stk::mesh::NgpField<int>;
 
-// define vector field typedef; however, what is the value of Cartesian?
-typedef stk::mesh::Field<double, stk::mesh::Cartesian> VectorFieldType;
-typedef stk::mesh::Field<double, stk::mesh::FullTensor> TensorFieldType;
+using VectorFieldType = stk::mesh::Field<double>;
+using TensorFieldType = stk::mesh::Field<double>;
 
-// define generic
-typedef stk::mesh::Field<double, stk::mesh::SimpleArrayTag> GenericFieldType;
+using GenericFieldType = stk::mesh::Field<double>;
 
-typedef stk::mesh::Field<int, stk::mesh::SimpleArrayTag> GenericIntFieldType;
+using GenericIntFieldType = stk::mesh::Field<int>;
 
-// field type for local ids
-typedef unsigned LocalId;
-typedef stk::mesh::Field<LocalId> LocalIdFieldType;
+using LocalId = unsigned;
+using LocalIdFieldType = stk::mesh::Field<LocalId>;
 
-// Hypre Integer types
 #ifdef NALU_USES_HYPRE
-typedef HYPRE_Int HypreIntType;
+using HypreIntType = HYPRE_Int;
 #else
-typedef int HypreIntType;
+using HypreIntType = int;
 #endif
 
 #ifdef NALU_USES_TRILINOS_SOLVERS
-typedef stk::mesh::Field<Tpetra::Details::DefaultTypes::global_ordinal_type>
-  TpetIDFieldType;
+using TpetIdType = Tpetra::Details::DefaultTypes::global_ordinal_type;
 #else
-typedef stk::mesh::Field<int64_t> TpetIDFieldType;
+using TpetIdType = int64_t;
 #endif
 
-typedef stk::mesh::Field<HypreIntType> HypreIDFieldType;
-typedef stk::mesh::NgpField<HypreIntType> NGPHypreIDFieldType;
+using TpetIDFieldType = stk::mesh::Field<TpetIdType>;
+using HypreIDFieldType = stk::mesh::Field<HypreIntType>;
+using NGPHypreIDFieldType = stk::mesh::NgpField<HypreIntType>;
 
 } // namespace nalu
 } // namespace sierra

include/MatrixFreeHeatCondEquationSystem.h

@@ -16,7 +16,6 @@
 #include "Kokkos_Array.hpp"
 
 #include "stk_mesh/base/Ngp.hpp"
-#include "stk_mesh/base/CoordinateSystems.hpp"
 
 #include "Realm.h"
 

include/ProjectedNodalGradientEquationSystem.h

@@ -16,7 +16,6 @@
 #include <NaluParsedTypes.h>
 
 #include <stk_mesh/base/FieldBase.hpp>
-#include <stk_mesh/base/CoordinateSystems.hpp>
 
 namespace stk {
 struct topology;

include/SmartField.h

@@ -401,9 +401,10 @@ struct MakeSmartField<LEGACY, ACCESS>
 {
   // use pointer since that is the common access type for stk::mesh::Field<T>
   template <typename T>
-  SmartField<T, LEGACY, ACCESS> operator()(T* field)
+  SmartField<stk::mesh::Field<T>, LEGACY, ACCESS>
+  operator()(stk::mesh::Field<T>* field)
   {
-    return SmartField<T, LEGACY, ACCESS>(*field);
+    return SmartField<stk::mesh::Field<T>, LEGACY, ACCESS>(*field);
   }
 };