Example usage of FieldManager in an NGP Kernel

include/FieldManager.h

@@ -15,7 +15,6 @@
 #include <stk_mesh/base/FieldState.hpp>
 #include "stk_mesh/base/GetNgpField.hpp"
 #include <string>
-#include <type_traits>
 
 namespace stk::mesh {
 class MetaData;
@@ -60,12 +59,12 @@ class FieldManager
   stk::mesh::Field<T>* register_field(
     const std::string& name,
     const stk::mesh::PartVector& parts,
-    const void* init_val = nullptr,
+    const void* initVal = nullptr,
     stk::mesh::FieldState state = stk::mesh::FieldState::StateNone) const
   {
     const int numStates = 0;
     const int numComponents = 0;
-    register_field(name, parts, numStates, numComponents, init_val);
+    register_field(name, parts, numStates, numComponents, initVal);
     return get_field_ptr<T>(name, state);
   }
 
@@ -82,10 +81,10 @@ class FieldManager
     const stk::mesh::PartVector& parts,
     const int numStates,
     const int numComponents,
-    const void* init_val = nullptr,
+    const void* initVal = nullptr,
     stk::mesh::FieldState state = stk::mesh::FieldState::StateNone) const
   {
-    register_field(name, parts, numStates, numComponents, init_val);
+    register_field(name, parts, numStates, numComponents, initVal);
     return get_field_ptr<GenericFieldType::value_type>(name, state);
   }
 
@@ -112,7 +111,7 @@ class FieldManager
   /// would otherwise be defined in the field Registry class.
   ///
   /// If numStates = 0 then the number of states comes from the
-  /// field Registry.  Same for numComponents = 0 and init_val = nullptr.
+  /// field Registry.  Same for numComponents = 0 and initVal = nullptr.
   ///
   /// This is useful for dynamic fields that depend on the input
   /// options to define the number of states or number of components since the
@@ -124,7 +123,8 @@ class FieldManager
     const stk::mesh::PartVector& parts,
     const int numStates = 0,
     const int numComponents = 0,
-    const void* init_val = nullptr) const;
+    const void* initVal = nullptr) const;
+
 
   /// Given the named field that has already been registered on the CPU
   /// return the GPU version of the same field.

include/SmartField.h

@@ -445,11 +445,11 @@ struct MakeSmartField<DEVICE, ACCESS>
 };
 
 template <typename T, typename ACCESS>
-using SmartDeviceField = SmartField<T, DEVICE, ACCESS>;
+using SmartDeviceField = SmartField<stk::mesh::NgpField<T>, DEVICE, ACCESS>;
 template <typename T, typename ACCESS>
-using SmartHostField = SmartField<T, HOST, ACCESS>;
+using SmartHostField = SmartField<stk::mesh::HostField<T>, HOST, ACCESS>;
 template <typename T, typename ACCESS>
-using SmartLegacyField = SmartField<T, LEGACY, ACCESS>;
+using SmartLegacyField = SmartField<stk::mesh::Field<T>, LEGACY, ACCESS>;
 } // namespace sierra::nalu
 
 #endif

include/node_kernels/TKESSTNodeKernel.h

@@ -10,12 +10,14 @@
 #ifndef TKESSTNODEKERNEL_H
 #define TKESSTNODEKERNEL_H
 
+#include "LinearSystem.h"
 #include "node_kernels/NodeKernel.h"
 #include "FieldTypeDef.h"
 #include "stk_mesh/base/BulkData.hpp"
 #include "stk_mesh/base/Ngp.hpp"
 #include "stk_mesh/base/NgpField.hpp"
 #include "stk_mesh/base/Types.hpp"
+#include <FieldManager.h>
 
 namespace sierra {
 namespace nalu {
@@ -25,7 +27,8 @@ class Realm;
 class TKESSTNodeKernel : public NGPNodeKernel<TKESSTNodeKernel>
 {
 public:
-  TKESSTNodeKernel(const stk::mesh::MetaData&);
+  TKESSTNodeKernel(
+    const stk::mesh::MetaData&, const FieldManager&, stk::mesh::PartVector& parts);
 
   TKESSTNodeKernel() = delete;
 
@@ -41,20 +44,14 @@ class TKESSTNodeKernel : public NGPNodeKernel<TKESSTNodeKernel>
     const stk::mesh::FastMeshIndex&) override;
 
 private:
+  // TODO can we make these SmartFields? Do we need to? review how kernels manage syncs
   stk::mesh::NgpField<double> tke_;
   stk::mesh::NgpField<double> sdr_;
   stk::mesh::NgpField<double> density_;
   stk::mesh::NgpField<double> tvisc_;
   stk::mesh::NgpField<double> dudx_;
   stk::mesh::NgpField<double> dualNodalVolume_;
 
-  unsigned tkeID_{stk::mesh::InvalidOrdinal};
-  unsigned sdrID_{stk::mesh::InvalidOrdinal};
-  unsigned densityID_{stk::mesh::InvalidOrdinal};
-  unsigned tviscID_{stk::mesh::InvalidOrdinal};
-  unsigned dudxID_{stk::mesh::InvalidOrdinal};
-  unsigned dualNodalVolumeID_{stk::mesh::InvalidOrdinal};
-
   NodeKernelTraits::DblType betaStar_;
   NodeKernelTraits::DblType tkeProdLimitRatio_;
   NodeKernelTraits::DblType tkeAmb_;

src/FieldManager.C

@@ -37,7 +37,7 @@ FieldManager::register_field(
   const stk::mesh::PartVector& parts,
   const int numStates,
   const int numComponents,
-  const void* init_val) const
+  const void* initVal) const
 {
   auto definition = FieldRegistry::query(numDimensions_, numStates_, name);
 
@@ -49,8 +49,9 @@ FieldManager::register_field(
         numComponents ? numComponents : def.num_components;
       const FieldLayout layout = def.layout;
 
-      const val_type* init = static_cast<const val_type*>(init_val);
+      const val_type* init = static_cast<const val_type*>(initVal);
       auto* id = &(meta_.declare_field<val_type>(def.rank, name, num_states));
+
       for (auto&& part : parts) {
         stk::mesh::put_field_on_mesh(*id, *part, num_components, init);
 

src/FieldRegistry.C

@@ -46,16 +46,17 @@ Registry()
     {"coordinates",               SingleStateNodalVector},
     {"coordinates_copy",          SingleStateNodalVector}, // Used in testing
     {"current_coordinates",       SingleStateNodalVector},
-    {"density",                   SingleStateNodalScalar},
+    {"density",                   MultiStateNodalScalar},
     {"dhdx",                      SingleStateNodalVector},
     {"diffFluxCoeff",             SingleStateNodalScalar},
     {"discreteLaplacian",         SingleStateNodalScalar},
     {"div_mesh_velocity",         SingleStateNodalScalar},
     {"dkdx",                      SingleStateNodalVector},
     {"dpdx",                      SingleStateNodalVector},
-    {"dual_nodal_volume",         SingleStateNodalScalar},
+    {"dual_nodal_volume",         MultiStateNodalScalar},
     {"dudx",                      SingleStateNodalTensor},
     {"dwdx",                      SingleStateNodalVector},
+    {"dynamic_pressure",          SingleStateNodalScalar}, // Used in testing
     {"edge_area_vector",          SingleStateEdgeVector},
     {"effective_viscosity" ,      SingleStateNodalScalar},
     {"elemCentroid",              SingleStateElemVector},
@@ -83,7 +84,7 @@ Registry()
     {"pressure",                  SingleStateNodalScalar},
     {"rans_time_scale" ,          SingleStateNodalScalar},
     {"scalarQ",                   SingleStateNodalScalar},
-    {"specific_dissipation_rate", MultiStateNodalScalar},
+    {"specific_dissipation_rate", SingleStateNodalScalar},
     {"thermal_conductivity",      SingleStateNodalScalar},
     {"specific_heat" ,            SingleStateNodalScalar},
     {"sst_f_one_blending"  ,      SingleStateNodalScalar},
@@ -94,6 +95,7 @@ Registry()
     {"turbulent_viscosity" ,      SingleStateNodalScalar},
     {"velocity",                  MultiStateNodalVector},
     {"velocity_rtm",              SingleStateNodalVector},
+    {"velocity_bc",               MultiStateNodalVector},  // Used in testing
     {"viscosity",                 SingleStateNodalScalar}
   };
   // clang-format on

src/Simulation.C

@@ -173,9 +173,12 @@ Simulation::setSerializedIOGroupSize(int siogs)
 void
 Simulation::breadboard()
 {
-  realms_->breadboard();
-  timeIntegrator_->breadboard();
-  transfers_->breadboard();
+  if (realms_)
+    realms_->breadboard();
+  if (timeIntegrator_)
+    timeIntegrator_->breadboard();
+  if (transfers_)
+    transfers_->breadboard();
 }
 
 void

src/TurbKineticEnergyEquationSystem.C

@@ -298,12 +298,18 @@ TurbKineticEnergyEquationSystem::register_interior_algorithm(
         if (!elementMassAlg)
           nodeAlg.add_kernel<ScalarMassBDFNodeKernel>(realm_.bulk_data(), tke_);
 
+        // hacky solution since updating everything to part vecs at once would be painful and error prone
+        stk::mesh::PartVector tempVec({part});
+
         switch (turbulenceModel_) {
         case TurbulenceModel::KSGS:
           nodeAlg.add_kernel<TKEKsgsNodeKernel>(realm_.meta_data());
           break;
         case TurbulenceModel::SST:
-          nodeAlg.add_kernel<TKESSTNodeKernel>(realm_.meta_data());
+          // should be able to keep creating actually since stk fields can be registered over and over as no-ops
+          // registration should still happen
+          nodeAlg.add_kernel<TKESSTNodeKernel>(
+            realm_.meta_data(), *(realm_.fieldManager_.get()), tempVec);
           break;
         case TurbulenceModel::SSTLR:
           nodeAlg.add_kernel<TKESSTLRNodeKernel>(realm_.meta_data());

src/node_kernels/TKESSTNodeKernel.C

@@ -19,29 +19,31 @@
 namespace sierra {
 namespace nalu {
 
-TKESSTNodeKernel::TKESSTNodeKernel(const stk::mesh::MetaData& meta)
-  : NGPNodeKernel<TKESSTNodeKernel>(),
-    tkeID_(get_field_ordinal(meta, "turbulent_ke")),
-    sdrID_(get_field_ordinal(meta, "specific_dissipation_rate")),
-    densityID_(get_field_ordinal(meta, "density")),
-    tviscID_(get_field_ordinal(meta, "turbulent_viscosity")),
-    dudxID_(get_field_ordinal(meta, "dudx")),
-    dualNodalVolumeID_(get_field_ordinal(meta, "dual_nodal_volume")),
-    nDim_(meta.spatial_dimension())
+TKESSTNodeKernel::TKESSTNodeKernel(
+  const stk::mesh::MetaData& meta,
+  const FieldManager& manager,
+  stk::mesh::PartVector& parts)
+  : NGPNodeKernel<TKESSTNodeKernel>(), nDim_(meta.spatial_dimension())
 {
+  manager.register_field("turbulent_ke", parts);
+  manager.register_field("specific_dissipation_rate", parts);
+  manager.register_field("density", parts);
+  manager.register_field("turbulent_viscosity", parts);
+  manager.register_field("dudx", parts);
+  manager.register_field("dual_nodal_volume", parts);
 }
 
 void
 TKESSTNodeKernel::setup(Realm& realm)
 {
-  const auto& fieldMgr = realm.ngp_field_manager();
+  const auto& fieldMgr = *(realm.fieldManager_.get());
 
-  tke_ = fieldMgr.get_field<double>(tkeID_);
-  sdr_ = fieldMgr.get_field<double>(sdrID_);
-  density_ = fieldMgr.get_field<double>(densityID_);
-  tvisc_ = fieldMgr.get_field<double>(tviscID_);
-  dudx_ = fieldMgr.get_field<double>(dudxID_);
-  dualNodalVolume_ = fieldMgr.get_field<double>(dualNodalVolumeID_);
+  tke_ = fieldMgr.get_ngp_field_ptr<double>("turbulent_ke");
+  sdr_ = fieldMgr.get_ngp_field_ptr<double>("specific_dissipation_rate");
+  density_ = fieldMgr.get_ngp_field_ptr<double>("density");
+  tvisc_ = fieldMgr.get_ngp_field_ptr<double>("turbulent_viscosity");
+  dudx_ = fieldMgr.get_ngp_field_ptr<double>("dudx");
+  dualNodalVolume_ = fieldMgr.get_ngp_field_ptr<double>("dual_nodal_volume");
 
   // Update turbulence model constants
   betaStar_ = realm.get_turb_model_constant(TM_betaStar);

unit_tests/UnitTestHelperObjects.h

@@ -16,6 +16,7 @@
 #include <stk_topology/topology.hpp>
 
 #include <memory>
+#include <FieldManager.h>
 
 namespace unit_test_utils {
 
@@ -33,6 +34,12 @@ struct HelperObjectsBase
       eqSystem(eqSystems)
   {
     realm.bulkData_ = bulk;
+    // TODO fix this! realm should not be getting time integrator this way!!
+    naluObj->sim_.breadboard();
+    assert(naluObj->sim_.timeIntegrator_ != NULL);
+    // bread board didn't work so set integrator manually
+    realm.timeIntegrator_ = naluObj->sim_.timeIntegrator_;
+    realm.setup_field_manager();
   }
 
   virtual ~HelperObjectsBase() { delete naluObj; }

unit_tests/kernels/UnitTestKernelUtils.h

@@ -256,40 +256,26 @@ class TestKernelHex8Mesh : public ::testing::Test
   TestKernelHex8Mesh()
     : comm_(MPI_COMM_WORLD), spatialDim_(3), solnOpts_(), coordinates_(nullptr)
   {
+    const int numStates = 3;
     stk::mesh::MeshBuilder meshBuilder(MPI_COMM_WORLD);
     meshBuilder.set_spatial_dimension(spatialDim_);
     bulk_ = meshBuilder.create();
     meta_ = &bulk_->mesh_meta_data();
     meta_->use_simple_fields();
-
-    naluGlobalId_ = &meta_->declare_field<stk::mesh::EntityId>(
-      stk::topology::NODE_RANK, "nalu_global_id", 1);
-    tpetGlobalId_ = &meta_->declare_field<sierra::nalu::TpetIdType>(
-      stk::topology::NODE_RANK, "tpet_global_id", 1);
-    dnvField_ = &meta_->declare_field<double>(
-      stk::topology::NODE_RANK, "dual_nodal_volume", 3);
-    divMeshVelField_ = &meta_->declare_field<double>(
-      stk::topology::NODE_RANK, "div_mesh_velocity");
-    edgeAreaVec_ = &meta_->declare_field<double>(
-      stk::topology::EDGE_RANK, "edge_area_vector");
-    elementVolume_ =
-      &meta_->declare_field<double>(stk::topology::ELEM_RANK, "element_volume");
+    fieldManager =
+      std::make_shared<sierra::nalu::FieldManager>(*meta_, numStates);
+
+    const stk::mesh::PartVector parts(1, &meta_->universal_part());
+
+    naluGlobalId_ = fieldManager->register_field<stk::mesh::EntityId>("nalu_global_id", parts);
+    tpetGlobalId_ = fieldManager->register_field<sierra::nalu::TpetIdType>("tpet_global_id", parts);
+    dnvField_ = fieldManager->register_field<double>("dual_nodal_volume", parts);
+    divMeshVelField_ = fieldManager->register_field<double>("div_mesh_velocity", parts);
+    edgeAreaVec_ = fieldManager->register_field<double>("edge_area_vector", parts);
+    elementVolume_ = fieldManager->register_field<double>("element_volume", parts);
+    // currently don't handle side ranks
     exposedAreaVec_ =
       &meta_->declare_field<double>(meta_->side_rank(), "exposed_area_vector");
-
-    stk::mesh::put_field_on_mesh(
-      *naluGlobalId_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(
-      *tpetGlobalId_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(*dnvField_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(
-      *divMeshVelField_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(
-      *edgeAreaVec_, meta_->universal_part(), spatialDim_, nullptr);
-    stk::io::set_field_output_type(
-      *edgeAreaVec_, stk::io::FieldOutputType::VECTOR_3D);
-    stk::mesh::put_field_on_mesh(
-      *elementVolume_, meta_->universal_part(), nullptr);
     stk::mesh::put_field_on_mesh(
       *exposedAreaVec_, meta_->universal_part(),
       spatialDim_ * sierra::nalu::AlgTraitsQuad4::numScsIp_, nullptr);
@@ -329,6 +315,7 @@ class TestKernelHex8Mesh : public ::testing::Test
   unsigned spatialDim_;
   stk::mesh::MetaData* meta_;
   std::shared_ptr<stk::mesh::BulkData> bulk_;
+  std::shared_ptr<sierra::nalu::FieldManager> fieldManager;
   stk::mesh::PartVector partVec_;
 
   sierra::nalu::SolutionOptions solnOpts_;
@@ -361,38 +348,21 @@ class LowMachKernelHex8Mesh : public TestKernelHex8Mesh
 {
 public:
   LowMachKernelHex8Mesh()
-    : TestKernelHex8Mesh(),
-      velocity_(
-        &meta_->declare_field<double>(stk::topology::NODE_RANK, "velocity", 2)),
-      dpdx_(&meta_->declare_field<double>(stk::topology::NODE_RANK, "dpdx", 2)),
-      density_(
-        &meta_->declare_field<double>(stk::topology::NODE_RANK, "density", 2)),
-      pressure_(
-        &meta_->declare_field<double>(stk::topology::NODE_RANK, "pressure", 2)),
-      Udiag_(&meta_->declare_field<double>(
-        stk::topology::NODE_RANK, "momentum_diag", 2)),
-      velocityBC_(
-        &meta_->declare_field<double>(stk::topology::NODE_RANK, "velocity_bc")),
-      dynP_(
-        &meta_->declare_field<double>(meta_->side_rank(), "dynamic_pressure"))
+    : TestKernelHex8Mesh()
   {
-    stk::mesh::put_field_on_mesh(
-      *velocity_, meta_->universal_part(), spatialDim_, nullptr);
-    stk::io::set_field_output_type(
-      *velocity_, stk::io::FieldOutputType::VECTOR_3D);
-    stk::mesh::put_field_on_mesh(
-      *dpdx_, meta_->universal_part(), spatialDim_, nullptr);
-    stk::io::set_field_output_type(*dpdx_, stk::io::FieldOutputType::VECTOR_3D);
-    stk::mesh::put_field_on_mesh(*density_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(*pressure_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(*Udiag_, meta_->universal_part(), nullptr);
-    stk::mesh::put_field_on_mesh(
-      *velocityBC_, meta_->universal_part(), spatialDim_, nullptr);
-    stk::io::set_field_output_type(
-      *velocityBC_, stk::io::FieldOutputType::VECTOR_3D);
-    stk::mesh::put_field_on_mesh(
-      *dynP_, meta_->universal_part(), sierra::nalu::AlgTraitsQuad4::numScsIp_,
-      nullptr);
+      const stk::mesh::PartVector parts(1, &meta_->universal_part());
+      velocity_ = fieldManager->register_field<double>("velocity", parts);
+      dpdx_ = fieldManager->register_field<double>("dpdx", parts);
+      density_ = fieldManager->register_field<double>("density", parts);
+      pressure_ = fieldManager->register_field<double>("pressure", parts);
+      Udiag_ = fieldManager->register_field<double>("momentum_diag", parts);
+      velocityBC_ = fieldManager->register_field<double>("velocity_bc", parts);
+
+      // currently don't handle side ranks
+      dynP_ =
+        &meta_->declare_field<double>(meta_->side_rank(), "dynamic_pressure");
+      stk::mesh::put_field_on_mesh(
+        *dynP_, meta_->universal_part(), sierra::nalu::AlgTraitsQuad4::numScsIp_, nullptr);
   }
 
   virtual ~LowMachKernelHex8Mesh() {}