diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0507404b6..1e37789af 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -11,13 +11,13 @@ jobs:
     runs-on: ubuntu-latest
     steps:
     - name: Clone
-      uses: actions/checkout@v3
+      uses: actions/checkout@v4
     - name: Check formatting
-      uses: DoozyX/clang-format-lint-action@v0.13
+      uses: DoozyX/clang-format-lint-action@v0.17
       with:
         source: 'nalu.C unit_tests.C ./include ./src ./unit_tests'
         extensions: 'H,h,cpp,C'
-        clangFormatVersion: 13
+        clangFormatVersion: 17
   CPU-Trilinos:
     needs: Formatting
     runs-on: ubuntu-latest
diff --git a/include/Algorithm.h b/include/Algorithm.h
index 4b19861ca..eb76e2153 100644
--- a/include/Algorithm.h
+++ b/include/Algorithm.h
@@ -25,6 +25,7 @@ class Realm;
 class MasterElement;
 class SupplementalAlgorithm;
 class Kernel;
+class FieldManager;
 
 class Algorithm
 {
@@ -43,6 +44,7 @@ class Algorithm
 
   Realm& realm_;
   stk::mesh::PartVector partVec_;
+  const FieldManager& fieldManager_;
   std::vector<SupplementalAlgorithm*> supplementalAlg_;
 
   std::vector<Kernel*> activeKernels_;
diff --git a/include/AssembleElemSolverAlgorithm.h b/include/AssembleElemSolverAlgorithm.h
index 3cf9ee589..c4ba55daa 100644
--- a/include/AssembleElemSolverAlgorithm.h
+++ b/include/AssembleElemSolverAlgorithm.h
@@ -87,7 +87,7 @@ class AssembleElemSolverAlgorithm : public SolverAlgorithm
         auto& b = ngpMesh.get_bucket(entityRank, bktId);
 
 #if !defined(KOKKOS_ENABLE_GPU)
-        ThrowAssertMsg(
+        STK_ThrowAssertMsg(
           b.topology().num_nodes() == (unsigned)nodesPerEntity_,
           "AssembleElemSolverAlgorithm expected nodesPerEntity_ = "
             << nodesPerEntity_
diff --git a/include/AssembleFaceElemSolverAlgorithm.h b/include/AssembleFaceElemSolverAlgorithm.h
index fbd9715db..2340e26ea 100644
--- a/include/AssembleFaceElemSolverAlgorithm.h
+++ b/include/AssembleFaceElemSolverAlgorithm.h
@@ -83,7 +83,7 @@ class AssembleFaceElemSolverAlgorithm : public SolverAlgorithm
         auto& b = ngpMesh.get_bucket(sideRank, bktId);
 
 #if !defined(KOKKOS_ENABLE_GPU)
-        ThrowAssertMsg(
+        STK_ThrowAssertMsg(
           b.topology().num_nodes() == (unsigned)nodesPerFace_,
           "AssembleFaceElemSolverAlgorithm expected nodesPerEntity_ = "
             << nodesPerFace_
@@ -111,7 +111,7 @@ class AssembleFaceElemSolverAlgorithm : public SolverAlgorithm
                 stk::mesh::Entity face =
                   b[bktIndex * simdLen + numFacesProcessed + simdFaceIndex];
                 const auto ngpFaceIndex = ngpMesh.fast_mesh_index(face);
-                // ThrowAssertMsg(
+                // STK_ThrowAssertMsg(
                 //   bulk.num_elements(face) == 1,
                 //   "Expecting just 1 element attached to face!");
                 int thisElemFaceOrdinal =
diff --git a/include/CopyAndInterleave.h b/include/CopyAndInterleave.h
index db656eefb..5905d937f 100644
--- a/include/CopyAndInterleave.h
+++ b/include/CopyAndInterleave.h
@@ -59,11 +59,11 @@ copy_and_interleave(
   for (unsigned viewIndex = 0; viewIndex < numViews; ++viewIndex) {
     for (int simdIndex = 0; simdIndex < simdElems; ++simdIndex) {
       src[simdIndex] = data[simdIndex]->get_1D_view_by_index(viewIndex).data();
-      NGP_ThrowAssert(
+      STK_NGP_ThrowAssert(
         data[simdIndex]->get_1D_view_by_index(viewIndex).size() ==
         simdData.get_1D_view_by_index(viewIndex).size());
-      NGP_ThrowAssert(src[simdIndex] != nullptr);
-      NGP_ThrowAssert(
+      STK_NGP_ThrowAssert(src[simdIndex] != nullptr);
+      STK_NGP_ThrowAssert(
         src[simdIndex][0] ==
         data[simdIndex]->get_1D_view_by_index(viewIndex).data()[0]);
     }
@@ -74,11 +74,11 @@ copy_and_interleave(
   for (unsigned viewIndex = 0; viewIndex < numViews; ++viewIndex) {
     for (int simdIndex = 0; simdIndex < simdElems; ++simdIndex) {
       src[simdIndex] = data[simdIndex]->get_2D_view_by_index(viewIndex).data();
-      NGP_ThrowAssert(
+      STK_NGP_ThrowAssert(
         data[simdIndex]->get_2D_view_by_index(viewIndex).size() ==
         simdData.get_2D_view_by_index(viewIndex).size());
-      NGP_ThrowAssert(src[simdIndex] != nullptr);
-      NGP_ThrowAssert(
+      STK_NGP_ThrowAssert(src[simdIndex] != nullptr);
+      STK_NGP_ThrowAssert(
         src[simdIndex][0] ==
         data[simdIndex]->get_2D_view_by_index(viewIndex).data()[0]);
     }
@@ -89,11 +89,11 @@ copy_and_interleave(
   for (unsigned viewIndex = 0; viewIndex < numViews; ++viewIndex) {
     for (int simdIndex = 0; simdIndex < simdElems; ++simdIndex) {
       src[simdIndex] = data[simdIndex]->get_3D_view_by_index(viewIndex).data();
-      NGP_ThrowAssert(
+      STK_NGP_ThrowAssert(
         data[simdIndex]->get_3D_view_by_index(viewIndex).size() ==
         simdData.get_3D_view_by_index(viewIndex).size());
-      NGP_ThrowAssert(src[simdIndex] != nullptr);
-      NGP_ThrowAssert(
+      STK_NGP_ThrowAssert(src[simdIndex] != nullptr);
+      STK_NGP_ThrowAssert(
         src[simdIndex][0] ==
         data[simdIndex]->get_3D_view_by_index(viewIndex).data()[0]);
     }
diff --git a/include/ElemDataRequests.h b/include/ElemDataRequests.h
index 89de49ee7..6e2cbed12 100644
--- a/include/ElemDataRequests.h
+++ b/include/ElemDataRequests.h
@@ -217,7 +217,7 @@ class ElemDataRequests
   get_coordinates_field(const COORDS_TYPES cType) const
   {
     auto it = coordsFields_.find(cType);
-    NGP_ThrowRequireMsg(
+    STK_NGP_ThrowRequireMsg(
       it != coordsFields_.end(),
       "ElemDataRequests:get_coordinates_field: Coordinates field "
       "must be registered to ElemDataRequests before access");
diff --git a/include/Enums.h b/include/Enums.h
index e54bc29cc..b12022330 100644
--- a/include/Enums.h
+++ b/include/Enums.h
@@ -105,6 +105,7 @@ static const std::string EquationTypeMap[] = {
   "PNG_U",
   "PNG_TKE",
   "Wall_Distance",
+  "Gamma_Transition",
   "Volume_of_Fluid"};
 
 enum UserDataType {
diff --git a/include/FieldManager.h b/include/FieldManager.h
index c510d5a6a..030bbda69 100644
--- a/include/FieldManager.h
+++ b/include/FieldManager.h
@@ -60,7 +60,7 @@ class FieldManager
     const std::string& name,
     const stk::mesh::PartVector& parts,
     const void* initVal = nullptr,
-    stk::mesh::FieldState state = stk::mesh::FieldState::StateNone)
+    stk::mesh::FieldState state = stk::mesh::FieldState::StateNone) const
   {
     const int numStates = 0;
     const int numComponents = 0;
@@ -69,7 +69,7 @@ class FieldManager
   }
 
   /// Check to see if the field has been registered.
-  bool field_exists(const std::string& name);
+  bool field_exists(const std::string& name) const;
 
   unsigned size() const { return meta_.get_fields().size(); }
   /// Register a Generic field.
@@ -82,9 +82,9 @@ class FieldManager
     const int numStates,
     const int numComponents,
     const void* initVal = nullptr,
-    stk::mesh::FieldState state = stk::mesh::FieldState::StateNone)
+    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);
   }
 
@@ -123,19 +123,13 @@ class FieldManager
     const stk::mesh::PartVector& parts,
     const int numStates = 0,
     const int numComponents = 0,
-    const void* initVal = nullptr) const;
-
-  FieldPointerTypes register_field(
-    const std::string& name,
-    stk::mesh::Part* part,
-    const int numStates = 0,
-    const int numComponents = 0,
     const void* initVal = nullptr) const
   {
     stk::mesh::PartVector parts;
     parts.push_back(part);
     register_field(name, parts, numStates, numComponents, initVal);
   }
+    
 
   /// Given the named field that has already been registered on the CPU
   /// return the GPU version of the same field.
diff --git a/include/MultiDimViews.h b/include/MultiDimViews.h
index f5353aaa1..c523907f6 100644
--- a/include/MultiDimViews.h
+++ b/include/MultiDimViews.h
@@ -63,7 +63,7 @@ class MultiDimViews
       views_3D_size(0),
       views_4D_size(0)
   {
-    NGP_ThrowRequireMsg(
+    STK_NGP_ThrowRequireMsg(
       (numNeededViews.num1DViews <= maxViewsPerDim ||
        numNeededViews.num2DViews <= maxViewsPerDim ||
        numNeededViews.num3DViews <= maxViewsPerDim ||
@@ -176,7 +176,7 @@ class MultiDimViews
   KOKKOS_FUNCTION
   void add_1D_view(unsigned ordinal, const SharedMemView1D& view)
   {
-    NGP_ThrowRequire(ordinal < maxFieldOrdinals);
+    STK_NGP_ThrowRequire(ordinal < maxFieldOrdinals);
 #if !defined(KOKKOS_ENABLE_GPU)
     views_1D[views_1D_size] = new SharedMemView1D;
     *views_1D[views_1D_size] = view;
@@ -190,7 +190,7 @@ class MultiDimViews
   KOKKOS_FUNCTION
   void add_2D_view(unsigned ordinal, const SharedMemView2D& view)
   {
-    NGP_ThrowRequire(ordinal < maxFieldOrdinals);
+    STK_NGP_ThrowRequire(ordinal < maxFieldOrdinals);
 #if !defined(KOKKOS_ENABLE_GPU)
     views_2D[views_2D_size] = new SharedMemView2D;
     *views_2D[views_2D_size] = view;
@@ -204,7 +204,7 @@ class MultiDimViews
   KOKKOS_FUNCTION
   void add_3D_view(unsigned ordinal, const SharedMemView3D& view)
   {
-    NGP_ThrowRequire(ordinal < maxFieldOrdinals);
+    STK_NGP_ThrowRequire(ordinal < maxFieldOrdinals);
 #if !defined(KOKKOS_ENABLE_GPU)
     views_3D[views_3D_size] = new SharedMemView3D;
     *views_3D[views_3D_size] = view;
@@ -218,7 +218,7 @@ class MultiDimViews
   KOKKOS_FUNCTION
   void add_4D_view(unsigned ordinal, const SharedMemView4D& view)
   {
-    NGP_ThrowRequire(ordinal < maxFieldOrdinals);
+    STK_NGP_ThrowRequire(ordinal < maxFieldOrdinals);
 #if !defined(KOKKOS_ENABLE_GPU)
     views_4D[views_4D_size] = new SharedMemView4D;
     *views_4D[views_4D_size] = view;
diff --git a/include/ScratchViews.h b/include/ScratchViews.h
index 8e863c900..3386303ad 100644
--- a/include/ScratchViews.h
+++ b/include/ScratchViews.h
@@ -61,7 +61,7 @@ count_needed_field_views(const FieldInfoViewType& neededFields)
         numNeededViews.num3DViews++;
       }
     } else {
-      NGP_ThrowRequireMsg(false, "Unknown stk-rank");
+      STK_NGP_ThrowRequireMsg(false, "Unknown stk-rank");
     }
   }
 
@@ -203,7 +203,7 @@ create_needed_field_views(
         numScalars += nodesPerEntity * scalarsDim1 * scalarsDim2;
       }
     } else {
-      NGP_ThrowRequireMsg(false, "Unknown stk-rank");
+      STK_NGP_ThrowRequireMsg(false, "Unknown stk-rank");
     }
   }
 
@@ -257,7 +257,7 @@ class ScratchViews
   MasterElementViews<T, TEAMHANDLETYPE, SHMEM>&
   get_me_views(const COORDS_TYPES cType)
   {
-    NGP_ThrowRequire(hasCoordField[cType] == true);
+    STK_NGP_ThrowRequire(hasCoordField[cType] == true);
     return meViews[cType];
   }
 
@@ -310,7 +310,7 @@ SharedMemView<T*, SHMEM>&
 ScratchViews<T, TEAMHANDLETYPE, SHMEM>::get_scratch_view_1D(
   const stk::mesh::FieldBase& field)
 {
-  //  ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
+  //  STK_ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
   //  "ScratchViews ERROR, trying to get 1D scratch-view for field
   //  "<<field.name()<<" which wasn't declared as pre-req field.");
   //  ViewT<SharedMemView<T*>>* vt =
@@ -323,7 +323,7 @@ SharedMemView<T**, SHMEM>&
 ScratchViews<T, TEAMHANDLETYPE, SHMEM>::get_scratch_view_2D(
   const stk::mesh::FieldBase& field)
 {
-  //  ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
+  //  STK_ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
   //  "ScratchViews ERROR, trying to get 2D scratch-view for field
   //  "<<field.name()<<" which wasn't declared as pre-req field.");
   //  ViewT<SharedMemView<T**>>* vt =
@@ -336,7 +336,7 @@ SharedMemView<T***, SHMEM>&
 ScratchViews<T, TEAMHANDLETYPE, SHMEM>::get_scratch_view_3D(
   const stk::mesh::FieldBase& field)
 {
-  //  ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
+  //  STK_ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
   //  "ScratchViews ERROR, trying to get 3D scratch-view for field
   //  "<<field.name()<<" which wasn't declared as pre-req field.");
   //  ViewT<SharedMemView<T***>>* vt =
@@ -349,7 +349,7 @@ SharedMemView<T****, SHMEM>&
 ScratchViews<T, TEAMHANDLETYPE, SHMEM>::get_scratch_view_4D(
   const stk::mesh::FieldBase& field)
 {
-  //  ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
+  //  STK_ThrowAssertMsg(fieldViews[field.mesh_meta_data_ordinal()] != nullptr,
   //  "ScratchViews ERROR, trying to get 4D scratch-view for field
   //  "<<field.name()<<" which wasn't declared as pre-req field.");
   //  ViewT<SharedMemView<T****>>* vt =
@@ -416,7 +416,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
   for (unsigned i = 0; i < dataEnums.size(); ++i) {
     switch (dataEnums(i)) {
     case FC_AREAV:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFaceIp > 0,
         "ERROR, meFC must be non-null if FC_AREAV is requested.");
       fc_areav =
@@ -424,7 +424,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       numScalars += numFaceIp * nDim;
       break;
     case FC_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFaceIp > 0,
         "ERROR, meFC must be non-null if FC_SHAPE_FCN is requested");
       fc_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -432,7 +432,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       numScalars += numFaceIp * nodesPerFace;
       break;
     case FC_SHIFTED_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFaceIp > 0,
         "ERROR, meFC must be non-null if FC_SHIFTED_SHAPE_FCN is requested");
       fc_shifted_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -440,7 +440,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       numScalars += numFaceIp * nodesPerFace;
       break;
     case SCS_FACE_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFaceIp > 0,
         "ERROR, meSCS must be non-null if SCS_FACE_GRAD_OP is requested.");
       dndx_fc_scs = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -450,7 +450,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       needDetjFC = true;
       break;
     case SCS_SHIFTED_FACE_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFaceIp > 0, "ERROR, meSCS must be non-null if "
                        "SCS_SHIFTED_FACE_GRAD_OP is requested.");
       dndx_shifted_fc_scs = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -460,7 +460,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       needDetjFC = true;
       break;
     case SCS_AREAV:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScsIp > 0,
         "ERROR, meSCS must be non-null if SCS_AREAV is requested.");
       scs_areav =
@@ -469,7 +469,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCS_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScsIp > 0,
         "ERROR, meSCS must be non-null if SCS_GRAD_OP is requested.");
       dndx = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -480,7 +480,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCS_SHIFTED_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScsIp > 0,
         "ERROR, meSCS must be non-null if SCS_SHIFTED_GRAD_OP is requested.");
       dndx_shifted = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -491,7 +491,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCS_GIJ:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScsIp > 0, "ERROR, meSCS must be non-null if SCS_GIJ is requested.");
       gijUpper =
         get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(team, numScsIp, nDim, nDim);
@@ -502,7 +502,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCV_MIJ:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScvIp > 0, "ERROR, meSCV must be non-null if SCV_MIJ is requested.");
       metric =
         get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(team, numScvIp, nDim, nDim);
@@ -511,7 +511,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCS_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScsIp > 0,
         "ERROR, meSCS must be non-null if SCS_SHAPE_FCN is requested");
       scs_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -520,7 +520,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCS_SHIFTED_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScsIp > 0,
         "ERROR, meSCS must be non-null if SCS_SHIFTED_SHAPE_FCN is requested");
       scs_shifted_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -529,7 +529,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCV_VOLUME:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScvIp > 0,
         "ERROR, meSCV must be non-null if SCV_VOLUME is requested.");
       scv_volume = get_shmem_view_1D<T, TEAMHANDLETYPE, SHMEM>(team, numScvIp);
@@ -537,7 +537,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCV_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScvIp > 0,
         "ERROR, meSCV must be non-null if SCV_GRAD_OP is requested.");
       dndx_scv = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -548,7 +548,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCV_SHIFTED_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScvIp > 0,
         "ERROR, meSCV must be non-null if SCV_SHIFTED_GRAD_OP is requested.");
       dndx_scv_shifted = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -559,7 +559,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCV_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScvIp > 0,
         "ERROR, meSCV must be non-null if SCV_SHAPE_FCN is requested");
       scv_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -568,7 +568,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case SCV_SHIFTED_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numScvIp > 0,
         "ERROR, meSCV must be non-null if SCV_SHIFTED_SHAPE_FCN is requested");
       scv_shifted_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -577,7 +577,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case FEM_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFemIp > 0,
         "ERROR, meFEM must be non-null if FEM_GRAD_OP is requested.");
       dndx_fem = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -589,7 +589,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case FEM_SHIFTED_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFemIp > 0,
         "ERROR, meFEM must be non-null if FEM_SHIFTED_GRAD_OP is requested.");
       dndx_fem = get_shmem_view_3D<T, TEAMHANDLETYPE, SHMEM>(
@@ -601,7 +601,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case FEM_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFemIp > 0,
         "ERROR, meFEM must be non-null if FEM_SHAPE_FCN is requested");
       fem_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -610,7 +610,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
       break;
 
     case FEM_SHIFTED_SHAPE_FCN:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         numFemIp > 0,
         "ERROR, meFEM must be non-null if FEM_SHIFTED_SHAPE_FCN is requested");
       fem_shifted_shape_fcn = get_shmem_view_2D<T, TEAMHANDLETYPE, SHMEM>(
@@ -669,7 +669,7 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::create_master_element_views(
 
   // error check
   if (femGradOp && femShiftedGradOp)
-    NGP_ThrowRequireMsg(
+    STK_NGP_ThrowRequireMsg(
       numFemIp > 0, "ERROR, femGradOp and femShiftedGradOp both requested.");
 
   return numScalars;
@@ -737,113 +737,113 @@ MasterElementViews<T, TEAMHANDLETYPE, SHMEM>::fill_master_element_views_new_me(
   for (unsigned i = 0; i < dataEnums.size(); ++i) {
     switch (dataEnums(i)) {
     case FC_AREAV:
-      NGP_ThrowRequireMsg(false, "FC_AREAV not implemented yet.");
+      STK_NGP_ThrowRequireMsg(false, "FC_AREAV not implemented yet.");
       break;
     case SCS_AREAV:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr,
         "ERROR, meSCS needs to be non-null if SCS_AREAV is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCS_AREAV requested.");
       meSCS->determinant(*coordsView, scs_areav);
       break;
     case SCS_FACE_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr,
         "ERROR, meSCS needs to be non-null if SCS_FACE_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr,
         "ERROR, coords null but SCS_FACE_GRAD_OP requested.");
       meSCS->face_grad_op(faceOrdinal, *coordsView, dndx_fc_scs, deriv_fc_scs);
       break;
     case SCS_SHIFTED_FACE_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr, "ERROR, meSCS needs to be non-null if "
                           "SCS_SHIFTED_FACE_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr,
         "ERROR, coords null but SCS_SHIFTED_FACE_GRAD_OP requested.");
       meSCS->shifted_face_grad_op(
         faceOrdinal, *coordsView, dndx_shifted_fc_scs, deriv_fc_scs);
       break;
     case SCS_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr,
         "ERROR, meSCS needs to be non-null if SCS_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCS_GRAD_OP requested.");
       meSCS->grad_op(*coordsView, dndx, deriv);
       break;
     case SCS_SHIFTED_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr,
         "ERROR, meSCS needs to be non-null if SCS_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCS_GRAD_OP requested.");
       meSCS->shifted_grad_op(*coordsView, dndx_shifted, deriv);
       break;
     case SCS_GIJ:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr,
         "ERROR, meSCS needs to be non-null if SCS_GIJ is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCS_GIJ requested.");
       meSCS->gij(*coordsView, gijUpper, gijLower, deriv);
       break;
     case SCS_MIJ:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCS != nullptr,
         "ERROR, meSCV needs to be non-null if SCS_MIJ is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCS_MIJ requested.");
       meSCS->Mij(*coordsView, metric, deriv);
       break;
     case SCV_MIJ:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCV != nullptr,
         "ERROR, meSCV needs to be non-null if SCV_MIJ is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCV_MIJ requested.");
       meSCV->Mij(*coordsView, metric, deriv_scv);
       break;
     case SCV_VOLUME:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCV != nullptr,
         "ERROR, meSCV needs to be non-null if SCV_VOLUME is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCV_VOLUME requested.");
       meSCV->determinant(*coordsView, scv_volume);
       break;
     case SCV_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCV != nullptr,
         "ERROR, meSCV needs to be non-null if SCV_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but SCV_GRAD_OP requested.");
       meSCV->grad_op(*coordsView, dndx_scv, deriv_scv);
       break;
     case SCV_SHIFTED_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meSCV != nullptr, "ERROR, meSCV needs to be non-null if "
                           "SCV_SHIFTED_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr,
         "ERROR, coords null but SCV_SHIFTED_GRAD_OP requested.");
       meSCV->shifted_grad_op(*coordsView, dndx_scv_shifted, deriv_scv);
       break;
     case FEM_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meFEM != nullptr,
         "ERROR, meFEM needs to be non-null if FEM_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but FEM_GRAD_OP requested.");
       meFEM->grad_op_fem(*coordsView, dndx_fem, deriv_fem, det_j_fem);
       break;
     case FEM_SHIFTED_GRAD_OP:
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         meFEM != nullptr, "ERROR, meFEM needs to be non-null if "
                           "FEM_SHIFTED_GRAD_OP is requested.");
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         coordsView != nullptr, "ERROR, coords null but FEM_GRAD_OP requested.");
       meFEM->shifted_grad_op_fem(*coordsView, dndx_fem, deriv_fem, det_j_fem);
       break;
diff --git a/include/SmartField.h b/include/SmartField.h
index 937d9ea62..af4f58005 100644
--- a/include/SmartField.h
+++ b/include/SmartField.h
@@ -110,6 +110,13 @@ class SmartField<FieldType, tags::LEGACY, ACCESS>
     return stk::mesh::field_data(stkField_, entity);
   }
 
+  template <typename A = ACCESS>
+  inline typename std::enable_if_t<!std::is_same_v<A, READ>, T>*
+  operator()(const stk::mesh::Bucket& bucket) const
+  {
+    return stk::mesh::field_data(stkField_, bucket);
+  }
+
   // --- Const Accessors
   template <typename A = ACCESS>
   inline const typename std::enable_if_t<std::is_same_v<A, READ>, T>*
@@ -125,6 +132,13 @@ class SmartField<FieldType, tags::LEGACY, ACCESS>
     return stk::mesh::field_data(stkField_, entity);
   }
 
+  template <typename A = ACCESS>
+  inline const typename std::enable_if_t<std::is_same_v<A, READ>, T>*
+  operator()(const stk::mesh::Bucket& bucket) const
+  {
+    return stk::mesh::field_data(stkField_, bucket);
+  }
+
   ~SmartField()
   {
     if (is_write_) {
diff --git a/include/TpetraLinearSystem.h b/include/TpetraLinearSystem.h
index b676df031..554013bbc 100644
--- a/include/TpetraLinearSystem.h
+++ b/include/TpetraLinearSystem.h
@@ -363,7 +363,7 @@ copy_kokkos_unordered_map(
       fail_count += insert_result.failed() ? 1 : 0;
     }
   }
-  ThrowRequire(fail_count == 0);
+  STK_ThrowRequire(fail_count == 0);
 }
 
 int getDofStatus_impl(stk::mesh::Entity node, const Realm& realm);
diff --git a/include/aero/actuator/ActuatorGenericSearchFunctor.h b/include/aero/actuator/ActuatorGenericSearchFunctor.h
index c1941d606..408a3fa5a 100644
--- a/include/aero/actuator/ActuatorGenericSearchFunctor.h
+++ b/include/aero/actuator/ActuatorGenericSearchFunctor.h
@@ -83,8 +83,8 @@ struct GenericLoopOverCoarseSearchResults
     const int numIp = meSCV->num_integration_points();
 
     // just allocate for largest expected size (hex27)
-    ThrowAssert(numIp <= 216);
-    ThrowAssert(numNodes <= 27);
+    STK_ThrowAssert(numIp <= 216);
+    STK_ThrowAssert(numNodes <= 27);
 
     double scvip[216];
     double elemcoords[27 * 3];
diff --git a/include/aero/actuator/ActuatorScalingFLLC.h b/include/aero/actuator/ActuatorScalingFLLC.h
index ddb80d5ac..62f8e2fc0 100644
--- a/include/aero/actuator/ActuatorScalingFLLC.h
+++ b/include/aero/actuator/ActuatorScalingFLLC.h
@@ -57,7 +57,7 @@ scale_lift_force(
   case (ActuatorType::ActLineFASTNGP):
   case (ActuatorType::ActDiskFASTNGP): {
 #ifndef NALU_USES_OPENFAST
-    ThrowErrorMsg("Actuator methods require OpenFAST");
+    STK_ThrowErrorMsg("Actuator methods require OpenFAST");
 #if !defined(KOKKOS_ENABLE_GPU)
     break;
 #endif
diff --git a/include/aero/actuator/UtilitiesActuator.h b/include/aero/actuator/UtilitiesActuator.h
index 1f9ef605b..58d12598f 100644
--- a/include/aero/actuator/UtilitiesActuator.h
+++ b/include/aero/actuator/UtilitiesActuator.h
@@ -87,8 +87,8 @@ template <typename T>
 inline void
 reduce_view_on_host(T view)
 {
-  ThrowAssert(view.size() > 0);
-  ThrowAssert(view.data());
+  STK_ThrowAssert(view.size() > 0);
+  STK_ThrowAssert(view.data());
   MPI_Datatype mpi_type;
   if (std::is_same<typename T::value_type, double>::value) {
     mpi_type = MPI_DOUBLE;
@@ -99,7 +99,7 @@ reduce_view_on_host(T view)
   } else if (std::is_same<typename T::value_type, uint64_t>::value) {
     mpi_type = MPI_LONG;
   } else {
-    ThrowErrorMsg("unsupported type to reduce view on host");
+    STK_ThrowErrorMsg("unsupported type to reduce view on host");
   }
 
   MPI_Allreduce(
diff --git a/include/kernel/KernelBuilder.h b/include/kernel/KernelBuilder.h
index 100770073..deaa4fda3 100644
--- a/include/kernel/KernelBuilder.h
+++ b/include/kernel/KernelBuilder.h
@@ -48,7 +48,7 @@ build_or_add_part_to_solver_alg(
       topo == stk::topology::QUADRILATERAL_4_2D ||
       topo == stk::topology::TRIANGLE_3_2D || topo == stk::topology::WEDGE_6 ||
       topo == stk::topology::TETRAHEDRON_4 || topo == stk::topology::PYRAMID_5);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     !isNotNGP, "Consolidated algorithm called on non-NGP MasterElement");
 
   auto itc = solverAlgs.find(algName);
@@ -57,7 +57,7 @@ build_or_add_part_to_solver_alg(
     auto* theSolverAlg = new AssembleElemSolverAlgorithm(
       eqSys.realm_, &part, &eqSys, stk::topology::ELEMENT_RANK,
       topo.num_nodes());
-    ThrowRequire(theSolverAlg != nullptr);
+    STK_ThrowRequire(theSolverAlg != nullptr);
 
     NaluEnv::self().naluOutputP0()
       << "Created the following interior elem alg: " << algName << std::endl;
@@ -71,7 +71,7 @@ build_or_add_part_to_solver_alg(
 
   auto* theSolverAlg =
     dynamic_cast<AssembleElemSolverAlgorithm*>(solverAlgs.at(algName));
-  ThrowRequire(theSolverAlg != nullptr);
+  STK_ThrowRequire(theSolverAlg != nullptr);
 
   return {theSolverAlg, createNewAlg};
 }
@@ -80,7 +80,7 @@ template <template <typename> class T, typename... Args>
 Kernel*
 build_fem_kernel(stk::topology topo, Args&&... args)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     topo == stk::topology::HEXAHEDRON_8,
     "FEM kernels only implemented for Hex8 topology");
   return new T<AlgTraitsHex8>(std::forward<Args>(args)...);
@@ -140,7 +140,7 @@ class KernelBuilder
       if (eqSys_.supp_alg_is_requested(name)) {
         Kernel* compKernel =
           build_topo_kernel<T>(part_.topology(), std::forward<Args>(args)...);
-        ThrowRequire(compKernel != nullptr);
+        STK_ThrowRequire(compKernel != nullptr);
         KernelBuilderLog::self().add_built_name(eqSys_.eqnTypeName_, name);
         solverAlg_->activeKernels_.push_back(compKernel);
         isCreated = true;
@@ -159,7 +159,7 @@ class KernelBuilder
       if (eqSys_.supp_alg_is_requested(name)) {
         Kernel* compKernel =
           build_fem_kernel<T>(part_.topology(), std::forward<Args>(args)...);
-        ThrowRequire(compKernel != nullptr);
+        STK_ThrowRequire(compKernel != nullptr);
         KernelBuilderLog::self().add_built_name(eqSys_.eqnTypeName_, name);
         solverAlg_->activeKernels_.push_back(compKernel);
         isCreated = true;
@@ -191,7 +191,7 @@ build_face_elem_topo_kernel(
     case stk::topology::WEDGE_6:
       return new T<AlgTraitsQuad4Wed6>(std::forward<Args>(args)...);
     default:
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         false, "Quad4 exposed face is not attached to either a hex8, pyr5, or "
                "wedge6.");
       return nullptr;
@@ -205,7 +205,7 @@ build_face_elem_topo_kernel(
     case stk::topology::WEDGE_6:
       return new T<AlgTraitsTri3Wed6>(std::forward<Args>(args)...);
     default:
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         false,
         "Tri3 exposed face is not attached to either a tet4, pyr5, or wedge6.");
       return nullptr;
@@ -252,7 +252,7 @@ build_topo_kernel_if_requested(
   if (eqSys.supp_alg_is_requested(name)) {
     Kernel* compKernel =
       build_topo_kernel<T>(topo, std::forward<Args>(args)...);
-    ThrowRequire(compKernel != nullptr);
+    STK_ThrowRequire(compKernel != nullptr);
     KernelBuilderLog::self().add_built_name(eqSys.eqnTypeName_, name);
     kernelVec.push_back(compKernel);
     isCreated = true;
@@ -273,7 +273,7 @@ build_face_elem_topo_kernel_automatic(
   KernelBuilderLog::self().add_valid_name(eqSys.eqnTypeName_, name);
   Kernel* compKernel = build_face_elem_topo_kernel<T>(
     faceTopo, elemTopo, std::forward<Args>(args)...);
-  ThrowRequire(compKernel != nullptr);
+  STK_ThrowRequire(compKernel != nullptr);
   KernelBuilderLog::self().add_built_name(eqSys.eqnTypeName_, name);
   kernelVec.push_back(compKernel);
   return true;
@@ -292,7 +292,7 @@ build_fem_kernel_if_requested(
   KernelBuilderLog::self().add_valid_name(eqSys.eqnTypeName_, name);
   if (eqSys.supp_alg_is_requested(name)) {
     Kernel* compKernel = build_fem_kernel<T>(topo, std::forward<Args>(args)...);
-    ThrowRequire(compKernel != nullptr);
+    STK_ThrowRequire(compKernel != nullptr);
     KernelBuilderLog::self().add_built_name(eqSys.eqnTypeName_, name);
     kernelVec.push_back(compKernel);
     isCreated = true;
@@ -312,7 +312,7 @@ build_face_topo_kernel_automatic(
   KernelBuilderLog::self().add_valid_name(eqSys.eqnTypeName_, name);
   Kernel* compKernel =
     build_face_topo_kernel<T>(topo, std::forward<Args>(args)...);
-  ThrowRequire(compKernel != nullptr);
+  STK_ThrowRequire(compKernel != nullptr);
   KernelBuilderLog::self().add_built_name(eqSys.eqnTypeName_, name);
   kernelVec.push_back(compKernel);
   return true;
@@ -339,7 +339,7 @@ build_or_add_part_to_face_elem_solver_alg(
       elemTopo == stk::topology::WEDGE_6 ||
       elemTopo == stk::topology::TETRAHEDRON_4 ||
       elemTopo == stk::topology::PYRAMID_5);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     !isNotNGP, "Consolidated algorithm called on non-NGP MasterElement");
 
   auto itc = solverAlgs.find(algName);
@@ -347,7 +347,7 @@ build_or_add_part_to_face_elem_solver_alg(
   if (createNewAlg) {
     auto* theSolverAlg = new AssembleFaceElemSolverAlgorithm(
       eqSys.realm_, &part, &eqSys, topo.num_nodes(), elemTopo.num_nodes());
-    ThrowRequire(theSolverAlg != nullptr);
+    STK_ThrowRequire(theSolverAlg != nullptr);
 
     NaluEnv::self().naluOutputP0()
       << "Created the following bc face/elem alg: " << algName << std::endl;
@@ -361,7 +361,7 @@ build_or_add_part_to_face_elem_solver_alg(
 
   auto* theSolverAlg =
     dynamic_cast<AssembleFaceElemSolverAlgorithm*>(solverAlgs.at(algName));
-  ThrowRequire(theSolverAlg != nullptr);
+  STK_ThrowRequire(theSolverAlg != nullptr);
 
   return {theSolverAlg, createNewAlg};
 }
@@ -381,7 +381,7 @@ build_or_add_part_to_face_bc_solver_alg(
     !(topo == stk::topology::QUAD_4 || topo == stk::topology::QUAD_9 ||
       topo == stk::topology::TRI_3 || topo == stk::topology::LINE_2 ||
       topo == stk::topology::LINE_3);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     !isNotNGP, "Consolidated algorithm called on non-NGP MasterElement");
 
   auto itc = solverAlgs.find(algName);
@@ -390,7 +390,7 @@ build_or_add_part_to_face_bc_solver_alg(
     auto* theSolverAlg = new AssembleElemSolverAlgorithm(
       eqSys.realm_, &part, &eqSys, eqSys.realm_.meta_data().side_rank(),
       topo.num_nodes());
-    ThrowRequire(theSolverAlg != nullptr);
+    STK_ThrowRequire(theSolverAlg != nullptr);
 
     NaluEnv::self().naluOutputP0()
       << "Created the following bc face alg: " << algName << std::endl;
@@ -404,7 +404,7 @@ build_or_add_part_to_face_bc_solver_alg(
 
   auto* theSolverAlg =
     dynamic_cast<AssembleElemSolverAlgorithm*>(solverAlgs.at(algName));
-  ThrowRequire(theSolverAlg != nullptr);
+  STK_ThrowRequire(theSolverAlg != nullptr);
 
   return {theSolverAlg, createNewAlg};
 }
diff --git a/include/master_element/Hex8FEM.h b/include/master_element/Hex8FEM.h
index 3eac756d7..e7c0b5831 100644
--- a/include/master_element/Hex8FEM.h
+++ b/include/master_element/Hex8FEM.h
@@ -84,7 +84,7 @@ class Hex8FEM : public MasterElement
                                     +glIP, +glIP, +glIP, -glIP, +glIP, +glIP};
 
   // standard integration location +/ sqrt(3)/3
-  const double gIP = std::sqrt(3.0) / 3.0;
+  const double gIP = 0.577350269189626; // std::sqrt(3.0) / 3.0;
   const double intgLoc_[numIntPoints_ * nDim_] = {
     -gIP, -gIP, -gIP, +gIP, -gIP, -gIP, +gIP, +gIP, -gIP, -gIP, +gIP, -gIP,
     -gIP, -gIP, +gIP, +gIP, -gIP, +gIP, +gIP, +gIP, +gIP, -gIP, +gIP, +gIP};
diff --git a/include/master_element/MasterElement.h b/include/master_element/MasterElement.h
index d8c5e32b3..21770a13d 100644
--- a/include/master_element/MasterElement.h
+++ b/include/master_element/MasterElement.h
@@ -58,7 +58,7 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* gradop */,
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */)
   {
-    NGP_ThrowErrorMsg("MasterElement::grad_op not implemented for element");
+    STK_NGP_ThrowErrorMsg("MasterElement::grad_op not implemented for element");
   }
 
   virtual void grad_op(
@@ -66,7 +66,7 @@ class MasterElement
     SharedMemView<double***>& /* gradop */,
     SharedMemView<double***>& /* deriv */)
   {
-    NGP_ThrowErrorMsg("MasterElement::grad_op not implemented for element");
+    STK_NGP_ThrowErrorMsg("MasterElement::grad_op not implemented for element");
   }
 
   KOKKOS_FUNCTION virtual void shifted_grad_op(
@@ -74,7 +74,7 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* gradop */,
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::shifted_grad_op not implemented for element");
   }
 
@@ -84,7 +84,7 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* gradop */,
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::face_grad_op not implemented for element");
   }
 
@@ -94,7 +94,7 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* gradop */,
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::shifted_face_grad_op not implemented for element");
   }
 
@@ -104,7 +104,8 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */,
     SharedMemView<DoubleType*, DeviceShmem>& /*det_j*/)
   {
-    NGP_ThrowErrorMsg("MasterElement::grad_op_fem not implemented for element");
+    STK_NGP_ThrowErrorMsg(
+      "MasterElement::grad_op_fem not implemented for element");
   }
 
   KOKKOS_FUNCTION virtual void shifted_grad_op_fem(
@@ -113,7 +114,7 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */,
     SharedMemView<DoubleType*, DeviceShmem>& /*det_j*/)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::shifted_grad_op_fem not implemented for element");
   }
 
@@ -121,15 +122,16 @@ class MasterElement
     const SharedMemView<DoubleType**, DeviceShmem>& /* coords */,
     SharedMemView<DoubleType**, DeviceShmem>& /* areav */)
   {
-    NGP_ThrowErrorMsg("MasterElement::determinant not implemented for element: "
-                      "DoubleType area");
+    STK_NGP_ThrowErrorMsg(
+      "MasterElement::determinant not implemented for element: "
+      "DoubleType area");
   }
 
   virtual void determinant(
     const SharedMemView<double**>& /* coords */,
     SharedMemView<double**>& /* areav */)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::determinant not implemented for element: double area");
   }
 
@@ -139,7 +141,7 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* glower */,
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */)
   {
-    NGP_ThrowErrorMsg("MasterElement::gij not implemented for element");
+    STK_NGP_ThrowErrorMsg("MasterElement::gij not implemented for element");
   }
 
   KOKKOS_FUNCTION virtual void Mij(
@@ -147,22 +149,23 @@ class MasterElement
     SharedMemView<DoubleType***, DeviceShmem>& /* metric */,
     SharedMemView<DoubleType***, DeviceShmem>& /* deriv */)
   {
-    NGP_ThrowErrorMsg("MasterElement::Mij not implemented for element");
+    STK_NGP_ThrowErrorMsg("MasterElement::Mij not implemented for element");
   }
 
   KOKKOS_FUNCTION virtual void determinant(
     const SharedMemView<DoubleType**, DeviceShmem>& /* coords */,
     SharedMemView<DoubleType*, DeviceShmem>& /* volume */)
   {
-    NGP_ThrowErrorMsg("MasterElement::determinant not implemented for element: "
-                      "DoubleType volume");
+    STK_NGP_ThrowErrorMsg(
+      "MasterElement::determinant not implemented for element: "
+      "DoubleType volume");
   }
 
   virtual void determinant(
     const SharedMemView<double**>& /* coords */,
     SharedMemView<double*>& /* volume */)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::determinant not implemented for element: double volume");
   }
 
@@ -180,33 +183,33 @@ class MasterElement
 
   KOKKOS_FUNCTION virtual const int* adjacentNodes()
   {
-    NGP_ThrowErrorMsg("MasterElement::adjacentNodes not implemented");
+    STK_NGP_ThrowErrorMsg("MasterElement::adjacentNodes not implemented");
     return nullptr;
   }
 
   KOKKOS_FUNCTION virtual const int* scsIpEdgeOrd()
   {
-    NGP_ThrowErrorMsg("MasterElement::scsIpEdgeOrd not implemented");
+    STK_NGP_ThrowErrorMsg("MasterElement::scsIpEdgeOrd not implemented");
     return nullptr;
   }
 
   KOKKOS_FUNCTION virtual const int* ipNodeMap(int /* ordinal */ = 0) const
   {
-    NGP_ThrowErrorMsg("MasterElement::ipNodeMap not implemented");
+    STK_NGP_ThrowErrorMsg("MasterElement::ipNodeMap not implemented");
     return nullptr;
   }
 
   KOKKOS_FUNCTION virtual int
   opposingNodes(const int /* ordinal */, const int /* node */)
   {
-    NGP_ThrowErrorMsg("opposingNodes not implemented");
+    STK_NGP_ThrowErrorMsg("opposingNodes not implemented");
     return -1;
   }
 
   KOKKOS_FUNCTION virtual int
   opposingFace(const int /* ordinal */, const int /* node */)
   {
-    NGP_ThrowErrorMsg("opposingFace not implemented");
+    STK_NGP_ThrowErrorMsg("opposingFace not implemented");
     return -1;
   }
 
@@ -293,7 +296,7 @@ class MasterElement
   }
   KOKKOS_FUNCTION virtual const int* side_node_ordinals(int) const
   {
-    NGP_ThrowErrorMsg("side_node_ordinals not implemented");
+    STK_NGP_ThrowErrorMsg("side_node_ordinals not implemented");
     return nullptr;
   }
 
@@ -306,22 +309,24 @@ class MasterElement
   KOKKOS_FUNCTION virtual void
   shape_fcn(SharedMemView<DoubleType**, DeviceShmem>&)
   {
-    NGP_ThrowErrorMsg("MasterElement::shape_fcn not implemented for element");
+    STK_NGP_ThrowErrorMsg(
+      "MasterElement::shape_fcn not implemented for element");
   }
   virtual void shape_fcn(SharedMemView<double**, HostShmem>&)
   {
-    NGP_ThrowErrorMsg("MasterElement::shape_fcn not implemented for element");
+    STK_NGP_ThrowErrorMsg(
+      "MasterElement::shape_fcn not implemented for element");
   }
 
   KOKKOS_FUNCTION virtual void
   shifted_shape_fcn(SharedMemView<DoubleType**, DeviceShmem>&)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::shifted_shape_fcn not implemented for element");
   }
   virtual void shifted_shape_fcn(SharedMemView<double**, HostShmem>&)
   {
-    NGP_ThrowErrorMsg(
+    STK_NGP_ThrowErrorMsg(
       "MasterElement::shifted_shape_fcn not implemented for element");
   }
 
diff --git a/include/master_element/MasterElementFunctions.h b/include/master_element/MasterElementFunctions.h
index cdb0c2bc9..66840e784 100644
--- a/include/master_element/MasterElementFunctions.h
+++ b/include/master_element/MasterElementFunctions.h
@@ -82,10 +82,11 @@ generic_grad_op(
     CoordViewType::rank == 2, "Coordinate view assumed to be rank 2");
   static_assert(OutputViewType::rank == 3, "Weight view assumed to be rank 3");
 
-  ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
-  ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
+  STK_ThrowAssert(
+    AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
+  STK_ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
   for (int i = 0; i < dim; ++i)
-    ThrowAssert(weights.extent(i) == referenceGradWeights.extent(i));
+    STK_ThrowAssert(weights.extent(i) == referenceGradWeights.extent(i));
 
   for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
     NALU_ALIGNED ftype jact[dim][dim];
@@ -111,7 +112,7 @@ generic_grad_op(
     NALU_ALIGNED ftype det = ftype(0.0);
     for (int i = 0; i < dim; ++i)
       det += jact[i][0] * adjJac[i][0];
-    ThrowAssertMsg(
+    STK_ThrowAssertMsg(
       stk::simd::are_any(det > tiny_positive_value()),
       "Problem with Jacobian determinant");
 
@@ -598,10 +599,11 @@ generic_determinant_3d(
   static_assert(OutputViewType::rank == 1, "Weight view assumed to be 1D");
   static_assert(AlgTraits::nDim_ == 3, "3D method");
 
-  ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
-  ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
+  STK_ThrowAssert(
+    AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
+  STK_ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
 
-  ThrowAssert(detj.extent(0) == referenceGradWeights.extent(0));
+  STK_ThrowAssert(detj.extent(0) == referenceGradWeights.extent(0));
 
   for (unsigned ip = 0; ip < referenceGradWeights.extent(0); ++ip) {
     NALU_ALIGNED ftype jac[3][3] = {
diff --git a/include/matrix_free/LinSysInfo.h b/include/matrix_free/LinSysInfo.h
index a7f23a17e..beaae7403 100644
--- a/include/matrix_free/LinSysInfo.h
+++ b/include/matrix_free/LinSysInfo.h
@@ -35,7 +35,7 @@ struct linsys_info
   static stk::mesh::NgpField<gid_type>
   get_gid_field(const stk::mesh::MetaData& meta)
   {
-    ThrowRequire(meta.get_field(stk::topology::NODE_RANK, gid_name));
+    STK_ThrowRequire(meta.get_field(stk::topology::NODE_RANK, gid_name));
     return stk::mesh::get_updated_ngp_field<gid_type>(
       *meta.get_field(stk::topology::NODE_RANK, gid_name));
   }
diff --git a/include/matrix_free/MomentumOperator.h b/include/matrix_free/MomentumOperator.h
index 7e2a12a18..84d6c8bea 100644
--- a/include/matrix_free/MomentumOperator.h
+++ b/include/matrix_free/MomentumOperator.h
@@ -57,8 +57,9 @@ class MomentumResidualOperator
     dirichlet_bc_active_ = dirichlet_offsets_in.extent_int(0) > 0;
     dirichlet_bc_offsets_ = dirichlet_offsets_in;
     bc_ = bc;
-    ThrowRequire(bc.up1.extent_int(0) == dirichlet_bc_offsets_.extent_int(0));
-    ThrowRequire(bc.ubc.extent_int(0) == bc.up1.extent_int(0));
+    STK_ThrowRequire(
+      bc.up1.extent_int(0) == dirichlet_bc_offsets_.extent_int(0));
+    STK_ThrowRequire(bc.ubc.extent_int(0) == bc.up1.extent_int(0));
   }
 
 private:
diff --git a/include/matrix_free/PolynomialOrders.h b/include/matrix_free/PolynomialOrders.h
index 140f8a6e8..9abc85e51 100644
--- a/include/matrix_free/PolynomialOrders.h
+++ b/include/matrix_free/PolynomialOrders.h
@@ -56,27 +56,27 @@ enum {
 #define P_INVOKEABLE(func)                                                     \
   template <int p, typename... Args>                                           \
   auto func(Args&&... args)                                                    \
-    ->decltype(IMPLNAME(func) < p > ::invoke(std::forward<Args>(args)...))     \
+    -> decltype(IMPLNAME(func) < p > ::invoke(std::forward<Args>(args)...))    \
   {                                                                            \
     return IMPLNAME(func)<p>::invoke(std::forward<Args>(args)...);             \
   }
 
 // can't return a value dependent on template parameter
-#define SWITCH_INVOKEABLE(func)                                                   \
-  template <typename... Args>                                                     \
-  auto func(int p, Args&&... args)                                                \
-    ->decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward<Args>(args)...)) \
-  {                                                                               \
-    switch (p) {                                                                  \
-    case inst::P2:                                                                \
-      return IMPLNAME(func)<inst::P2>::invoke(std::forward<Args>(args)...);       \
-    case inst::P3:                                                                \
-      return IMPLNAME(func)<inst::P3>::invoke(std::forward<Args>(args)...);       \
-    case inst::P4:                                                                \
-      return IMPLNAME(func)<inst::P4>::invoke(std::forward<Args>(args)...);       \
-    default:                                                                      \
-      return IMPLNAME(func)<inst::P1>::invoke(std::forward<Args>(args)...);       \
-    }                                                                             \
+#define SWITCH_INVOKEABLE(func)                                                    \
+  template <typename... Args>                                                      \
+  auto func(int p, Args&&... args)                                                 \
+    -> decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward<Args>(args)...)) \
+  {                                                                                \
+    switch (p) {                                                                   \
+    case inst::P2:                                                                 \
+      return IMPLNAME(func)<inst::P2>::invoke(std::forward<Args>(args)...);        \
+    case inst::P3:                                                                 \
+      return IMPLNAME(func)<inst::P3>::invoke(std::forward<Args>(args)...);        \
+    case inst::P4:                                                                 \
+      return IMPLNAME(func)<inst::P4>::invoke(std::forward<Args>(args)...);        \
+    default:                                                                       \
+      return IMPLNAME(func)<inst::P1>::invoke(std::forward<Args>(args)...);        \
+    }                                                                              \
   }
 
 } // namespace matrix_free
diff --git a/include/ngp_utils/NgpFieldManager.h b/include/ngp_utils/NgpFieldManager.h
index abdfa1095..1b449f7f4 100644
--- a/include/ngp_utils/NgpFieldManager.h
+++ b/include/ngp_utils/NgpFieldManager.h
@@ -47,7 +47,7 @@ class FieldManager
   template <typename T>
   stk::mesh::NgpField<T>& get_field(unsigned fieldOrdinal) const
   {
-    ThrowAssertMsg(
+    STK_ThrowAssertMsg(
       m_meta.get_fields().size() > fieldOrdinal, "Invalid field ordinal.");
     stk::mesh::FieldBase* stkField = m_meta.get_fields()[fieldOrdinal];
     stk::mesh::NgpField<T>& tmp =
diff --git a/include/ngp_utils/NgpFieldOps.h b/include/ngp_utils/NgpFieldOps.h
index 5225877f6..e35121972 100644
--- a/include/ngp_utils/NgpFieldOps.h
+++ b/include/ngp_utils/NgpFieldOps.h
@@ -361,7 +361,7 @@ template <typename Mesh, typename Field>
 KOKKOS_INLINE_FUNCTION impl::NodeFieldOp<Mesh, Field, ElemSimdData<Mesh>>
 simd_elem_nodal_field_updater(const Mesh& mesh, const Field& fld)
 {
-  NGP_ThrowAssert(fld.get_rank() == stk::topology::NODE_RANK);
+  STK_NGP_ThrowAssert(fld.get_rank() == stk::topology::NODE_RANK);
   return impl::NodeFieldOp<Mesh, Field, ElemSimdData<Mesh>>{mesh, fld};
 }
 
@@ -369,7 +369,7 @@ template <typename Mesh, typename Field>
 KOKKOS_INLINE_FUNCTION impl::NodeFieldOp<Mesh, Field, FaceElemSimdData<Mesh>>
 simd_face_elem_nodal_field_updater(const Mesh& mesh, const Field& fld)
 {
-  NGP_ThrowAssert(fld.get_rank() == stk::topology::NODE_RANK);
+  STK_NGP_ThrowAssert(fld.get_rank() == stk::topology::NODE_RANK);
   return impl::NodeFieldOp<Mesh, Field, FaceElemSimdData<Mesh>>{mesh, fld};
 }
 
@@ -379,7 +379,7 @@ template <typename Mesh, typename Field>
 KOKKOS_INLINE_FUNCTION impl::ElemFieldOp<Mesh, Field, ElemSimdData<Mesh>>
 simd_elem_field_updater(const Mesh&, const Field& fld)
 {
-  NGP_ThrowAssert(
+  STK_NGP_ThrowAssert(
     (fld.get_rank() == stk::topology::ELEM_RANK) ||
     (fld.get_rank() == stk::topology::FACE_RANK) ||
     (fld.get_rank() == stk::topology::EDGE_RANK));
@@ -392,7 +392,7 @@ template <typename Mesh, typename Field>
 KOKKOS_INLINE_FUNCTION impl::ElemFieldOp<Mesh, Field, FaceElemSimdData<Mesh>>
 simd_face_elem_field_updater(const Mesh&, const Field& fld)
 {
-  NGP_ThrowAssert(
+  STK_NGP_ThrowAssert(
     (fld.get_rank() == stk::topology::FACE_RANK) ||
     (fld.get_rank() == stk::topology::EDGE_RANK));
   return impl::ElemFieldOp<Mesh, Field, FaceElemSimdData<Mesh>>{fld};
@@ -402,7 +402,7 @@ template <typename Mesh, typename Field>
 KOKKOS_INLINE_FUNCTION impl::SimpleNodeFieldOp<Mesh, Field>
 edge_nodal_field_updater(const Mesh& mesh, const Field& fld)
 {
-  NGP_ThrowAssert(fld.get_rank() == stk::topology::NODE_RANK);
+  STK_NGP_ThrowAssert(fld.get_rank() == stk::topology::NODE_RANK);
   return impl::SimpleNodeFieldOp<Mesh, Field>{mesh, fld};
 }
 
diff --git a/include/ngp_utils/NgpLoopUtils.h b/include/ngp_utils/NgpLoopUtils.h
index 9d862f176..f537db38e 100644
--- a/include/ngp_utils/NgpLoopUtils.h
+++ b/include/ngp_utils/NgpLoopUtils.h
@@ -340,7 +340,7 @@ run_elem_algorithm(
   ElemDataRequestsGPU dataReqNGP(fieldMgr, dataReqs);
 
   const int nodesPerElement = nodes_per_entity(dataReqNGP);
-  NGP_ThrowRequire(nodesPerElement != 0);
+  STK_NGP_ThrowRequire(nodesPerElement != 0);
 
   const auto reqType =
     (rank == stk::topology::ELEM_RANK) ? ElemReqType::ELEM : ElemReqType::FACE;
@@ -436,7 +436,7 @@ run_elem_par_reduce(
   ElemDataRequestsGPU dataReqNGP(fieldMgr, dataReqs);
 
   const int nodesPerElement = nodes_per_entity(dataReqNGP);
-  NGP_ThrowRequire(nodesPerElement != 0);
+  STK_NGP_ThrowRequire(nodesPerElement != 0);
 
   const auto reqType =
     (rank == stk::topology::ELEM_RANK) ? ElemReqType::ELEM : ElemReqType::FACE;
@@ -524,8 +524,8 @@ run_face_elem_algorithm(
 
   const int nodesPerElement = nodes_per_entity(elemDataNGP);
   const int nodesPerFace = nodes_per_entity(faceDataNGP, METype::FACE);
-  NGP_ThrowRequire(nodesPerElement != 0);
-  NGP_ThrowRequire(nodesPerFace != 0);
+  STK_NGP_ThrowRequire(nodesPerElement != 0);
+  STK_NGP_ThrowRequire(nodesPerFace != 0);
 
   const int bytes_per_team = 0;
   const int bytes_per_thread =
@@ -648,8 +648,8 @@ run_face_elem_par_reduce(
 
   const int nodesPerElement = nodes_per_entity(elemDataNGP);
   const int nodesPerFace = nodes_per_entity(faceDataNGP, METype::FACE);
-  NGP_ThrowRequire(nodesPerElement != 0);
-  NGP_ThrowRequire(nodesPerFace != 0);
+  STK_NGP_ThrowRequire(nodesPerElement != 0);
+  STK_NGP_ThrowRequire(nodesPerFace != 0);
 
   const int bytes_per_team = 0;
   const int bytes_per_thread =
@@ -776,8 +776,8 @@ run_face_elem_algorithm_nosimd(
 
   const int nodesPerElement = nodes_per_entity(elemDataNGP);
   const int nodesPerFace = nodes_per_entity(faceDataNGP, METype::FACE);
-  NGP_ThrowRequire(nodesPerElement != 0);
-  NGP_ThrowRequire(nodesPerFace != 0);
+  STK_NGP_ThrowRequire(nodesPerElement != 0);
+  STK_NGP_ThrowRequire(nodesPerFace != 0);
 
   const int bytes_per_team = 0;
   const int bytes_per_thread =
@@ -807,7 +807,7 @@ run_face_elem_algorithm_nosimd(
           const auto faceIdx = ngpMesh.fast_mesh_index(face);
           const auto elements = ngpMesh.get_elements(sideRank, faceIdx);
 
-          NGP_ThrowAssert(elements.size() == 1);
+          STK_NGP_ThrowAssert(elements.size() == 1);
           const auto faceOrd =
             ngpMesh.get_element_ordinals(sideRank, faceIdx)[0];
           const auto elem = elements[0];
diff --git a/include/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.h b/include/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.h
index 3fe5910cb..1b1f6791a 100644
--- a/include/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.h
+++ b/include/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.h
@@ -30,21 +30,11 @@ class ThermalConductivityFromPrandtlPropAlgorithm : public Algorithm
 {
 public:
   ThermalConductivityFromPrandtlPropAlgorithm(
-    Realm& realm,
-    const stk::mesh::PartVector& part_vec,
-    ScalarFieldType* thermalCond,
-    ScalarFieldType* specificHeat,
-    ScalarFieldType* viscosity,
-    const double Pr);
+    Realm& realm, const stk::mesh::PartVector& part_vec, const double Pr);
 
   virtual ~ThermalConductivityFromPrandtlPropAlgorithm() {}
 
   virtual void execute();
-
-  ScalarFieldType* thermalCond_;
-  ScalarFieldType* specHeat_;
-  ScalarFieldType* viscosity_;
-
   const double Pr_;
 };
 
diff --git a/include/tabular_props/LU.h b/include/tabular_props/LU.h
index 3188503c9..12e6ae6a7 100644
--- a/include/tabular_props/LU.h
+++ b/include/tabular_props/LU.h
@@ -22,8 +22,8 @@ class LU
   // Writable element access
   inline double& operator()(int row, int col)
   {
-    ThrowRequire(row <= dim_);
-    ThrowRequire(col <= dim_);
+    STK_ThrowRequire(row <= dim_);
+    STK_ThrowRequire(col <= dim_);
     isReady_ = false;
     return AA_(row, col);
   };
@@ -31,8 +31,8 @@ class LU
   // Read-only element access
   inline double operator()(int row, int col) const
   {
-    ThrowRequire(row <= dim_);
-    ThrowRequire(col <= dim_);
+    STK_ThrowRequire(row <= dim_);
+    STK_ThrowRequire(col <= dim_);
     return AA_(row, col);
   };
 
@@ -40,8 +40,8 @@ class LU
   // the overloaded const operator() when possible for some reason.
   inline double value(int row, int col) const
   {
-    ThrowRequire(row <= dim_);
-    ThrowRequire(col <= dim_);
+    STK_ThrowRequire(row <= dim_);
+    STK_ThrowRequire(col <= dim_);
     return AA_(row, col);
   };
 
diff --git a/include/utils/AMSUtils.h b/include/utils/AMSUtils.h
index 6c4a43926..603d8a66d 100644
--- a/include/utils/AMSUtils.h
+++ b/include/utils/AMSUtils.h
@@ -21,7 +21,8 @@ template <class T, int dim = 3>
 KOKKOS_FUNCTION T
 get_M43_constant(T D[dim][dim], const double CMdeg)
 {
-  NGP_ThrowRequireMsg(dim == 3, "Compute of M43 constant requires 3D problem");
+  STK_NGP_ThrowRequireMsg(
+    dim == 3, "Compute of M43 constant requires 3D problem");
 
   // Coefficients for the polynomial
   double c[15] = {0.971903113666644,  0.065591700544879,  0.071103489538998,
diff --git a/include/utils/StkHelpers.h b/include/utils/StkHelpers.h
index 394c20f7a..600e2d20f 100644
--- a/include/utils/StkHelpers.h
+++ b/include/utils/StkHelpers.h
@@ -60,7 +60,7 @@ get_field_ordinal(
   const stk::mesh::EntityRank entity_rank = stk::topology::NODE_RANK)
 {
   stk::mesh::FieldBase* field = meta.get_field(entity_rank, fieldName);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     (field != nullptr), "Requested field does not exist: " + fieldName);
   return field->mesh_meta_data_ordinal();
 }
@@ -76,9 +76,9 @@ get_field_ordinal(
   const stk::mesh::EntityRank entity_rank = stk::topology::NODE_RANK)
 {
   const auto* field = meta.get_field(entity_rank, fieldName);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     (field != nullptr), "Requested field does not exist: " + fieldName);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     (field->is_state_valid(state)), "Requested invalid state: " + fieldName);
 
   const auto* fState = field->field_state(state);
@@ -92,8 +92,8 @@ get_node_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
   return stk::mesh::get_updated_ngp_field<T>(
     *meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
diff --git a/src/Algorithm.C b/src/Algorithm.C
index 04b39e296..80be27ec2 100644
--- a/src/Algorithm.C
+++ b/src/Algorithm.C
@@ -9,6 +9,7 @@
 
 #include <Algorithm.h>
 #include <SupplementalAlgorithm.h>
+#include <Realm.h>
 #include <kernel/Kernel.h>
 
 namespace sierra {
@@ -23,14 +24,18 @@ namespace nalu {
 //-------- constructor -----------------------------------------------------
 //--------------------------------------------------------------------------
 Algorithm::Algorithm(Realm& realm, stk::mesh::Part* part)
-  : realm_(realm), partVec_(1, part)
+  : realm_(realm),
+    partVec_(1, part),
+    fieldManager_(*(realm.fieldManager_.get()))
 {
   // nothing to do
 }
 
 // alternative; provide full partVec
 Algorithm::Algorithm(Realm& realm, const stk::mesh::PartVector& partVec)
-  : realm_(realm), partVec_(partVec)
+  : realm_(realm),
+    partVec_(partVec),
+    fieldManager_(*(realm.fieldManager_.get()))
 {
   // nothing to do
 }
diff --git a/src/AssemblePNGBoundarySolverAlgorithm.C b/src/AssemblePNGBoundarySolverAlgorithm.C
index 428d66c95..d8db9935b 100644
--- a/src/AssemblePNGBoundarySolverAlgorithm.C
+++ b/src/AssemblePNGBoundarySolverAlgorithm.C
@@ -146,7 +146,7 @@ AssemblePNGBoundarySolverAlgorithm::execute()
       stk::mesh::Entity const* face_node_rels = b.begin_nodes(k);
       int num_face_nodes = b.num_nodes(k);
       // sanity check on num nodes
-      ThrowAssert(num_face_nodes == nodesPerFace);
+      STK_ThrowAssert(num_face_nodes == nodesPerFace);
       for (int ni = 0; ni < num_face_nodes; ++ni) {
         // get the node and form connected_node
         stk::mesh::Entity node = face_node_rels[ni];
diff --git a/src/AssemblePNGElemSolverAlgorithm.C b/src/AssemblePNGElemSolverAlgorithm.C
index 191815538..afb9cf822 100644
--- a/src/AssemblePNGElemSolverAlgorithm.C
+++ b/src/AssemblePNGElemSolverAlgorithm.C
@@ -177,7 +177,7 @@ AssemblePNGElemSolverAlgorithm::execute()
       int num_nodes = b.num_nodes(k);
 
       // sanity check on num nodes
-      ThrowAssert(num_nodes == nodesPerElement);
+      STK_ThrowAssert(num_nodes == nodesPerElement);
 
       for (int ni = 0; ni < num_nodes; ++ni) {
         stk::mesh::Entity node = node_rels[ni];
diff --git a/src/AssembleScalarEigenEdgeSolverAlgorithm.C b/src/AssembleScalarEigenEdgeSolverAlgorithm.C
index 00321b85f..86147ebd5 100644
--- a/src/AssembleScalarEigenEdgeSolverAlgorithm.C
+++ b/src/AssembleScalarEigenEdgeSolverAlgorithm.C
@@ -214,7 +214,7 @@ AssembleScalarEigenEdgeSolverAlgorithm::execute()
       stk::mesh::Entity const* edge_node_rels = bulk_data.begin_nodes(edge);
 
       // sanity check on number or nodes
-      ThrowAssert(bulk_data.num_nodes(edge) == 2);
+      STK_ThrowAssert(bulk_data.num_nodes(edge) == 2);
 
       // pointer to edge area vector
       for (int j = 0; j < nDim; ++j)
diff --git a/src/AuxFunctionAlgorithm.C b/src/AuxFunctionAlgorithm.C
index bf34b2fb4..4572cc47f 100644
--- a/src/AuxFunctionAlgorithm.C
+++ b/src/AuxFunctionAlgorithm.C
@@ -47,7 +47,7 @@ AuxFunctionAlgorithm::execute()
 {
 
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   stk::mesh::MetaData& meta_data = realm_.meta_data();
 
diff --git a/src/ChienKEpsilonEquationSystem.C b/src/ChienKEpsilonEquationSystem.C
index 8506e79d7..cbd02995f 100644
--- a/src/ChienKEpsilonEquationSystem.C
+++ b/src/ChienKEpsilonEquationSystem.C
@@ -268,7 +268,7 @@ ChienKEpsilonEquationSystem::post_external_data_transfer_work()
   auto tdrBCField = meta.get_field<double>(stk::topology::NODE_RANK, "tdr_bc");
   auto tkeBCField = meta.get_field<double>(stk::topology::NODE_RANK, "tke_bc");
   if (tdrBCField != nullptr) {
-    ThrowRequire(tkeBCField);
+    STK_ThrowRequire(tkeBCField);
     auto bc_sel = owned_and_shared & stk::mesh::selectField(*tdrBCField);
     auto ngpTkeBC =
       fieldMgr.get_field<double>(tkeBCField->mesh_meta_data_ordinal());
diff --git a/src/ComputeHeatTransferEdgeWallAlgorithm.C b/src/ComputeHeatTransferEdgeWallAlgorithm.C
index b2e61bc17..e9768c666 100644
--- a/src/ComputeHeatTransferEdgeWallAlgorithm.C
+++ b/src/ComputeHeatTransferEdgeWallAlgorithm.C
@@ -109,7 +109,7 @@ ComputeHeatTransferEdgeWallAlgorithm::execute()
 
     // extract connected element topology
     b.parent_topology(stk::topology::ELEMENT_RANK, parentTopo);
-    ThrowAssert(parentTopo.size() == 1);
+    STK_ThrowAssert(parentTopo.size() == 1);
     stk::topology theElemTopo = parentTopo[0];
     MasterElement* meSCS =
       sierra::nalu::MasterElementRepo::get_surface_master_element_on_host(
@@ -128,7 +128,7 @@ ComputeHeatTransferEdgeWallAlgorithm::execute()
       // extract the connected element to this exposed face; should be single in
       // size!
       stk::mesh::Entity const* face_elem_rels = b.begin_elements(k);
-      ThrowAssert(b.num_elements(k) == 1);
+      STK_ThrowAssert(b.num_elements(k) == 1);
 
       // get element; its face ordinal number and populate face_node_ordinals
       stk::mesh::Entity element = face_elem_rels[0];
diff --git a/src/ComputeWallFrictionVelocityAlgorithm.C b/src/ComputeWallFrictionVelocityAlgorithm.C
index f089acf6d..aba7cb220 100644
--- a/src/ComputeWallFrictionVelocityAlgorithm.C
+++ b/src/ComputeWallFrictionVelocityAlgorithm.C
@@ -153,7 +153,7 @@ ComputeWallFrictionVelocityAlgorithm::execute()
 
     // extract connected element topology
     b.parent_topology(stk::topology::ELEMENT_RANK, parentTopo);
-    ThrowAssert(parentTopo.size() == 1);
+    STK_ThrowAssert(parentTopo.size() == 1);
     stk::topology theElemTopo = parentTopo[0];
 
     // extract master element
@@ -204,7 +204,7 @@ ComputeWallFrictionVelocityAlgorithm::execute()
       stk::mesh::Entity const* face_node_rels = bulk_data.begin_nodes(face);
       int num_face_nodes = bulk_data.num_nodes(face);
       // sanity check on num nodes
-      ThrowAssert(num_face_nodes == nodesPerFace);
+      STK_ThrowAssert(num_face_nodes == nodesPerFace);
       for (int ni = 0; ni < num_face_nodes; ++ni) {
         stk::mesh::Entity node = face_node_rels[ni];
 
@@ -232,7 +232,7 @@ ComputeWallFrictionVelocityAlgorithm::execute()
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulk_data.begin_elements(face);
-      ThrowAssert(bulk_data.num_elements(face) == 1);
+      STK_ThrowAssert(bulk_data.num_elements(face) == 1);
 
       // get element; its face ordinal number
       stk::mesh::Entity element = face_elem_rels[0];
diff --git a/src/ConstantAuxFunction.C b/src/ConstantAuxFunction.C
index 42e0c980d..ccc4d25b2 100644
--- a/src/ConstantAuxFunction.C
+++ b/src/ConstantAuxFunction.C
@@ -20,7 +20,7 @@ ConstantAuxFunction::ConstantAuxFunction(
   const std::vector<double>& values)
   : AuxFunction(beginPos, endPos), values_(values)
 {
-  ThrowRequire(endPos_ <= values_.size());
+  STK_ThrowRequire(endPos_ <= values_.size());
 }
 
 void
diff --git a/src/CoriolisSrc.C b/src/CoriolisSrc.C
index 027602834..b7ada45f2 100644
--- a/src/CoriolisSrc.C
+++ b/src/CoriolisSrc.C
@@ -35,8 +35,8 @@ CoriolisSrc::CoriolisSrc(const SolutionOptions& solnOpts)
   earthAngularVelocity_ = solnOpts.earthAngularVelocity_;
   latitude_ = solnOpts.latitude_ * pi_ / 180.0;
 
-  ThrowRequire(solnOpts.eastVector_.size() == nDim_);
-  ThrowRequire(solnOpts.northVector_.size() == nDim_);
+  STK_ThrowRequire(solnOpts.eastVector_.size() == nDim_);
+  STK_ThrowRequire(solnOpts.northVector_.size() == nDim_);
 
   for (int i = 0; i < nDim_; ++i) {
     eastVector_[i] = solnOpts.eastVector_[i];
diff --git a/src/DataProbePostProcessing.C b/src/DataProbePostProcessing.C
index 747ef74eb..1d86d2f21 100644
--- a/src/DataProbePostProcessing.C
+++ b/src/DataProbePostProcessing.C
@@ -808,7 +808,7 @@ DataProbePostProcessing::create_exodus()
   for (const auto* probeSpec : dataProbeSpecInfo_) {
     for (const auto& fieldInfo : probeSpec->fieldInfo_) {
       const auto& meta = realm_.meta_data();
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         meta.get_field(stk::topology::NODE_RANK, fieldInfo.first) != nullptr,
         "No field named `" + fieldInfo.first + "' of node rank");
       io->add_field(
diff --git a/src/ElemDataRequests.C b/src/ElemDataRequests.C
index d2efc5fde..001364a7b 100644
--- a/src/ElemDataRequests.C
+++ b/src/ElemDataRequests.C
@@ -16,7 +16,7 @@ void
 ElemDataRequests::add_gathered_nodal_field(
   const stk::mesh::FieldBase& field, unsigned scalarsPerNode)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     field.entity_rank() == stk::topology::NODE_RANK,
     "ElemDataRequests ERROR, add_gathered_nodal_field called with field "
       << field.name() << " which has entity-rank==" << field.entity_rank()
@@ -27,7 +27,7 @@ ElemDataRequests::add_gathered_nodal_field(
   if (iter == fields.end()) {
     fields.insert(fieldInfo);
   } else {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       iter->scalarsDim1 == scalarsPerNode,
       "ElemDataRequests ERROR, gathered-nodal-field "
         << field.name() << " requested with scalarsPerNode==" << scalarsPerNode
@@ -40,7 +40,7 @@ void
 ElemDataRequests::add_gathered_nodal_field(
   const stk::mesh::FieldBase& field, unsigned tensorDim1, unsigned tensorDim2)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     field.entity_rank() == stk::topology::NODE_RANK,
     "ElemDataRequests ERROR, add_gathered_nodal_field called with field "
       << field.name() << " which has entity-rank==" << field.entity_rank()
@@ -51,7 +51,7 @@ ElemDataRequests::add_gathered_nodal_field(
   if (iter == fields.end()) {
     fields.insert(fieldInfo);
   } else {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       iter->scalarsDim1 == tensorDim1 && iter->scalarsDim2 == tensorDim2,
       "ElemDataRequests ERROR, gathered-nodal-field "
         << field.name() << " requested with tensorDim1==" << tensorDim1
@@ -70,7 +70,7 @@ ElemDataRequests::add_ip_field(
   if (iter == fields.end()) {
     fields.insert(fieldInfo);
   } else {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       iter->scalarsDim1 == scalarsPerElement,
       "ElemDataRequests ERROR, ip-field "
         << field.name()
@@ -84,7 +84,7 @@ void
 ElemDataRequests::add_face_field(
   const stk::mesh::FieldBase& field, unsigned scalarsPerFace)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     field.entity_rank() == stk::topology::FACE_RANK ||
       field.entity_rank() == stk::topology::EDGE_RANK,
     "ElemDataRequests ERROR, add_face_field called with field "
@@ -98,7 +98,7 @@ void
 ElemDataRequests::add_element_field(
   const stk::mesh::FieldBase& field, unsigned scalarsPerElement)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     field.entity_rank() == stk::topology::ELEM_RANK,
     "ElemDataRequests ERROR, add_element_field called with field "
       << field.name() << " which has entity-rank==" << field.entity_rank()
@@ -116,7 +116,7 @@ ElemDataRequests::add_ip_field(
   if (iter == fields.end()) {
     fields.insert(fieldInfo);
   } else {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       iter->scalarsDim1 == tensorDim1 && iter->scalarsDim2 == tensorDim2,
       "ElemDataRequests ERROR, ip-field "
         << field.name() << " requested with tensorDim1==" << tensorDim1
@@ -130,7 +130,7 @@ void
 ElemDataRequests::add_face_field(
   const stk::mesh::FieldBase& field, unsigned tensorDim1, unsigned tensorDim2)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     field.entity_rank() == stk::topology::FACE_RANK ||
       field.entity_rank() == stk::topology::EDGE_RANK,
     "ElemDataRequests ERROR, add_face_field called with field "
@@ -144,7 +144,7 @@ void
 ElemDataRequests::add_element_field(
   const stk::mesh::FieldBase& field, unsigned tensorDim1, unsigned tensorDim2)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     field.entity_rank() == stk::topology::ELEM_RANK,
     "ElemDataRequests ERROR, add_element_field called with field "
       << field.name() << " which has entity-rank==" << field.entity_rank()
@@ -168,7 +168,7 @@ ElemDataRequests::add_master_element_call(
   ELEM_DATA_NEEDED data, COORDS_TYPES cType)
 {
   auto it = coordsFields_.find(cType);
-  NGP_ThrowRequireMsg(
+  STK_NGP_ThrowRequireMsg(
     it != coordsFields_.end(),
     "ElemDataRequests:add_master_element_call: Coordinates field "
     "must be registered to ElemDataRequests before registering MasterElement "
@@ -176,13 +176,13 @@ ElemDataRequests::add_master_element_call(
 
   // Check that the appropriate MasterElement has been registered
   if ((data >= BEGIN_FC) && (data <= END_FC)) {
-    NGP_ThrowRequire(meFC_ != nullptr);
+    STK_NGP_ThrowRequire(meFC_ != nullptr);
   } else if ((data >= BEGIN_SCS) && (data <= END_SCS)) {
-    NGP_ThrowRequire(meSCS_ != nullptr);
+    STK_NGP_ThrowRequire(meSCS_ != nullptr);
   } else if ((data >= BEGIN_SCV) && (data <= END_SCV)) {
-    NGP_ThrowRequire(meSCV_ != nullptr);
+    STK_NGP_ThrowRequire(meSCV_ != nullptr);
   } else if ((data >= BEGIN_FEM) && (data <= END_FEM)) {
-    NGP_ThrowRequire(meFEM_ != nullptr);
+    STK_NGP_ThrowRequire(meFEM_ != nullptr);
   }
 
   dataEnums[cType].insert(data);
diff --git a/src/EnthalpyEquationSystem.C b/src/EnthalpyEquationSystem.C
index 01cbffbb9..24f9982ed 100644
--- a/src/EnthalpyEquationSystem.C
+++ b/src/EnthalpyEquationSystem.C
@@ -300,7 +300,7 @@ EnthalpyEquationSystem::register_nodal_fields(
          "this constant value"
       << std::endl;
     Algorithm* propAlg = new ThermalConductivityFromPrandtlPropAlgorithm(
-      realm_, part_vec, thermalCond_, specHeat_, visc_, providedPr);
+      realm_, part_vec, providedPr);
     propertyAlg_.push_back(propAlg);
   } else {
     // no Pr provided, simply augment property map and expect lambda to be
@@ -442,7 +442,7 @@ EnthalpyEquationSystem::register_interior_algorithm(stk::mesh::Part* part)
       [&](AssembleNGPNodeSolverAlgorithm& nodeAlg, std::string& srcName) {
         bool added = true;
         if (srcName == "abl_forcing") {
-          ThrowRequireMsg(
+          STK_ThrowRequireMsg(
             ((NULL != realm_.ablForcingAlg_) &&
              (realm_.ablForcingAlg_->temperatureForcingOn())),
             "ERROR! ABL Forcing parameters not "
diff --git a/src/EquationSystems.C b/src/EquationSystems.C
index 9deac309c..3854ed75c 100644
--- a/src/EquationSystems.C
+++ b/src/EquationSystems.C
@@ -380,7 +380,7 @@ EquationSystems::register_wall_bc(
     const std::vector<stk::mesh::Part*>& mesh_parts = targetPart->subsets();
     for (std::vector<stk::mesh::Part*>::const_iterator i = mesh_parts.begin();
          i != mesh_parts.end(); ++i) {
-      ThrowRequire(*i != nullptr);
+      STK_ThrowRequire(*i != nullptr);
       stk::mesh::Part* const part = *i;
       const stk::topology the_topo = part->topology();
 
diff --git a/src/FieldFunctions.C b/src/FieldFunctions.C
index ff61b81fc..58b8d29d0 100644
--- a/src/FieldFunctions.C
+++ b/src/FieldFunctions.C
@@ -48,7 +48,7 @@ field_axpby(
     stk::mesh::Bucket& b = *buckets[i];
     const stk::mesh::Bucket::size_type length = b.size();
     const size_t fieldSize = field_bytes_per_entity(xField, b) / sizeof(double);
-    ThrowAssert(
+    STK_ThrowAssert(
       fieldSize == field_bytes_per_entity(yField, b) / sizeof(double));
     const unsigned kmax = length * fieldSize;
     const double* x = (double*)stk::mesh::field_data(xField, b);
@@ -145,7 +145,7 @@ field_copy(
     stk::mesh::Bucket& b = *buckets[i];
     const stk::mesh::Bucket::size_type length = b.size();
     const size_t fieldSize = field_bytes_per_entity(xField, b) / sizeof(double);
-    ThrowAssert(
+    STK_ThrowAssert(
       fieldSize == field_bytes_per_entity(yField, b) / sizeof(double));
     const unsigned kmax = length * fieldSize;
     const double* x = (double*)stk::mesh::field_data(xField, b);
diff --git a/src/FieldManager.C b/src/FieldManager.C
index ee52abf26..007666f53 100644
--- a/src/FieldManager.C
+++ b/src/FieldManager.C
@@ -19,7 +19,7 @@ FieldManager::FieldManager(stk::mesh::MetaData& meta, const int numStates)
 }
 
 bool
-FieldManager::field_exists(const std::string& name)
+FieldManager::field_exists(const std::string& name) const
 {
   auto definition = FieldRegistry::query(numDimensions_, numStates_, name);
 
diff --git a/src/FieldRegistry.C b/src/FieldRegistry.C
index 6459181a4..e89efb7f9 100644
--- a/src/FieldRegistry.C
+++ b/src/FieldRegistry.C
@@ -83,7 +83,8 @@ Registry()
     {"pressure",                  SingleStateNodalScalar},
     {"rans_time_scale" ,          SingleStateNodalScalar},
     {"scalarQ",                   SingleStateNodalScalar},
-    {"specific_dissipation_rate", SingleStateNodalScalar},
+    {"specific_dissipation_rate", MultiStateNodalScalar},
+    {"thermal_conductivity",      SingleStateNodalScalar},
     {"specific_heat" ,            SingleStateNodalScalar},
     {"sst_f_one_blending"  ,      SingleStateNodalScalar},
     {"sst_max_length_scale",      SingleStateNodalScalar},
diff --git a/src/FixPressureAtNodeAlgorithm.C b/src/FixPressureAtNodeAlgorithm.C
index 2cb5ca2ab..28feae4b7 100644
--- a/src/FixPressureAtNodeAlgorithm.C
+++ b/src/FixPressureAtNodeAlgorithm.C
@@ -58,7 +58,7 @@ FixPressureAtNodeAlgorithm::execute()
 {
 
   int numNodes = refNodeList_.size();
-  ThrowAssertMsg(
+  STK_ThrowAssertMsg(
     numNodes <= 1,
     "Invalid number of nodes encountered in FixPressureAtNodeAlgorithm");
 
diff --git a/src/HypreLinearSystem.C b/src/HypreLinearSystem.C
index fae0c44be..9b99b986c 100644
--- a/src/HypreLinearSystem.C
+++ b/src/HypreLinearSystem.C
@@ -574,7 +574,7 @@ HypreLinearSystem::buildFaceElemToNodeGraph(const stk::mesh::PartVector& parts)
         // extract the connected element to this exposed face; should be single
         // in size!
         const stk::mesh::Entity* face_elem_rels = bulkData.begin_elements(face);
-        ThrowAssert(bulkData.num_elements(face) == 1);
+        STK_ThrowAssert(bulkData.num_elements(face) == 1);
 
         // get connected element and nodal relations
         stk::mesh::Entity element = face_elem_rels[0];
@@ -603,7 +603,7 @@ HypreLinearSystem::buildFaceElemToNodeGraph(const stk::mesh::PartVector& parts)
         // extract the connected element to this exposed face; should be single
         // in size!
         const stk::mesh::Entity* face_elem_rels = bulkData.begin_elements(face);
-        ThrowAssert(bulkData.num_elements(face) == 1);
+        STK_ThrowAssert(bulkData.num_elements(face) == 1);
 
         // get connected element and nodal relations
         stk::mesh::Entity element = face_elem_rels[0];
@@ -823,7 +823,7 @@ HypreLinearSystem::finalizeLinearSystem()
   gettimeofday(&_start, NULL);
 #endif
 
-  ThrowRequire(inConstruction_);
+  STK_ThrowRequire(inConstruction_);
   inConstruction_ = false;
 
 #ifdef HYPRE_LINEAR_SYSTEM_DEBUG
diff --git a/src/HypreUVWLinearSystem.C b/src/HypreUVWLinearSystem.C
index 0efd8bf44..a03f38065 100644
--- a/src/HypreUVWLinearSystem.C
+++ b/src/HypreUVWLinearSystem.C
@@ -52,7 +52,7 @@ HypreUVWLinearSystem::finalizeLinearSystem()
   gettimeofday(&_start, NULL);
 #endif
 
-  ThrowRequire(inConstruction_);
+  STK_ThrowRequire(inConstruction_);
   inConstruction_ = false;
 
 #ifdef HYPRE_LINEAR_SYSTEM_DEBUG
@@ -1089,7 +1089,7 @@ HypreUVWLinearSystem::buildFaceElemToNodeGraph(
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulkData.begin_elements(face);
-      ThrowAssert(bulkData.num_elements(face) == 1);
+      STK_ThrowAssert(bulkData.num_elements(face) == 1);
 
       // get connected element and nodal relations
       stk::mesh::Entity element = face_elem_rels[0];
diff --git a/src/LinearSolver.C b/src/LinearSolver.C
index 553c032d1..6de400676 100644
--- a/src/LinearSolver.C
+++ b/src/LinearSolver.C
@@ -72,8 +72,8 @@ void
 TpetraLinearSolver::setSystemObjects(
   Teuchos::RCP<LinSys::Matrix> matrix, Teuchos::RCP<LinSys::MultiVector> rhs)
 {
-  ThrowRequire(!matrix.is_null());
-  ThrowRequire(!rhs.is_null());
+  STK_ThrowRequire(!matrix.is_null());
+  STK_ThrowRequire(!rhs.is_null());
 
   matrix_ = matrix;
   rhs_ = rhs;
@@ -169,7 +169,7 @@ TpetraLinearSolver::residual_norm(
 {
   const size_t numVecs = sln->getNumVectors();
   LinSys::MultiVector resid(rhs_->getMap(), numVecs);
-  ThrowRequire(!(sln.is_null() || rhs_.is_null()));
+  STK_ThrowRequire(!(sln.is_null() || rhs_.is_null()));
 
   if (matrix_->isFillActive()) {
     // FIXME
@@ -213,7 +213,7 @@ TpetraLinearSolver::solve(
   double& finalResidNrm,
   bool isFinalOuterIter)
 {
-  ThrowRequire(!sln.is_null());
+  STK_ThrowRequire(!sln.is_null());
 
   const int status = 0;
   int whichNorm = 2;
diff --git a/src/LinearSystem.C b/src/LinearSystem.C
index e02c3bbfd..9e4867f26 100644
--- a/src/LinearSystem.C
+++ b/src/LinearSystem.C
@@ -109,7 +109,7 @@ LinearSystem::useSegregatedSolver() const
 const LinearSolverConfig&
 LinearSystem::config() const
 {
-  ThrowAssert(linearSolver_ != nullptr);
+  STK_ThrowAssert(linearSolver_ != nullptr);
   return *(linearSolver_->getConfig());
 }
 
diff --git a/src/LowMachEquationSystem.C b/src/LowMachEquationSystem.C
index fcfec16fd..6711d3824 100644
--- a/src/LowMachEquationSystem.C
+++ b/src/LowMachEquationSystem.C
@@ -1417,7 +1417,7 @@ MomentumEquationSystem::register_interior_algorithm(stk::mesh::Part* part)
               "MomentumEQS::body_force_box: No force vector found");
           }
         } else if (srcName == "abl_forcing") {
-          ThrowRequireMsg(
+          STK_ThrowRequireMsg(
             ((NULL != realm_.ablForcingAlg_) &&
              (realm_.ablForcingAlg_->momentumForcingOn())),
             "ERROR! ABL Forcing parameters not "
diff --git a/src/MatrixFreeHeatCondEquationSystem.C b/src/MatrixFreeHeatCondEquationSystem.C
index 978aebc1d..117f01cb8 100644
--- a/src/MatrixFreeHeatCondEquationSystem.C
+++ b/src/MatrixFreeHeatCondEquationSystem.C
@@ -37,10 +37,10 @@ MatrixFreeHeatCondEquationSystem::MatrixFreeHeatCondEquationSystem(
     meta_(realm_.meta_data())
 {
   realm_.push_equation_to_systems(this);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     realm_.spatialDimension_ == dim,
     "Only 3D supported for matrix free heat conduction");
-  ThrowRequireMsg(realm_.matrixFree_, "Only matrix free supported");
+  STK_ThrowRequireMsg(realm_.matrixFree_, "Only matrix free supported");
 }
 
 MatrixFreeHeatCondEquationSystem::~MatrixFreeHeatCondEquationSystem() = default;
@@ -53,8 +53,8 @@ get_node_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
   return stk::mesh::get_updated_ngp_field<T>(
     *meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
@@ -126,7 +126,7 @@ void
 MatrixFreeHeatCondEquationSystem::register_interior_algorithm(
   stk::mesh::Part* part)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     matrix_free::part_is_valid_for_matrix_free(polynomial_order_, *part),
     "part " + part->name() + " has invalid topology " +
       part->topology().name() + ". Only hex8/hex27 supported");
@@ -139,7 +139,7 @@ MatrixFreeHeatCondEquationSystem::register_wall_bc(
   const stk::topology&,
   const WallBoundaryConditionData& wallBCData)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     matrix_free::part_is_valid_for_matrix_free(polynomial_order_, *part),
     "part " + part->name() + " has invalid topology " +
       part->topology().name() + ". Only Quad4 and Quad9 supported");
diff --git a/src/MatrixFreeLowMachEquationSystem.C b/src/MatrixFreeLowMachEquationSystem.C
index b068bd19e..416dac115 100644
--- a/src/MatrixFreeLowMachEquationSystem.C
+++ b/src/MatrixFreeLowMachEquationSystem.C
@@ -70,10 +70,10 @@ MatrixFreeLowMachEquationSystem::MatrixFreeLowMachEquationSystem(
     meta_(realm_.meta_data())
 {
   realm_.push_equation_to_systems(this);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     realm_.spatialDimension_ == dim,
     "Only 3D supported for matrix free heat conduction");
-  ThrowRequireMsg(realm_.matrixFree_, "Only matrix free supported");
+  STK_ThrowRequireMsg(realm_.matrixFree_, "Only matrix free supported");
   realm_.hasFluids_ = true;
 }
 
@@ -116,8 +116,8 @@ MatrixFreeLowMachEquationSystem::check_part_is_valid(
   const stk::mesh::PartVector& part_vec)
 {
   for (auto part : part_vec) {
-    ThrowRequire(part);
-    ThrowRequireMsg(
+    STK_ThrowRequire(part);
+    STK_ThrowRequireMsg(
       matrix_free::part_is_valid_for_matrix_free(polynomial_order_, *part),
       "part " + part->name() + " has invalid topology " +
         part->topology().name() + ". Only hex8/hex27 supported");
@@ -146,7 +146,7 @@ MatrixFreeLowMachEquationSystem::register_nodal_fields(
   const std::array<double, 3> x{{0, 0, 0}};
   const int dim = meta_.spatial_dimension();
   check_part_is_valid(part_vec);
-  ThrowRequire(realm_.number_of_states() == 3);
+  STK_ThrowRequire(realm_.number_of_states() == 3);
   constexpr int one_state = 1;
   constexpr int three_states = 3;
   stk::mesh::Selector selector = stk::mesh::selectUnion(part_vec);
@@ -226,7 +226,7 @@ MatrixFreeLowMachEquationSystem::register_wall_bc(
   check_part_is_valid(stk::mesh::PartVector(1, part));
 
   auto data = bc.userData_;
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     !(data.wallFunctionApproach_ || data.ablWallFunctionApproach_),
     "Wall function not implemented");
 
@@ -241,7 +241,8 @@ MatrixFreeLowMachEquationSystem::register_wall_bc(
 
   auto velocity_name = std::string(names::velocity);
   auto bc_data_type = get_bc_data_type(data, velocity_name);
-  ThrowRequireMsg(bc_data_type != FUNCTION_UD, "No user functions yet enabled");
+  STK_ThrowRequireMsg(
+    bc_data_type != FUNCTION_UD, "No user functions yet enabled");
 
   auto* bc_field =
     meta_.get_field(stk::topology::NODE_RANK, names::velocity_bc);
@@ -306,7 +307,7 @@ MatrixFreeLowMachEquationSystem::compute_filter_scale() const
     stk::mesh::for_each_entity_run(
       realm_.ngp_mesh(), stk::topology::NODE_RANK, interior_selector_,
       KOKKOS_LAMBDA(stk::mesh::FastMeshIndex mi) {
-        NGP_ThrowAssert(filter_scale.get(mi, 0) > 0);
+        STK_NGP_ThrowAssert(filter_scale.get(mi, 0) > 0);
         filter_scale.get(mi, 0) =
           scaling * stk::math::cbrt(filter_scale.get(mi, 0));
       });
@@ -324,7 +325,7 @@ MatrixFreeLowMachEquationSystem::register_initial_condition_fcn(
 
   auto it = names.find(names::velocity);
   if (it != names.end()) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       (it->second == "TaylorGreen" || it->second == "SinProfileChannelFlow"),
       "Only TaylorGreen/SinProfileChannelFlow currently implemented for "
       "matrix-free");
@@ -332,7 +333,7 @@ MatrixFreeLowMachEquationSystem::register_initial_condition_fcn(
     if (it->second == "TaylorGreen") {
       auto* velocity_field =
         meta_.get_field<double>(stk::topology::NODE_RANK, names::velocity);
-      ThrowRequire(velocity_field);
+      STK_ThrowRequire(velocity_field);
       auto* vel_func = new TaylorGreenVelocityAuxFunction(0, dim);
       auto* vel_aux_alg = new AuxFunctionAlgorithm(
         realm_, part, velocity_field, vel_func, stk::topology::NODE_RANK);
@@ -340,7 +341,7 @@ MatrixFreeLowMachEquationSystem::register_initial_condition_fcn(
     } else if (it->second == "SinProfileChannelFlow") {
       auto* velocity_field =
         meta_.get_field<double>(stk::topology::NODE_RANK, names::velocity);
-      ThrowRequire(velocity_field);
+      STK_ThrowRequire(velocity_field);
       auto* vel_func = new SinProfileChannelFlowVelocityAuxFunction(0, dim);
       auto* vel_aux_alg = new AuxFunctionAlgorithm(
         realm_, part, velocity_field, vel_func, stk::topology::NODE_RANK);
@@ -350,7 +351,7 @@ MatrixFreeLowMachEquationSystem::register_initial_condition_fcn(
     if (it->second == "TaylorGreen") {
       auto pressure_field =
         meta_.get_field<double>(stk::topology::NODE_RANK, names::pressure);
-      ThrowRequire(pressure_field);
+      STK_ThrowRequire(pressure_field);
       auto* pressure_func = new TaylorGreenPressureAuxFunction();
       auto* pressure_aux_alg = new AuxFunctionAlgorithm(
         realm_, part, pressure_field, pressure_func, stk::topology::NODE_RANK);
@@ -438,8 +439,8 @@ namespace {
 Kokkos::Array<double, 3>
 compute_scaled_gammas(const TimeIntegrator& ti)
 {
-  ThrowRequire(ti.get_time_step() > 0);
-  ThrowRequire(ti.get_gamma1() > 0);
+  STK_ThrowRequire(ti.get_time_step() > 0);
+  STK_ThrowRequire(ti.get_gamma1() > 0);
   return Kokkos::Array<double, 3>{
     {ti.get_gamma1() / ti.get_time_step(), ti.get_gamma2() / ti.get_time_step(),
      ti.get_gamma3() / ti.get_time_step()}};
@@ -448,8 +449,8 @@ compute_scaled_gammas(const TimeIntegrator& ti)
 double
 compute_projected_timescale(const TimeIntegrator& ti)
 {
-  ThrowRequire(ti.get_time_step() > 0);
-  ThrowRequire(ti.get_gamma1() > 0);
+  STK_ThrowRequire(ti.get_time_step() > 0);
+  STK_ThrowRequire(ti.get_gamma1() > 0);
   return ti.get_time_step() / ti.get_gamma1();
 }
 
@@ -501,10 +502,10 @@ MatrixFreeLowMachEquationSystem::get_muelu_xml_file_name()
   const auto block_name =
     equationSystems_.get_solver_block_name(names::pressure);
   auto it = solver_config_map.find(block_name);
-  ThrowRequire(it != solver_config_map.end());
+  STK_ThrowRequire(it != solver_config_map.end());
   auto precond_params = it->second->paramsPrecond();
-  ThrowRequire(precond_params);
-  ThrowRequire(precond_params->isParameter("xml parameter file"));
+  STK_ThrowRequire(precond_params);
+  STK_ThrowRequire(precond_params->isParameter("xml parameter file"));
   return precond_params->get<std::string>("xml parameter file");
 }
 
@@ -695,7 +696,7 @@ MatrixFreeLowMachEquationSystem::compute_body_force() const
     const auto it = realm_.solutionOptions_->srcTermParamMap_.find("momentum");
     if (it != realm_.solutionOptions_->srcTermParamMap_.end()) {
       const auto& force_vec = it->second;
-      ThrowRequireMsg(force_vec.size() == 3u, "Only 3d body force");
+      STK_ThrowRequireMsg(force_vec.size() == 3u, "Only 3d body force");
 
       for (int d = 0; d < 3; ++d) {
         constant_force[d] = force_vec[d];
diff --git a/src/MomentumBoussinesqRASrcNodeSuppAlg.C b/src/MomentumBoussinesqRASrcNodeSuppAlg.C
index c8b880253..4b8517e97 100644
--- a/src/MomentumBoussinesqRASrcNodeSuppAlg.C
+++ b/src/MomentumBoussinesqRASrcNodeSuppAlg.C
@@ -51,11 +51,11 @@ MomentumBoussinesqRASrcNodeSuppAlg::MomentumBoussinesqRASrcNodeSuppAlg(
 
   temperature_ =
     meta_data.get_field<double>(stk::topology::NODE_RANK, "temperature");
-  ThrowRequire(temperature_ != nullptr);
+  STK_ThrowRequire(temperature_ != nullptr);
 
   dualNodalVolume_ =
     meta_data.get_field<double>(stk::topology::NODE_RANK, "dual_nodal_volume");
-  ThrowRequire(dualNodalVolume_ != nullptr);
+  STK_ThrowRequire(dualNodalVolume_ != nullptr);
 
   rhoRef_ = realm_.solutionOptions_->referenceDensity_;
   beta_ = realm_.solutionOptions_->thermalExpansionCoeff_;
@@ -73,7 +73,7 @@ MomentumBoussinesqRASrcNodeSuppAlg::setup()
   raTemperature_ = realm_.meta_data().get_field<double>(
     stk::topology::NODE_RANK,
     MovingAveragePostProcessor::filtered_field_name("temperature"));
-  ThrowRequire(raTemperature_ != nullptr);
+  STK_ThrowRequire(raTemperature_ != nullptr);
 }
 
 //--------------------------------------------------------------------------
diff --git a/src/MovingAveragePostProcessor.C b/src/MovingAveragePostProcessor.C
index 4d2ef8757..6b356751e 100644
--- a/src/MovingAveragePostProcessor.C
+++ b/src/MovingAveragePostProcessor.C
@@ -42,7 +42,7 @@ ExponentialMovingAverager::ExponentialMovingAverager(
 void
 ExponentialMovingAverager::compute_and_set_alpha(double delta_t)
 {
-  ThrowAssert(delta_t > 1.0e6 * std::numeric_limits<double>::min());
+  STK_ThrowAssert(delta_t > 1.0e6 * std::numeric_limits<double>::min());
   alpha_ = (isInit_) ? 1.0 : std::min(1.0, delta_t / timeScale_);
 }
 //--------------------------------------------------------------------------
@@ -71,15 +71,15 @@ MovingAveragePostProcessor::add_fields(std::vector<std::string> fieldNames)
   for (const auto& fieldName : fieldNames) {
     auto& meta = bulk_.mesh_meta_data();
     auto* field = meta.get_field(stk::topology::NODE_RANK, fieldName);
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       field != nullptr,
       "Requested field `" + fieldName + "' not available for averaging");
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       field->type_is<double>(), "Only double precision-typed fields allowed");
 
     stk::mesh::FieldBase* avgField = meta.get_field(
       stk::topology::NODE_RANK, filtered_field_name(field->name()));
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       avgField != nullptr,
       filtered_field_name(field->name()) + " field not registered");
     fieldMap_.insert({field, avgField});
@@ -121,7 +121,7 @@ MovingAveragePostProcessor::execute()
       double* avgFieldVals =
         static_cast<double*>(stk::mesh::field_data(avgField, b));
       const unsigned fieldSize = stk::mesh::field_scalars_per_entity(field, b);
-      ThrowAssert(
+      STK_ThrowAssert(
         fieldSize == stk::mesh::field_scalars_per_entity(avgField, b));
       for (size_t k = 0u; k < b.size(); ++k) {
         const size_t offset = k * fieldSize;
diff --git a/src/NonConformalInfo.C b/src/NonConformalInfo.C
index bc9e1349c..2a82414f9 100644
--- a/src/NonConformalInfo.C
+++ b/src/NonConformalInfo.C
@@ -189,7 +189,7 @@ NonConformalInfo::construct_dgInfo()
 
     // extract connected element topology
     b.parent_topology(stk::topology::ELEMENT_RANK, parentTopo);
-    ThrowAssert(parentTopo.size() == 1);
+    STK_ThrowAssert(parentTopo.size() == 1);
     stk::topology currentElemTopo = parentTopo[0];
 
     // volume and surface master element
@@ -218,7 +218,7 @@ NonConformalInfo::construct_dgInfo()
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulk_data.begin_elements(face);
-      ThrowAssert(bulk_data.num_elements(face) == 1);
+      STK_ThrowAssert(bulk_data.num_elements(face) == 1);
 
       // get element; its face ordinal number
       stk::mesh::Entity element = face_elem_rels[0];
@@ -325,7 +325,7 @@ NonConformalInfo::construct_bounding_points()
     const int num_face_nodes = bulk_data.num_nodes(firstDgInfo->currentFace_);
 
     // sanity check on num nodes
-    ThrowAssert(num_face_nodes == nodesPerFace);
+    STK_ThrowAssert(num_face_nodes == nodesPerFace);
     for (int ni = 0; ni < num_face_nodes; ++ni) {
       stk::mesh::Entity node = face_node_rels[ni];
       double* coords = stk::mesh::field_data(*coordinates, node);
@@ -516,7 +516,7 @@ NonConformalInfo::determine_elems_to_ghost()
 
       // extract the connected element
       const stk::mesh::Entity* face_elem_rels = bulk_data.begin_elements(face);
-      ThrowAssert(bulk_data.num_elements(face) == 1);
+      STK_ThrowAssert(bulk_data.num_elements(face) == 1);
       stk::mesh::Entity element = face_elem_rels[0];
 
       // deal with elements to push back to be ghosted; downward relations come
@@ -624,7 +624,7 @@ NonConformalInfo::complete_search()
             // extract the connected element to the opposing face
             const stk::mesh::Entity* face_elem_rels =
               bulk_data.begin_elements(opposingFace);
-            ThrowAssert(bulk_data.num_elements(opposingFace) == 1);
+            STK_ThrowAssert(bulk_data.num_elements(opposingFace) == 1);
             stk::mesh::Entity opposingElement = face_elem_rels[0];
 
             // extract the opposing element topo and associated master element
diff --git a/src/PeriodicManager.C b/src/PeriodicManager.C
index c91b82cca..388669f1d 100644
--- a/src/PeriodicManager.C
+++ b/src/PeriodicManager.C
@@ -817,7 +817,7 @@ PeriodicManager::ngp_periodic_parallel_communicate_field(
     }
 #endif
     else {
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         false, "Error, unsupported field type in "
                "PeriodicManager::periodic_parallel_communicate_field.");
     }
@@ -882,7 +882,7 @@ PeriodicManager::ngp_parallel_communicate_field(
     }
 #endif
     else {
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         false, "Error, field "
                  << theField->name() << " with unsupported type ("
                  << theField->data_traits().type_info.name()
@@ -1070,7 +1070,7 @@ PeriodicManager::ngp_add_slave_to_master(
     ngp_periodic_parallel_communicate_field(theField);
   }
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     theField->type_is<double>(),
     "Error in PeriodicManager::add_slave_to_master, theField ("
       << theField->name() << ") is required to be double.");
@@ -1202,7 +1202,7 @@ PeriodicManager::ngp_set_slave_to_master(
     ngp_periodic_parallel_communicate_field(theField);
   }
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     theField->type_is<double>(),
     "Argh, theField (" << theField->name() << ") is not double.");
 
diff --git a/src/Realm.C b/src/Realm.C
index 42519e33b..aefef81e7 100644
--- a/src/Realm.C
+++ b/src/Realm.C
@@ -1153,7 +1153,7 @@ Realm::setup_bc()
         name, *reinterpret_cast<const ABLTopBoundaryConditionData*>(bc.get()));
       break;
     case PERIODIC_BC: {
-      ThrowAssert(
+      STK_ThrowAssert(
         reinterpret_cast<const PeriodicBoundaryConditionData*>(bc.get()) !=
         nullptr);
       const auto& pbc =
@@ -1295,13 +1295,13 @@ Realm::setup_initial_conditions()
         const ConstantInitialConditionData& genIC =
           *reinterpret_cast<const ConstantInitialConditionData*>(
             initCond.get());
-        ThrowAssert(genIC.data_.size() == genIC.fieldNames_.size());
+        STK_ThrowAssert(genIC.data_.size() == genIC.fieldNames_.size());
         for (size_t ifield = 0; ifield < genIC.fieldNames_.size(); ++ifield) {
 
           std::vector<double> genSpec = genIC.data_[ifield];
           stk::mesh::FieldBase* field = stk::mesh::get_field_by_name(
             genIC.fieldNames_[ifield], meta_data());
-          ThrowAssert(field);
+          STK_ThrowAssert(field);
 
           stk::mesh::FieldBase* fieldWithState =
             (field->number_of_states() > 1)
@@ -1842,7 +1842,7 @@ Realm::makeSureNodesHaveValidTopology()
   // now we require all nodes are in proper node part
   if (nodes_vector.size())
     std::cout << "nodes_vector= " << nodes_vector.size() << std::endl;
-  ThrowRequire(0 == nodes_vector.size());
+  STK_ThrowRequire(0 == nodes_vector.size());
 }
 
 void
@@ -2217,7 +2217,7 @@ Realm::create_restart_mesh()
     }
 
     // set max size for restart data base
-    ioBroker_->get_output_io_region(restartFileIndex_)
+    ioBroker_->get_output_ioss_region(restartFileIndex_)
       ->get_database()
       ->set_cycle_count(outputInfo_->restartMaxDataBaseStepSize_);
   }
@@ -2991,7 +2991,7 @@ void
 Realm::register_non_conformal_bc(
   stk::mesh::Part* part, const stk::topology& theTopo)
 {
-  ThrowRequire(!matrixFree_);
+  STK_ThrowRequire(!matrixFree_);
 
   // push back the part for book keeping and, later, skin mesh
   bcPartVec_.push_back(part);
@@ -3708,7 +3708,7 @@ Realm::check_job(bool get_node_count)
   // set number of nodes, check job run size
   if (get_node_count) {
     size_t localNodeCount =
-      ioBroker_->get_input_io_region()->get_property("node_count").get_int();
+      ioBroker_->get_input_ioss_region()->get_property("node_count").get_int();
     stk::all_reduce_sum(
       NaluEnv::self().parallel_comm(), &localNodeCount, &nodeCount_, 1);
     NaluEnv::self().naluOutputP0()
@@ -3782,7 +3782,7 @@ Realm::check_job(bool get_node_count)
   if (meta_data().is_commit()) {
     std::vector<size_t> counts;
     stk::mesh::comm_mesh_counts(*bulkData_, counts);
-    ThrowRequire(counts.size() >= 4);
+    STK_ThrowRequire(counts.size() >= 4);
     size_t nodeCount = counts[stk::topology::NODE_RANK];
     size_t edgeCount = counts[stk::topology::EDGE_RANK];
     size_t faceCount = counts[stk::topology::FACE_RANK];
@@ -3792,10 +3792,10 @@ Realm::check_job(bool get_node_count)
     unsigned nfields = fields.size();
     for (unsigned ifld = 0; ifld < nfields; ++ifld) {
       stk::mesh::FieldBase* field = fields[ifld];
-      unsigned fszNode = field->max_size(stk::topology::NODE_RANK);
-      unsigned fszEdge = field->max_size(stk::topology::EDGE_RANK);
-      unsigned fszFace = field->max_size(stk::topology::FACE_RANK);
-      unsigned fszElem = field->max_size(stk::topology::ELEM_RANK);
+      unsigned fszNode = field->max_size();
+      unsigned fszEdge = field->max_size();
+      unsigned fszFace = field->max_size();
+      unsigned fszElem = field->max_size();
 
       memoryEstimateFields += (nodeCount * fszNode + edgeCount * fszEdge +
                                faceCount * fszFace + elemCount * fszElem) *
@@ -5164,7 +5164,7 @@ Realm::handle_all_element_part_alias(
   if (names.size() == 1u && names.at(0) == allElementPartAlias) {
     std::vector<std::string> new_names;
     for (const auto* part : meta_data().get_mesh_parts()) {
-      ThrowRequire(part);
+      STK_ThrowRequire(part);
       if (part->topology().rank() == stk::topology::ELEMENT_RANK) {
         new_names.push_back(part->name());
       }
diff --git a/src/ScratchViews.C b/src/ScratchViews.C
index 5bb3c2d0c..60e5e0264 100644
--- a/src/ScratchViews.C
+++ b/src/ScratchViews.C
@@ -41,7 +41,7 @@ gather_elem_node_tensor_field(
   const stk::mesh::NgpMesh::ConnectedNodes& elemNodes,
   ViewType& shmemView)
 {
-  NGP_ThrowRequireMsg(
+  STK_NGP_ThrowRequireMsg(
     numNodes == (int)elemNodes.size(),
     "gather_elem_node_tensor_field, numNodes = mismatch with elemNodes.size()");
   for (int i = 0; i < numNodes; ++i) {
@@ -218,20 +218,20 @@ get_num_scalars_pre_req_data(
 
   switch (reqType) {
   case ElemReqType::ELEM:
-    NGP_ThrowRequireMsg(
+    STK_NGP_ThrowRequireMsg(
       hasElemME, "Requesting ELEM data, but no ELEM_RANK master element has "
                  "been registered");
     break;
 
   case ElemReqType::FACE:
-    NGP_ThrowRequireMsg(
+    STK_NGP_ThrowRequireMsg(
       hasFaceME || hasSCS, "Request SIDE_RANK data, but no SIDE_RANK master "
                            "element has been registered");
     break;
 
   case ElemReqType::FACE_ELEM:
     // In case of FACE_ELEM register meFC so that numFaceIp can be queried
-    NGP_ThrowRequireMsg(
+    STK_NGP_ThrowRequireMsg(
       (hasSCS && hasFaceME),
       "Requesting FACE_ELEM data but does not have necessary MasterElements");
     break;
@@ -440,7 +440,7 @@ fill_pre_req_data(
         }
       }
     } else {
-      NGP_ThrowRequireMsg(
+      STK_NGP_ThrowRequireMsg(
         false, "Unknown stk-rank in ScratchViewsNGP.C::fill_pre_req_data");
     }
   }
diff --git a/src/ShearStressTransportEquationSystem.C b/src/ShearStressTransportEquationSystem.C
index 24aa93354..c3bef8992 100644
--- a/src/ShearStressTransportEquationSystem.C
+++ b/src/ShearStressTransportEquationSystem.C
@@ -368,7 +368,7 @@ ShearStressTransportEquationSystem::post_external_data_transfer_work()
   auto tkeBCField = meta.get_field<double>(stk::topology::NODE_RANK, "tke_bc");
 
   if (sdrBCField != nullptr) {
-    ThrowRequire(tkeBCField);
+    STK_ThrowRequire(tkeBCField);
     auto bc_sel = owned_and_shared & stk::mesh::selectField(*sdrBCField);
     auto ngpTkeBC =
       fieldMgr.get_field<double>(tkeBCField->mesh_meta_data_ordinal());
diff --git a/src/SideWriter.C b/src/SideWriter.C
index a32a3d138..bdbfa0eb0 100644
--- a/src/SideWriter.C
+++ b/src/SideWriter.C
@@ -76,7 +76,8 @@ get_coordinate_field(const stk::mesh::BulkData& bulk)
 {
   auto* coord = dynamic_cast<const stk::mesh::Field<double>*>(
     bulk.mesh_meta_data().coordinate_field());
-  ThrowRequireMsg(coord, "Model coordinates must be a double-valued vector");
+  STK_ThrowRequireMsg(
+    coord, "Model coordinates must be a double-valued vector");
   return *coord;
 }
 
@@ -92,7 +93,7 @@ write_element_block_definition(
       auto block = std::make_unique<Ioss::ElementBlock>(
         &io, subset->name(), subset->topology().name(),
         count_faces(bulk, *subset));
-      ThrowRequire(block);
+      STK_ThrowRequire(block);
       region.add(block.release());
     }
   }
@@ -111,7 +112,7 @@ write_node_block_definition(
   const auto nnodes = count_ent(buckets);
   const int dim = get_dimension(bulk);
   auto block = std::make_unique<Ioss::NodeBlock>(&io, block_name, nnodes, dim);
-  ThrowRequire(block);
+  STK_ThrowRequire(block);
   region.add(block.release());
 }
 
@@ -192,7 +193,7 @@ write_element_connectivity(
         }
       }
       auto block = region.get_element_block(subset->name());
-      ThrowRequire(block);
+      STK_ThrowRequire(block);
       block->put_field_data("ids", ids);
       block->put_field_data("connectivity", connectivity);
     }
@@ -207,14 +208,14 @@ put_data_on_node_block(
   const std::vector<int64_t>& ids,
   const stk::mesh::FieldBase& field)
 {
-  ThrowRequire(field.type_is<T>());
-  const int max_size = field.max_size(stk::topology::NODE_RANK);
+  STK_ThrowRequire(field.type_is<T>());
+  const int max_size = field.max_size();
   std::vector<T> flat_array(ids.size() * max_size, 0);
   for (decltype(ids.size()) k = 0; k < ids.size(); ++k) {
     const auto node = bulk.get_entity(stk::topology::NODE_RANK, ids[k]);
     const T* field_data = static_cast<T*>(stk::mesh::field_data(field, node));
     if (field_data) {
-      ThrowRequire(
+      STK_ThrowRequire(
         stk::mesh::field_scalars_per_entity(field, node) ==
         static_cast<unsigned>(max_size));
       for (int j = 0; j < max_size; ++j) {
@@ -250,7 +251,7 @@ SideWriter::SideWriter(
 
   auto database = Ioss::IOFactory::create(
     "exodus", fname, Ioss::WRITE_RESULTS, bulk.parallel(), prop);
-  ThrowRequire(database != nullptr && database->ok(true));
+  STK_ThrowRequire(database != nullptr && database->ok(true));
   output_ = std::make_unique<Ioss::Region>(database, "SideOutput");
 
   const std::string node_block_name("side_nodes");
@@ -287,7 +288,7 @@ SideWriter::write_database_data(double time)
     output_->begin_state(current_output_step);
     {
       for (auto* block : output_->get_node_blocks()) {
-        ThrowRequire(block);
+        STK_ThrowRequire(block);
         std::vector<int64_t> ids;
         block->get_field_data("ids", ids);
         for (const auto* field : fields_) {
@@ -309,9 +310,10 @@ SideWriter::add_fields(std::vector<const stk::mesh::FieldBase*> fields)
 
   for (auto* block : output_->get_node_blocks()) {
     for (const auto* field : fields_) {
-      ThrowRequireMsg(field->type_is<double>(), "only double fields supported");
+      STK_ThrowRequireMsg(
+        field->type_is<double>(), "only double fields supported");
       const size_t nb_size = block->get_property("entity_count").get_int();
-      switch (field->max_size(stk::topology::NODE_RANK)) {
+      switch (field->max_size()) {
       case 1: {
         Ioss::Field ioss_field(
           field->name(), Ioss::Field::DOUBLE, "scalar", Ioss::Field::TRANSIENT,
diff --git a/src/SolutionNormPostProcessing.C b/src/SolutionNormPostProcessing.C
index 9d70f7dda..ba9f21a9b 100644
--- a/src/SolutionNormPostProcessing.C
+++ b/src/SolutionNormPostProcessing.C
@@ -181,7 +181,7 @@ SolutionNormPostProcessing::setup()
 
     // find the size; is there a better wat to determine a field size on a given
     // part?
-    const int dofSize = dofField->max_size(stk::topology::NODE_RANK);
+    const int dofSize = dofField->max_size();
 
     // increment total dof component size
     totalDofCompSize_ += dofSize;
diff --git a/src/SurfaceForceAndMomentAlgorithm.C b/src/SurfaceForceAndMomentAlgorithm.C
index c685cd151..c1e509b65 100644
--- a/src/SurfaceForceAndMomentAlgorithm.C
+++ b/src/SurfaceForceAndMomentAlgorithm.C
@@ -223,7 +223,7 @@ SurfaceForceAndMomentAlgorithm::execute()
 
     // extract connected element topology
     b.parent_topology(stk::topology::ELEMENT_RANK, parentTopo);
-    ThrowAssert(parentTopo.size() == 1);
+    STK_ThrowAssert(parentTopo.size() == 1);
     stk::topology theElemTopo = parentTopo[0];
 
     // extract master element for this element topo
@@ -277,7 +277,7 @@ SurfaceForceAndMomentAlgorithm::execute()
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulk_data.begin_elements(face);
-      ThrowAssert(bulk_data.num_elements(face) == 1);
+      STK_ThrowAssert(bulk_data.num_elements(face) == 1);
 
       // get element; its face ordinal number
       stk::mesh::Entity element = face_elem_rels[0];
diff --git a/src/TpetraLinearSystem.C b/src/TpetraLinearSystem.C
index 17236a5fa..1d998ff26 100644
--- a/src/TpetraLinearSystem.C
+++ b/src/TpetraLinearSystem.C
@@ -180,7 +180,7 @@ TpetraLinearSystem::beginLinearSystemConstruction()
   if (inConstruction_)
     return;
   inConstruction_ = true;
-  ThrowRequire(ownedGraph_.is_null());
+  STK_ThrowRequire(ownedGraph_.is_null());
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
   stk::mesh::MetaData& metaData = realm_.meta_data();
 
@@ -293,7 +293,7 @@ TpetraLinearSystem::beginLinearSystemConstruction()
       ownedGids.push_back(basegid + idof);
     gstart++;
   }
-  ThrowRequire(localId == numOwnedNodes);
+  STK_ThrowRequire(localId == numOwnedNodes);
   // communicate the newly stored GID's.
 
   std::vector<const stk::mesh::FieldBase*> fVec{realm_.tpetGlobalId_};
@@ -391,7 +391,7 @@ TpetraLinearSystem::insert_connection(stk::mesh::Entity a, stk::mesh::Entity b)
 {
   size_t idx = entityToLIDHost_[a.local_offset()] / numDof_;
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     idx < ownedAndSharedNodes_.size(),
     "Error, insert_connection got index out of range.");
 
@@ -406,7 +406,7 @@ TpetraLinearSystem::insert_connection(stk::mesh::Entity a, stk::mesh::Entity b)
     correctEntity =
       ownedAndSharedNodes_[idx] == master || naluid_a == naluid_master;
   }
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     correctEntity,
     "Error, indexing of rowEntities to connections isn't right.");
 
@@ -589,11 +589,11 @@ TpetraLinearSystem::buildSparsifiedEdgeElemToNodeGraph(
   for (const auto* ib : buckets) {
     const auto& b = *ib;
     if (poly == 1) {
-      ThrowRequire(b.topology() == stk::topology::HEX_8);
+      STK_ThrowRequire(b.topology() == stk::topology::HEX_8);
     } else if (poly == 2) {
-      ThrowRequire(b.topology() == stk::topology::HEX_27);
+      STK_ThrowRequire(b.topology() == stk::topology::HEX_27);
     } else {
-      ThrowRequire(b.topology().is_superelement());
+      STK_ThrowRequire(b.topology().is_superelement());
     }
     for (size_t k = 0u; k < b.size(); ++k) {
       stk::mesh::Entity const* elem_nodes = b.begin_nodes(k);
@@ -649,7 +649,7 @@ TpetraLinearSystem::buildFaceElemToNodeGraph(const stk::mesh::PartVector& parts)
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulkData.begin_elements(face);
-      ThrowAssert(bulkData.num_elements(face) == 1);
+      STK_ThrowAssert(bulkData.num_elements(face) == 1);
 
       // get connected element and nodal relations
       stk::mesh::Entity element = face_elem_rels[0];
@@ -765,8 +765,8 @@ TpetraLinearSystem::copy_stk_to_tpetra(
   const stk::mesh::FieldBase* stkField,
   const Teuchos::RCP<LinSys::MultiVector> tpetraField)
 {
-  ThrowAssert(!tpetraField.is_null());
-  ThrowAssert(stkField);
+  STK_ThrowAssert(!tpetraField.is_null());
+  STK_ThrowAssert(stkField);
   const int numVectors = tpetraField->getNumVectors();
 
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
@@ -786,7 +786,7 @@ TpetraLinearSystem::copy_stk_to_tpetra(
     const int fieldSize =
       field_bytes_per_entity(*stkField, b) / (sizeof(double));
 
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       numVectors == fieldSize, "TpetraLinearSystem::copy_stk_to_tpetra");
 
     const stk::mesh::Bucket::size_type length = b.size();
@@ -802,7 +802,7 @@ TpetraLinearSystem::copy_stk_to_tpetra(
 
       const stk::mesh::EntityId nodeTpetGID =
         *stk::mesh::field_data(*realm_.tpetGlobalId_, node);
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         nodeTpetGID != 0 &&
           nodeTpetGID != static_cast<stk::mesh::EntityId>(
                            std::numeric_limits<LinSys::GlobalOrdinal>::max()),
@@ -979,7 +979,7 @@ TpetraLinearSystem::insert_graph_connections(
 
     const stk::mesh::Entity entity_a = rowEntities[i];
     int dofStatus_a = getDofStatus(entity_a);
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       entityToColLIDHost_[entity_a.local_offset()] != -1,
       "insert_graph_connections bad lid ");
     localDofs_a[0] = entityToColLIDHost_[entity_a.local_offset()];
@@ -987,7 +987,7 @@ TpetraLinearSystem::insert_graph_connections(
       const stk::mesh::Entity entity_b = entities_b[j];
       dofStatus[j] = getDofStatus(entity_b);
 
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         entityToColLIDHost_[entity_b.local_offset()] != -1,
         "insert_graph_connections bad lid #2 ");
 
@@ -1123,7 +1123,7 @@ TpetraLinearSystem::storeOwnersForShared()
         stk::mesh::Entity master = get_entity_master(bulkData, node, naluId);
         GlobalOrdinal gidbase =
           *stk::mesh::field_data(*realm_.tpetGlobalId_, master);
-        ThrowRequire(gidbase != 0);
+        STK_ThrowRequire(gidbase != 0);
         for (unsigned idof = 0; idof < numDof_; ++idof) {
           GlobalOrdinal gid = gidbase + idof;
           ownersAndGids_.insert(
@@ -1137,7 +1137,7 @@ TpetraLinearSystem::storeOwnersForShared()
 void
 TpetraLinearSystem::finalizeLinearSystem()
 {
-  ThrowRequire(inConstruction_);
+  STK_ThrowRequire(inConstruction_);
   inConstruction_ = false;
 
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
@@ -1274,10 +1274,10 @@ TpetraLinearSystem::finalizeLinearSystem()
 void
 TpetraLinearSystem::zeroSystem()
 {
-  ThrowRequire(!ownedMatrix_.is_null());
-  ThrowRequire(!sharedNotOwnedMatrix_.is_null());
-  ThrowRequire(!sharedNotOwnedRhs_.is_null());
-  ThrowRequire(!ownedRhs_.is_null());
+  STK_ThrowRequire(!ownedMatrix_.is_null());
+  STK_ThrowRequire(!sharedNotOwnedMatrix_.is_null());
+  STK_ThrowRequire(!sharedNotOwnedRhs_.is_null());
+  STK_ThrowRequire(!ownedRhs_.is_null());
 
   sharedNotOwnedMatrix_->resumeFill();
   ownedMatrix_->resumeFill();
@@ -1367,7 +1367,8 @@ sum_into_row(
     }
 
     if (offset < length) {
-      ThrowAssertMsg(std::isfinite(input_values[perm_index]), "Inf or NAN lhs");
+      STK_ThrowAssertMsg(
+        std::isfinite(input_values[perm_index]), "Inf or NAN lhs");
       if (forceAtomic) {
         Kokkos::atomic_add(&(row_view.value(offset)), input_values[perm_index]);
       } else {
@@ -1431,7 +1432,7 @@ sum_into(
     const LocalOrdinal cur_perm_index = sortPermutation[r];
     const double* const cur_lhs = &lhs(cur_perm_index, 0);
     const double cur_rhs = rhs[cur_perm_index];
-    //    ThrowAssertMsg(std::isfinite(cur_rhs), "Inf or NAN rhs");
+    //    STK_ThrowAssertMsg(std::isfinite(cur_rhs), "Inf or NAN rhs");
 
     if (rowLid < maxOwnedRowId) {
       sum_into_row(
@@ -1507,7 +1508,7 @@ reset_rows(
       const LocalOrdinal actualLocalId =
         useOwned ? localId : (localId - maxOwnedRowId);
 
-      NGP_ThrowRequire(localId <= maxSharedNotOwnedRowId);
+      STK_NGP_ThrowRequireMsg(localId <= maxSharedNotOwnedRowId, "Error");
 
       // Adjust the LHS; zero out all entries (including diagonal)
       reset_row(localMatrix.row(actualLocalId), actualLocalId, diag_value);
@@ -1595,10 +1596,11 @@ TpetraLinearSystem::sumInto(
   const SharedMemView<int*, DeviceShmem>& sortPermutation,
   const char* /* trace_tag */)
 {
-  ThrowAssertMsg(lhs.span_is_contiguous(), "LHS assumed contiguous");
-  ThrowAssertMsg(rhs.span_is_contiguous(), "RHS assumed contiguous");
-  ThrowAssertMsg(localIds.span_is_contiguous(), "localIds assumed contiguous");
-  ThrowAssertMsg(
+  STK_ThrowAssertMsg(lhs.span_is_contiguous(), "LHS assumed contiguous");
+  STK_ThrowAssertMsg(rhs.span_is_contiguous(), "RHS assumed contiguous");
+  STK_ThrowAssertMsg(
+    localIds.span_is_contiguous(), "localIds assumed contiguous");
+  STK_ThrowAssertMsg(
     sortPermutation.span_is_contiguous(), "sortPermutation assumed contiguous");
 
   sum_into(
@@ -1620,15 +1622,15 @@ TpetraLinearSystem::sumInto(
   const size_t n_obj = entities.size();
   const unsigned numRows = n_obj * numDof_;
 
-  ThrowAssert(numRows == rhs.size());
-  ThrowAssert(numRows * numRows == lhs.size());
+  STK_ThrowAssert(numRows == rhs.size());
+  STK_ThrowAssert(numRows * numRows == lhs.size());
 
   scratchIds.resize(numRows);
   sortPermutation_.resize(numRows);
   for (size_t i = 0; i < n_obj; i++) {
     const stk::mesh::Entity entity = entities[i];
     const LocalOrdinal localOffset = entityToColLIDHost_[entity.local_offset()];
-    ThrowRequireMsg(localOffset != -1, "sumInto bad lid #2 ");
+    STK_ThrowRequireMsg(localOffset != -1, "sumInto bad lid #2 ");
     for (size_t d = 0; d < numDof_; ++d) {
       size_t lid = i * numDof_ + d;
       scratchIds[lid] = localOffset + d;
@@ -1648,7 +1650,7 @@ TpetraLinearSystem::sumInto(
     const LocalOrdinal cur_perm_index = sortPermutation_[r];
     const double* const cur_lhs = &lhs[cur_perm_index * numRows];
     const double cur_rhs = rhs[cur_perm_index];
-    ThrowAssertMsg(std::isfinite(cur_rhs), "Invalid rhs");
+    STK_ThrowAssertMsg(std::isfinite(cur_rhs), "Invalid rhs");
 
     if (rowLid < maxOwnedRowId_) {
       sum_into_row(
@@ -1738,7 +1740,7 @@ TpetraLinearSystem::applyDirichletBCs(
         const LocalOrdinal actualLocalId =
           useOwned ? localId : localId - maxOwnedRowId;
 
-        NGP_ThrowAssert(localId <= maxSharedNotOwnedRowId);
+        STK_NGP_ThrowAssert(localId <= maxSharedNotOwnedRowId);
 
         adjust_lhs_row(
           local_matrix.row(actualLocalId), actualLocalId, diagonalValue);
@@ -2186,8 +2188,8 @@ TpetraLinearSystem::copy_tpetra_to_stk(
 
   const stk::mesh::MetaData& metaData = realm_.meta_data();
 
-  ThrowAssert(!tpetraField.is_null());
-  ThrowAssert(stkField);
+  STK_ThrowAssert(!tpetraField.is_null());
+  STK_ThrowAssert(stkField);
   const auto deviceVector =
     tpetraField->getLocalViewDevice(Tpetra::Access::ReadWrite);
 
@@ -2212,7 +2214,7 @@ TpetraLinearSystem::copy_tpetra_to_stk(
       const LocalOrdinal localIdOffset = entityToLID[node.local_offset()];
       for (unsigned d = 0; d < numDof; ++d) {
         const LocalOrdinal localId = localIdOffset + d;
-        NGP_ThrowRequire(localId < maxOwnedRowId);
+        STK_NGP_ThrowRequireMsg(localId < maxOwnedRowId, "Error");
 
         ngpField.get(meshIdx, d) = deviceVector(localId, 0);
       }
diff --git a/src/TpetraLinearSystemHelpers.C b/src/TpetraLinearSystemHelpers.C
index 57672826b..eaf2a93b7 100644
--- a/src/TpetraLinearSystemHelpers.C
+++ b/src/TpetraLinearSystemHelpers.C
@@ -86,7 +86,7 @@ fill_owned_and_shared_then_nonowned_ordered_by_proc(
   for (unsigned i = 0; i < sharedIndices.size(); ++i) {
     totalGids.push_back(sharedIndices[i]);
     srcPids.push_back(sharedPids[i]);
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       sharedPids[i] != localProc && sharedPids[i] >= 0,
       "Error, bad sharedPid = " << sharedPids[i] << ", localProc = "
                                 << localProc << ", gid = " << sharedIndices[i]);
@@ -101,13 +101,13 @@ fill_owned_and_shared_then_nonowned_ordered_by_proc(
       !sharedNotOwnedRowsMap->isNodeGlobalElement(gid)) {
       totalGids.push_back(gid);
       srcPids.push_back(procAndGid.first);
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         procAndGid.first != localProc && procAndGid.first >= 0,
         "Error, bad remote proc = " << procAndGid.first);
     }
   }
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     srcPids.size() == (totalGids.size() - ownedIndices.size()),
     "Error, bad srcPids.size() = " << srcPids.size());
 }
@@ -134,7 +134,7 @@ get_entity_master(
          << bulk.bucket(elems[i]).topology()
          << ",owned=" << bulk.bucket(elems[i]).owned() << "}";
     }
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       bulk.is_valid(master),
       "get_entity_master, P"
         << bulk.parallel_rank() << " failed to get entity for naluId=" << naluId
@@ -150,7 +150,7 @@ get_neighbor_index(const std::vector<int>& neighborProcs, int proc)
 {
   std::vector<int>::const_iterator neighbor =
     std::find(neighborProcs.begin(), neighborProcs.end(), proc);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     neighbor != neighborProcs.end(),
     "Error, failed to find p=" << proc << " in neighborProcs.");
 
@@ -224,7 +224,7 @@ add_lengths_to_comm_tpet(
   stk::CommBufferV& sbuf = commNeighbors.send_buffer(owner);
   const auto node = bulk.get_entity(stk::topology::NODE_RANK, entityId_a);
   LinSys::GlobalOrdinal rowGid = *stk::mesh::field_data(*tpetGID_label, node);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     rowGid != 0 && rowGid != std::numeric_limits<LinSys::GlobalOrdinal>::max(),
     "add_lengths_to_comm_tpet");
   sbuf.pack(rowGid);
@@ -234,7 +234,7 @@ add_lengths_to_comm_tpet(
       bulk.get_entity(stk::topology::NODE_RANK, colEntityIds[c]);
     LinSys::GlobalOrdinal colGid0 =
       *stk::mesh::field_data(*tpetGID_label, centity);
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       colGid0 != 0 &&
         colGid0 != std::numeric_limits<LinSys::GlobalOrdinal>::max(),
       "add_lengths_to_comm_tpet");
@@ -261,7 +261,7 @@ communicate_remote_columns(
       LinSys::GlobalOrdinal rowGid = 0;
       rbuf.unpack(rowGid);
 
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         rowGid != 0 &&
           rowGid != std::numeric_limits<LinSys::GlobalOrdinal>::max(),
         "communicate_remote_columns");
@@ -282,7 +282,7 @@ communicate_remote_columns(
         LinSys::GlobalOrdinal colGid = 0;
         rbuf.unpack(colGid);
 
-        ThrowRequireMsg(
+        STK_ThrowRequireMsg(
           colGid != 0 &&
             colGid != std::numeric_limits<LinSys::GlobalOrdinal>::max(),
           "communicate_remote_columns");
@@ -329,7 +329,7 @@ insert_communicated_col_indices(
       stk::mesh::EntityId rowGid = 0;
       rbuf.unpack(rowGid);
 
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         rowGid != 0 &&
           rowGid != static_cast<stk::mesh::EntityId>(
                       std::numeric_limits<LinSys::GlobalOrdinal>::max()),
@@ -344,7 +344,7 @@ insert_communicated_col_indices(
         GlobalOrdinal colGid = 0;
         rbuf.unpack(colGid);
 
-        ThrowRequireMsg(
+        STK_ThrowRequireMsg(
           colGid != 0 &&
             colGid != std::numeric_limits<LinSys::GlobalOrdinal>::max(),
           " insert_communicated_col_indices");
diff --git a/src/TpetraSegregatedLinearSystem.C b/src/TpetraSegregatedLinearSystem.C
index 023388ada..fade166aa 100644
--- a/src/TpetraSegregatedLinearSystem.C
+++ b/src/TpetraSegregatedLinearSystem.C
@@ -176,7 +176,7 @@ TpetraSegregatedLinearSystem::beginLinearSystemConstruction()
   if (inConstruction_)
     return;
   inConstruction_ = true;
-  ThrowRequire(ownedGraph_.is_null());
+  STK_ThrowRequire(ownedGraph_.is_null());
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
   stk::mesh::MetaData& metaData = realm_.meta_data();
 
@@ -275,7 +275,7 @@ TpetraSegregatedLinearSystem::beginLinearSystemConstruction()
     myLIDs_[entityId] = localId++;
     ownedGids.push_back(static_cast<GlobalOrdinal>(entityId));
   }
-  ThrowRequire(localId == numOwnedNodes);
+  STK_ThrowRequire(localId == numOwnedNodes);
 
   // now sharedNotOwned:
   for (const stk::mesh::Bucket* bptr : buckets) {
@@ -336,7 +336,7 @@ TpetraSegregatedLinearSystem::insert_connection(
 {
   const size_t idx = entityToLIDHost_[a.local_offset()];
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     idx < ownedAndSharedNodes_.size(),
     "Error, insert_connection got index out of range.");
 
@@ -351,7 +351,7 @@ TpetraSegregatedLinearSystem::insert_connection(
     correctEntity =
       ownedAndSharedNodes_[idx] == master || naluid_a == naluid_master;
   }
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     correctEntity,
     "Error, indexing of rowEntities to connections isn't right.");
 
@@ -524,7 +524,7 @@ TpetraSegregatedLinearSystem::buildFaceElemToNodeGraph(
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulkData.begin_elements(face);
-      ThrowAssert(bulkData.num_elements(face) == 1);
+      STK_ThrowAssert(bulkData.num_elements(face) == 1);
 
       // get connected element and nodal relations
       stk::mesh::Entity element = face_elem_rels[0];
@@ -642,8 +642,8 @@ TpetraSegregatedLinearSystem::copy_stk_to_tpetra(
   stk::mesh::FieldBase* stkField,
   const Teuchos::RCP<LinSys::MultiVector> tpetraField)
 {
-  ThrowAssert(!tpetraField.is_null());
-  ThrowAssert(stkField);
+  STK_ThrowAssert(!tpetraField.is_null());
+  STK_ThrowAssert(stkField);
   const int numVectors = tpetraField->getNumVectors();
 
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
@@ -663,7 +663,7 @@ TpetraSegregatedLinearSystem::copy_stk_to_tpetra(
     const int fieldSize =
       field_bytes_per_entity(*stkField, b) / (sizeof(double));
 
-    ThrowRequire(numVectors == fieldSize);
+    STK_ThrowRequire(numVectors == fieldSize);
 
     const stk::mesh::Bucket::size_type length = b.size();
 
@@ -966,7 +966,7 @@ TpetraSegregatedLinearSystem::storeOwnersForShared()
 void
 TpetraSegregatedLinearSystem::finalizeLinearSystem()
 {
-  ThrowRequire(inConstruction_);
+  STK_ThrowRequire(inConstruction_);
   inConstruction_ = false;
 
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
@@ -1100,10 +1100,10 @@ TpetraSegregatedLinearSystem::finalizeLinearSystem()
 void
 TpetraSegregatedLinearSystem::zeroSystem()
 {
-  ThrowRequire(!ownedMatrix_.is_null());
-  ThrowRequire(!sharedNotOwnedMatrix_.is_null());
-  ThrowRequire(!sharedNotOwnedRhs_.is_null());
-  ThrowRequire(!ownedRhs_.is_null());
+  STK_ThrowRequire(!ownedMatrix_.is_null());
+  STK_ThrowRequire(!sharedNotOwnedMatrix_.is_null());
+  STK_ThrowRequire(!sharedNotOwnedRhs_.is_null());
+  STK_ThrowRequire(!ownedRhs_.is_null());
 
   sharedNotOwnedMatrix_->resumeFill();
   ownedMatrix_->resumeFill();
@@ -1144,7 +1144,7 @@ segregated_sum_into_row(
     }
 
     if (offset < length) {
-      ThrowAssertMsg(
+      STK_ThrowAssertMsg(
         std::isfinite(input_values[perm_index * numDof]), "Inf or NAN lhs");
       if (forceAtomic) {
         Kokkos::atomic_add(
@@ -1283,7 +1283,7 @@ reset_rows(
       const LocalOrdinal actualLocalId =
         useOwned ? localId : (localId - maxOwnedRowId);
 
-      NGP_ThrowRequire(localId <= maxSharedNotOwnedRowId);
+      STK_NGP_ThrowRequireMsg(localId <= maxSharedNotOwnedRowId, "Error");
 
       // Adjust the LHS; zero out all entries (including diagonal)
       reset_row(localMatrix.row(actualLocalId), actualLocalId, diag_value);
@@ -1371,10 +1371,11 @@ TpetraSegregatedLinearSystem::sumInto(
   const SharedMemView<int*, DeviceShmem>& sortPermutation,
   const char* /* trace_tag */)
 {
-  ThrowAssertMsg(lhs.span_is_contiguous(), "LHS assumed contiguous");
-  ThrowAssertMsg(rhs.span_is_contiguous(), "RHS assumed contiguous");
-  ThrowAssertMsg(localIds.span_is_contiguous(), "localIds assumed contiguous");
-  ThrowAssertMsg(
+  STK_ThrowAssertMsg(lhs.span_is_contiguous(), "LHS assumed contiguous");
+  STK_ThrowAssertMsg(rhs.span_is_contiguous(), "RHS assumed contiguous");
+  STK_ThrowAssertMsg(
+    localIds.span_is_contiguous(), "localIds assumed contiguous");
+  STK_ThrowAssertMsg(
     sortPermutation.span_is_contiguous(), "sortPermutation assumed contiguous");
 
   segregated_sum_into(
@@ -1396,8 +1397,8 @@ TpetraSegregatedLinearSystem::sumInto(
   const size_t n_obj = entities.size();
   const unsigned numRows = n_obj;
 
-  ThrowAssert(numRows * numDof_ == rhs.size());
-  ThrowAssert(numRows * numDof_ * numRows * numDof_ == lhs.size());
+  STK_ThrowAssert(numRows * numDof_ == rhs.size());
+  STK_ThrowAssert(numRows * numDof_ * numRows * numDof_ == lhs.size());
 
   scratchIds.resize(numRows);
   sortPermutation_.resize(numRows);
@@ -1422,7 +1423,7 @@ TpetraSegregatedLinearSystem::sumInto(
 
       for (unsigned dofIdx = 0; dofIdx < numDof_; ++dofIdx) {
         const double cur_rhs = rhs[cur_perm_index * numDof_ + dofIdx];
-        ThrowAssertMsg(std::isfinite(cur_rhs), "Invalid rhs");
+        STK_ThrowAssertMsg(std::isfinite(cur_rhs), "Invalid rhs");
         getOwnedLocalRhs()(rowLid, dofIdx) += cur_rhs;
       }
     } else if (rowLid < maxSharedNotOwnedRowId_) {
@@ -1433,7 +1434,7 @@ TpetraSegregatedLinearSystem::sumInto(
 
       for (unsigned dofIdx = 0; dofIdx < numDof_; ++dofIdx) {
         const double cur_rhs = rhs[cur_perm_index * numDof_ + dofIdx];
-        ThrowAssertMsg(std::isfinite(cur_rhs), "Invalid rhs");
+        STK_ThrowAssertMsg(std::isfinite(cur_rhs), "Invalid rhs");
         getSharedNotOwnedLocalRhs()(actualLocalId, dofIdx) += cur_rhs;
       }
     }
@@ -1464,7 +1465,7 @@ TpetraSegregatedLinearSystem::applyDirichletBCs(
 
     const unsigned fieldSize =
       field_bytes_per_entity(*solutionField, b) / sizeof(double);
-    ThrowRequire(fieldSize == numDof_);
+    STK_ThrowRequire(fieldSize == numDof_);
 
     const stk::mesh::Bucket::size_type length = b.size();
     const double* solution =
@@ -1966,8 +1967,8 @@ TpetraSegregatedLinearSystem::copy_tpetra_to_stk(
   stk::mesh::BulkData& bulkData = realm_.bulk_data();
   stk::mesh::MetaData& metaData = realm_.meta_data();
 
-  ThrowAssert(!tpetraField.is_null());
-  ThrowAssert(stkField);
+  STK_ThrowAssert(!tpetraField.is_null());
+  STK_ThrowAssert(stkField);
   const LinSys::ConstOneDVector& tpetraVector = tpetraField->get1dView();
   const size_t numNodes = tpetraField->getLocalLength();
 
@@ -1986,7 +1987,7 @@ TpetraSegregatedLinearSystem::copy_tpetra_to_stk(
 
     const unsigned fieldSize =
       field_bytes_per_entity(*stkField, b) / sizeof(double);
-    ThrowRequire(fieldSize == numDof_);
+    STK_ThrowRequire(fieldSize == numDof_);
 
     const stk::mesh::Bucket::size_type length = b.size();
     double* stkFieldPtr = (double*)stk::mesh::field_data(*stkField, *b.begin());
@@ -2015,7 +2016,7 @@ TpetraSegregatedLinearSystem::copy_tpetra_to_stk(
                     << " stkId= " << stkId << " naluId= " << naluId
                     << std::endl;
         }
-        ThrowRequire(useOwned);
+        STK_ThrowRequire(useOwned);
 
         const size_t stkIndex = k * numDof_ + dofIdx;
         if (useOwned) {
diff --git a/src/TurbulenceAveragingPostProcessing.C b/src/TurbulenceAveragingPostProcessing.C
index f33872491..2de0bc059 100644
--- a/src/TurbulenceAveragingPostProcessing.C
+++ b/src/TurbulenceAveragingPostProcessing.C
@@ -291,7 +291,7 @@ TurbulenceAveragingPostProcessing::setup()
 
     auto* tempField =
       metaData.get_field(stk::topology::NODE_RANK, "temperature");
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       tempField != nullptr, "Temperature field must be registered");
 
     auto& field =
@@ -522,8 +522,7 @@ TurbulenceAveragingPostProcessing::register_field_from_primitive(
       primitiveName);
 
   // extract size (would love to do this by part), however, not yet a use case
-  const unsigned fieldSizePrimitive =
-    primitiveField->max_size(stk::topology::NODE_RANK);
+  const unsigned fieldSizePrimitive = primitiveField->max_size();
 
   // register the averaged field with this size; treat velocity as a special
   // case to retain the vector aspect
@@ -564,8 +563,7 @@ TurbulenceAveragingPostProcessing::construct_pair(
     metaData.get_field(stk::topology::NODE_RANK, averagedName);
 
   // the size; guaranteed to be the same based on the field registration
-  const unsigned fieldSizeAveraged =
-    averagedField->max_size(stk::topology::NODE_RANK);
+  const unsigned fieldSizeAveraged = averagedField->max_size();
   fieldSizeVec.push_back(fieldSizeAveraged);
 
   // construct pairs
diff --git a/src/WilcoxKOmegaEquationSystem.C b/src/WilcoxKOmegaEquationSystem.C
index 32a116544..98df0d228 100644
--- a/src/WilcoxKOmegaEquationSystem.C
+++ b/src/WilcoxKOmegaEquationSystem.C
@@ -264,7 +264,7 @@ WilcoxKOmegaEquationSystem::post_external_data_transfer_work()
   auto sdrBCField = meta.get_field<double>(stk::topology::NODE_RANK, "sdr_bc");
   auto tkeBCField = meta.get_field<double>(stk::topology::NODE_RANK, "tke_bc");
   if (sdrBCField != nullptr) {
-    ThrowRequire(tkeBCField);
+    STK_ThrowRequire(tkeBCField);
     auto bc_sel = owned_and_shared & stk::mesh::selectField(*sdrBCField);
     auto ngpTkeBC =
       fieldMgr.get_field<double>(tkeBCField->mesh_meta_data_ordinal());
diff --git a/src/aero/actuator/ActuatorBladeDistributor.C b/src/aero/actuator/ActuatorBladeDistributor.C
index faf9ea005..dbc691542 100644
--- a/src/aero/actuator/ActuatorBladeDistributor.C
+++ b/src/aero/actuator/ActuatorBladeDistributor.C
@@ -26,7 +26,7 @@ blade_belongs_on_this_rank(
   const bool isInDivisionIncrement =
     globBladeNum >= (div * rank) && globBladeNum < div * (rank + 1);
   const bool isInRemainderIncrement = (globBladeNum - div * numRanks) == rank;
-  ThrowAssert(!(isInDivisionIncrement && isInRemainderIncrement));
+  STK_ThrowAssert(!(isInDivisionIncrement && isInRemainderIncrement));
   return isInDivisionIncrement || isInRemainderIncrement;
 }
 
diff --git a/src/aero/actuator/ActuatorBulkDiskFAST.C b/src/aero/actuator/ActuatorBulkDiskFAST.C
index 074d3003f..34a114bb3 100644
--- a/src/aero/actuator/ActuatorBulkDiskFAST.C
+++ b/src/aero/actuator/ActuatorBulkDiskFAST.C
@@ -25,7 +25,7 @@ ActuatorBulkDiskFAST::ActuatorBulkDiskFAST(
       "numSweptOffset", actMeta.numberOfActuators_, actMeta.maxNumPntsPerBlade_)
 {
 
-  ThrowErrorIf(!actMeta.is_disk());
+  STK_ThrowErrorIf(!actMeta.is_disk());
   compute_swept_point_count(actMeta);
   resize_arrays(actMeta);
   Kokkos::parallel_for(
diff --git a/src/aero/actuator/ActuatorBulkFAST.C b/src/aero/actuator/ActuatorBulkFAST.C
index b98755787..02c2915db 100644
--- a/src/aero/actuator/ActuatorBulkFAST.C
+++ b/src/aero/actuator/ActuatorBulkFAST.C
@@ -89,7 +89,7 @@ ActuatorBulkFAST::init_openfast(
   const int remainder = actMeta.numberOfActuators_ % nProcs;
   const int nOffset = intDivision * nProcs;
 
-  ThrowErrorMsgIf(
+  STK_ThrowErrorMsgIf(
     remainder && intDivision,
     "nalu-wind can't process more turbines than ranks.");
 
@@ -120,7 +120,7 @@ ActuatorBulkFAST::init_openfast(
 
   for (int i = 0; i < nTurb; ++i) {
     if (localTurbineId_ == openFast_.get_procNo(i)) {
-      ThrowErrorMsgIf(
+      STK_ThrowErrorMsgIf(
         actMeta.nBlades_(i) != openFast_.get_numBlades(i),
         "Mismatch in number of blades between OpenFAST and input deck."
         " InputDeck: " +
@@ -146,7 +146,8 @@ ActuatorBulkFAST::init_epsilon(const ActuatorMetaFAST& actMeta)
 
   for (int iTurb = 0; iTurb < nTurb; iTurb++) {
     if (openFast_.get_procNo(iTurb) == NaluEnv::self().parallel_rank()) {
-      ThrowAssert(actMeta.numPointsTotal_ >= openFast_.get_numForcePts(iTurb));
+      STK_ThrowAssert(
+        actMeta.numPointsTotal_ >= openFast_.get_numForcePts(iTurb));
       const int numForcePts = actMeta.numPointsTurbine_.h_view(iTurb);
       const int offset = turbIdOffset_.h_view(iTurb);
       auto epsilonChord =
@@ -215,7 +216,7 @@ ActuatorBulkFAST::init_epsilon(const ActuatorMetaFAST& actMeta)
         }
 
         for (int i = 0; i < 3; ++i) {
-          ThrowAssertMsg(
+          STK_ThrowAssertMsg(
             epsilonLocal(i) > 0.0,
             "Epsilon zero for point: " + std::to_string(np) + " index " +
               std::to_string(i));
diff --git a/src/aero/actuator/ActuatorFLLC.C b/src/aero/actuator/ActuatorFLLC.C
index edd22a214..cdf7b4630 100644
--- a/src/aero/actuator/ActuatorFLLC.C
+++ b/src/aero/actuator/ActuatorFLLC.C
@@ -56,7 +56,7 @@ FilteredLiftingLineCorrection::compute_lift_force_distribution()
     const auto offset = info.offset_;
     const auto nPoints = info.nPoints_;
     // debug size checks
-    ThrowAssert(offset + nPoints < static_cast<int>(G.size()));
+    STK_ThrowAssert(offset + nPoints < static_cast<int>(G.size()));
 
     auto range_policy =
       Kokkos::RangePolicy<exec_space>(offset, offset + nPoints);
@@ -97,8 +97,8 @@ FilteredLiftingLineCorrection::grad_lift_force_distribution()
 
     const auto offset = info.offset_;
     const auto nPoints = info.nPoints_;
-    ThrowAssert(offset + nPoints < static_cast<int>(G.size()));
-    ThrowAssert(offset + nPoints < static_cast<int>(deltaG.size()));
+    STK_ThrowAssert(offset + nPoints < static_cast<int>(G.size()));
+    STK_ThrowAssert(offset + nPoints < static_cast<int>(deltaG.size()));
 
     auto range_policy =
       Kokkos::RangePolicy<exec_space>(offset, offset + nPoints);
diff --git a/src/aero/actuator/ActuatorFunctorsFAST.C b/src/aero/actuator/ActuatorFunctorsFAST.C
index 0d4cd10e9..00d486138 100644
--- a/src/aero/actuator/ActuatorFunctorsFAST.C
+++ b/src/aero/actuator/ActuatorFunctorsFAST.C
@@ -54,7 +54,7 @@ void
 ActFastUpdatePoints::operator()(int index) const
 {
 
-  ThrowAssert(turbId_ >= 0);
+  STK_ThrowAssert(turbId_ >= 0);
   const int pointId = index - offsets_(turbId_);
   auto point = Kokkos::subview(points_, index, Kokkos::ALL);
 
diff --git a/src/aero/actuator/ActuatorFunctorsSimple.C b/src/aero/actuator/ActuatorFunctorsSimple.C
index 0b9835412..b87b5a2e8 100644
--- a/src/aero/actuator/ActuatorFunctorsSimple.C
+++ b/src/aero/actuator/ActuatorFunctorsSimple.C
@@ -52,7 +52,7 @@ InterpActuatorDensity::operator()(int index) const
     double ws_coordinates[max_size], ws_density[max_size];
 
     // Check to make sure the size is sufficient
-    ThrowAssert(max_size >= 3 * nodesPerElem);
+    STK_ThrowAssert(max_size >= 3 * nodesPerElem);
 
     actuator_utils::gather_field(
       3, &ws_coordinates[0], *coordinates_, stkBulk_.begin_nodes(elem),
@@ -87,7 +87,7 @@ void
 ActSimpleUpdatePoints::operator()(int index) const
 {
 
-  ThrowAssert(turbId_ >= 0);
+  STK_ThrowAssert(turbId_ >= 0);
   const int pointId = index - offsets_(turbId_);
   auto point = Kokkos::subview(points_, index, Kokkos::ALL);
 
@@ -118,7 +118,7 @@ ActSimpleWriteToFile(
 
   if (actBulk.localTurbineId_ == NaluEnv::self().parallel_rank()) {
     std::ofstream outFile;
-    // ThrowErrorIf(NaluEnv::self().parallel_rank()!=0);
+    // STK_ThrowErrorIf(NaluEnv::self().parallel_rank()!=0);
 
     outFile.open(filename, std::ios_base::app);
     const int stop =
diff --git a/src/aero/actuator/ActuatorModel.C b/src/aero/actuator/ActuatorModel.C
index 106e0c1be..056b26d1e 100644
--- a/src/aero/actuator/ActuatorModel.C
+++ b/src/aero/actuator/ActuatorModel.C
@@ -71,12 +71,12 @@ ActuatorModel::setup(double timeStep, stk::mesh::BulkData& stkBulk)
     return;
 
   // hack to surpress Wunused-parameter on non-openfast builds
-  ThrowErrorIf(timeStep <= 0.0);
+  STK_ThrowErrorIf(timeStep <= 0.0);
 
   switch (actMeta_->actuatorType_) {
   case (ActuatorType::ActLineFASTNGP): {
 #ifndef NALU_USES_OPENFAST
-    ThrowErrorMsg("Actuator methods require OpenFAST");
+    STK_ThrowErrorMsg("Actuator methods require OpenFAST");
 #if !defined(KOKKOS_ENABLE_GPU)
     break;
 #endif
@@ -94,7 +94,7 @@ ActuatorModel::setup(double timeStep, stk::mesh::BulkData& stkBulk)
   }
   case (ActuatorType::ActDiskFASTNGP): {
 #ifndef NALU_USES_OPENFAST
-    ThrowErrorMsg("Actuator methods require OpenFAST");
+    STK_ThrowErrorMsg("Actuator methods require OpenFAST");
 #if !defined(KOKKOS_ENABLE_GPU)
     break;
 #endif
@@ -122,7 +122,7 @@ ActuatorModel::setup(double timeStep, stk::mesh::BulkData& stkBulk)
     break;
   }
   default: {
-    ThrowErrorMsg("Unsupported actuator type");
+    STK_ThrowErrorMsg("Unsupported actuator type");
   }
   }
 }
@@ -140,7 +140,7 @@ ActuatorModel::init(stk::mesh::BulkData& stkBulk)
   case (ActuatorType::ActLineFASTNGP):
   case (ActuatorType::ActDiskFASTNGP): {
 #ifndef NALU_USES_OPENFAST
-    ThrowErrorMsg("Actuator methods require OpenFAST");
+    STK_ThrowErrorMsg("Actuator methods require OpenFAST");
 #if !defined(KOKKOS_ENABLE_GPU)
     break;
 #endif
@@ -154,7 +154,7 @@ ActuatorModel::init(stk::mesh::BulkData& stkBulk)
     break;
   }
   default: {
-    ThrowErrorMsg("Unsupported actuator type");
+    STK_ThrowErrorMsg("Unsupported actuator type");
   }
   }
 }
diff --git a/src/aero/actuator/ActuatorParsing.C b/src/aero/actuator/ActuatorParsing.C
index 9878c3bd5..b6f8d8856 100644
--- a/src/aero/actuator/ActuatorParsing.C
+++ b/src/aero/actuator/ActuatorParsing.C
@@ -98,7 +98,7 @@ ActuatorMeta
 actuator_parse(const YAML::Node& y_node)
 {
   const YAML::Node y_actuator = y_node["actuator"];
-  ThrowErrorMsgIf(
+  STK_ThrowErrorMsgIf(
     !y_actuator, "actuator argument is "
                  "missing from yaml node passed to actuator_parse");
   int nTurbines = 0;
diff --git a/src/aero/actuator/ActuatorParsingFAST.C b/src/aero/actuator/ActuatorParsingFAST.C
index 00c103fc1..e83016cd5 100644
--- a/src/aero/actuator/ActuatorParsingFAST.C
+++ b/src/aero/actuator/ActuatorParsingFAST.C
@@ -110,7 +110,7 @@ readTurbineData(int iTurb, ActuatorMetaFAST& actMetaFAST, YAML::Node turbNode)
   int* numBlades = &(actMetaFAST.nBlades_(iTurb));
   *numBlades = 3;
   get_if_present_no_default(turbNode, "num_blades", *numBlades);
-  ThrowErrorMsgIf(
+  STK_ThrowErrorMsgIf(
     *numBlades != 3 && *numBlades != 2,
     "ERROR::ActuatorParsingFAST::Currently only 2 and 3 bladed turbines are "
     "supported.");
@@ -120,7 +120,7 @@ readTurbineData(int iTurb, ActuatorMetaFAST& actMetaFAST, YAML::Node turbNode)
       turbNode, "num_swept_pts", actMetaFAST.nPointsSwept_(iTurb));
     actMetaFAST.useUniformAziSampling_(iTurb) =
       actMetaFAST.nPointsSwept_(iTurb) != 0;
-    ThrowErrorMsgIf(
+    STK_ThrowErrorMsgIf(
       *numBlades != 3, "The ActuatorDisk model requires a base 3 bladed "
                        "turbine, but a 2 bladed turbine was supplied.");
   }
@@ -142,7 +142,7 @@ actuator_FAST_parse(const YAML::Node& y_node, const ActuatorMeta& actMeta)
   fi.nTurbinesGlob = actMetaFAST.numberOfActuators_;
 
   const YAML::Node y_actuator = y_node["actuator"];
-  ThrowErrorMsgIf(
+  STK_ThrowErrorMsgIf(
     !y_actuator, "actuator argument is "
                  "missing from yaml node passed to actuator_FAST_parse");
   if (fi.nTurbinesGlob > 0) {
diff --git a/src/aero/actuator/ActuatorParsingSimple.C b/src/aero/actuator/ActuatorParsingSimple.C
index a980733ba..a486d9858 100644
--- a/src/aero/actuator/ActuatorParsingSimple.C
+++ b/src/aero/actuator/ActuatorParsingSimple.C
@@ -55,7 +55,7 @@ actuator_Simple_parse(const YAML::Node& y_node, const ActuatorMeta& actMeta)
     << "In actuator_Simple_parse() " << std::endl; // LCCOUT
 
   const YAML::Node y_actuator = y_node["actuator"];
-  ThrowErrorMsgIf(
+  STK_ThrowErrorMsgIf(
     !y_actuator, "actuator argument is "
                  "missing from yaml node passed to actuator_Simple_parse");
 
diff --git a/src/aero/actuator/ActuatorSearch.C b/src/aero/actuator/ActuatorSearch.C
index 168f4aee1..360bfc5f5 100644
--- a/src/aero/actuator/ActuatorSearch.C
+++ b/src/aero/actuator/ActuatorSearch.C
@@ -161,8 +161,8 @@ ExecuteFineSearch(
 {
   const int nDim = 3;
 
-  ThrowAssert(isLocalPoint.extent(0) == points.extent(0));
-  ThrowAssert(coarsePointIds.extent(0) == coarseElemIds.extent(0));
+  STK_ThrowAssert(isLocalPoint.extent(0) == points.extent(0));
+  STK_ThrowAssert(coarsePointIds.extent(0) == coarseElemIds.extent(0));
 
   // extract fields
   stk::mesh::MetaData& stkMeta = stkBulk.mesh_meta_data();
diff --git a/src/aero/actuator/UtilitiesActuator.C b/src/aero/actuator/UtilitiesActuator.C
index e7fe0c235..c6dbf23b1 100644
--- a/src/aero/actuator/UtilitiesActuator.C
+++ b/src/aero/actuator/UtilitiesActuator.C
@@ -85,7 +85,7 @@ get_fast_point_index(
     break;
   }
   default: {
-    ThrowErrorMsg("Invalid fast type");
+    STK_ThrowErrorMsg("Invalid fast type");
     return -1;
     break;
   }
diff --git a/src/aero/fsi/CalcLoads.C b/src/aero/fsi/CalcLoads.C
index 0119860d3..3af7f9833 100644
--- a/src/aero/fsi/CalcLoads.C
+++ b/src/aero/fsi/CalcLoads.C
@@ -142,7 +142,7 @@ CalcLoads::execute()
 
     // extract connected element topology
     b->parent_topology(stk::topology::ELEMENT_RANK, parentTopo);
-    ThrowAssert(parentTopo.size() == 1);
+    STK_ThrowAssert(parentTopo.size() == 1);
     stk::topology theElemTopo = parentTopo[0];
 
     // extract master element for this element topo
@@ -197,7 +197,7 @@ CalcLoads::execute()
       // extract the connected element to this exposed face; should be single in
       // size!
       const stk::mesh::Entity* face_elem_rels = bulk_->begin_elements(face);
-      ThrowAssert(bulk_->num_elements(face) == 1);
+      STK_ThrowAssert(bulk_->num_elements(face) == 1);
 
       // get element; its face ordinal number
       stk::mesh::Entity element = face_elem_rels[0];
diff --git a/src/aero/fsi/FSIturbine.C b/src/aero/fsi/FSIturbine.C
index 72b27de7f..46ce95454 100644
--- a/src/aero/fsi/FSIturbine.C
+++ b/src/aero/fsi/FSIturbine.C
@@ -88,7 +88,7 @@ fsiTurbine::fsiTurbine(int iTurb, const YAML::Node& node)
     // Displacement Ramping
     // --------------------------------------------------------------------------
     const YAML::Node defNode = node["deflection_ramping"];
-    ThrowErrorMsgIf(
+    STK_ThrowErrorMsgIf(
       !defNode,
       "defleciton_ramping inputs are required for FSI Turbines with blades");
     DeflectionRampingParams& defParams = deflectionRampParams_;
@@ -653,7 +653,7 @@ fsiTurbine::write_nc_def_loads(const size_t tStep, const double curTime)
         (0.5 * (brFSIdata_.bld_rloc[i + iStart + 1] -
                 brFSIdata_.bld_rloc[i + iStart - 1]));
     }
-    brFSIdata_.bld_ld[(iStart)*6 + 4] =
+    brFSIdata_.bld_ld[(iStart) * 6 + 4] =
       (0.5 * (brFSIdata_.bld_rloc[iStart + 1] - brFSIdata_.bld_rloc[iStart]));
     brFSIdata_.bld_ld[(iStart + nBldPts - 1) * 6 + 4] =
       (0.5 * (brFSIdata_.bld_rloc[iStart + nBldPts - 1] -
@@ -1654,7 +1654,8 @@ fsiTurbine::computeMapping()
 
   auto& meta = bulk_->mesh_meta_data();
   const int ndim = meta.spatial_dimension();
-  ThrowRequireMsg(ndim == 3, "fsiTurbine: spatial dim is required to be 3.");
+  STK_ThrowRequireMsg(
+    ndim == 3, "fsiTurbine: spatial dim is required to be 3.");
   const VectorFieldType* modelCoords =
     meta.get_field<double>(stk::topology::NODE_RANK, "coordinates");
   modelCoords->sync_to_host();
diff --git a/src/aero/fsi/MapLoad.C b/src/aero/fsi/MapLoad.C
index 8d11fd14a..acf320194 100644
--- a/src/aero/fsi/MapLoad.C
+++ b/src/aero/fsi/MapLoad.C
@@ -110,11 +110,11 @@ mapTowerLoad(
 
         // Find the interpolated reference position first
         linInterpVec(
-          &twrRefPos[(loadMap_bip)*6], &twrRefPos[(loadMap_bip + 1) * 6],
+          &twrRefPos[(loadMap_bip) * 6], &twrRefPos[(loadMap_bip + 1) * 6],
           loadMapInterp_bip, tmpNodePos.data());
         // Find the interpolated linear displacement
         linInterpVec(
-          &twrDef[(loadMap_bip)*6], &twrDef[(loadMap_bip + 1) * 6],
+          &twrDef[(loadMap_bip) * 6], &twrDef[(loadMap_bip + 1) * 6],
           loadMapInterp_bip, tmpNodeDisp.data());
         // Add displacement to find actual position
         for (auto idim = 0; idim < 3; idim++)
@@ -130,8 +130,8 @@ mapTowerLoad(
         // Split the force and moment into the two surrounding nodes in a
         // variationally consistent manner using the interpolation factor
         fsi::splitForceMoment(
-          tmpForceMoment.data(), loadMapInterp_bip, &(twrLoad[(loadMap_bip)*6]),
-          &(twrLoad[(loadMap_bip + 1) * 6]));
+          tmpForceMoment.data(), loadMapInterp_bip,
+          &(twrLoad[(loadMap_bip) * 6]), &(twrLoad[(loadMap_bip + 1) * 6]));
       }
     }
   }
diff --git a/src/aero/fsi/OpenfastFSI.C b/src/aero/fsi/OpenfastFSI.C
index 38ed14a69..1b53996ba 100644
--- a/src/aero/fsi/OpenfastFSI.C
+++ b/src/aero/fsi/OpenfastFSI.C
@@ -450,8 +450,8 @@ OpenfastFSI::send_loads(const double curTime)
           (fsiTurbineData_[i]->params_.nBRfsiPtsTwr) * 6, MPI_DOUBLE, MPI_SUM,
           turbProc, bulk_->parallel());
         iError = MPI_Reduce(
-          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL, (nTotBldNodes)*6,
-          MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
+          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL,
+          (nTotBldNodes) * 6, MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
       }
     }
   }
@@ -664,8 +664,8 @@ OpenfastFSI::map_loads(const int tStep, const double curTime)
           (fsiTurbineData_[i]->params_.nBRfsiPtsTwr) * 6, MPI_DOUBLE, MPI_SUM,
           turbProc, bulk_->parallel());
         iError = MPI_Reduce(
-          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL, (nTotBldNodes)*6,
-          MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
+          fsiTurbineData_[i]->brFSIdata_.bld_ld.data(), NULL,
+          (nTotBldNodes) * 6, MPI_DOUBLE, MPI_SUM, turbProc, bulk_->parallel());
       }
 
       fsiTurbineData_[i]->write_nc_def_loads(tStep, curTime);
diff --git a/src/edge_kernels/MomentumOpenEdgeKernel.C b/src/edge_kernels/MomentumOpenEdgeKernel.C
index e0627eb08..8973b2f16 100644
--- a/src/edge_kernels/MomentumOpenEdgeKernel.C
+++ b/src/edge_kernels/MomentumOpenEdgeKernel.C
@@ -268,7 +268,7 @@ MomentumOpenEdgeKernel<BcAlgTraits>::execute(
       break;
     }
     default:
-      NGP_ThrowErrorMsg("invalid entrainment method");
+      STK_NGP_ThrowErrorMsg("invalid entrainment method");
     }
   }
 }
diff --git a/src/element_promotion/PromoteElement.C b/src/element_promotion/PromoteElement.C
index 2baf2bc02..283284a05 100644
--- a/src/element_promotion/PromoteElement.C
+++ b/src/element_promotion/PromoteElement.C
@@ -50,7 +50,7 @@ create_tensor_product_hex_elements(
   const VectorFieldType& coordField,
   const stk::mesh::PartVector& partsToBePromoted)
 {
-  ThrowRequire(check_parts_for_promotion(partsToBePromoted));
+  STK_ThrowRequire(check_parts_for_promotion(partsToBePromoted));
   return impl::promote_elements_hex(
     nodeLocs1D, bulk, coordField, partsToBePromoted);
 }
diff --git a/src/element_promotion/PromoteElementImpl.C b/src/element_promotion/PromoteElementImpl.C
index dd59cd856..e4d8a7d65 100644
--- a/src/element_promotion/PromoteElementImpl.C
+++ b/src/element_promotion/PromoteElementImpl.C
@@ -181,7 +181,7 @@ create_boundary_elements(
         subset->topology().rank() == side_rank &&
         !subset->topology().is_super_topology()) {
         soloFacePart[0] = super_subset_part(*subset);
-        ThrowRequire(soloFacePart[0] != nullptr);
+        STK_ThrowRequire(soloFacePart[0] != nullptr);
 
         const auto& buckets = bulk.get_buckets(side_rank, *subset);
         bucket_loop(buckets, [&](const stk::mesh::Entity side) {
@@ -189,7 +189,7 @@ create_boundary_elements(
             bulk.declare_solo_side(availableFaceIds[faceIdIndex], soloFacePart);
           stk::mesh::Entity superElem = sideToSuperElemMap.at(side);
 
-          ThrowRequireMsg(
+          STK_ThrowRequireMsg(
             bulk.num_elements(side) == 1u,
             "Multiple elements attached to boundary side");
           const auto sideOrdinal = bulk.begin_element_ordinals(side)[0];
@@ -245,7 +245,7 @@ perform_parallel_consolidation_of_node_ids(
     for (const auto& pair : connectivityMap) {
       auto entKey = bulk.entity_key(pair.first);
       const auto& nodeIds = pair.second;
-      ThrowRequire(entKey.rank() == domainTopoRank);
+      STK_ThrowRequire(entKey.rank() == domainTopoRank);
 
       std::vector<int> procs;
       bulk.comm_shared_procs(entKey, procs);
@@ -407,7 +407,7 @@ index_face_nodes(int i, int j, int len1D, stk::mesh::Permutation perm)
     break;
   }
   default: {
-    ThrowRequireMsg(false, "Invalid permutation of quad face");
+    STK_ThrowRequireMsg(false, "Invalid permutation of quad face");
   }
   }
 
@@ -433,8 +433,8 @@ add_edge_nodes_to_elem_connectivity(
       edgeConnectivity.at(edge_rels[edge_ord]);
     const std::vector<int>& ords = desc.edge_node_connectivities(edge_ord);
 
-    ThrowAssert(nodeIds.size() == ords.size());
-    ThrowAssert(static_cast<int>(ords.size()) == newNodesPerEdge);
+    STK_ThrowAssert(nodeIds.size() == ords.size());
+    STK_ThrowAssert(static_cast<int>(ords.size()) == newNodesPerEdge);
     for (int i = 0; i < newNodesPerEdge; ++i) {
       allNodes.at(ords.at(i)) =
         nodeIds.at(index_edge_nodes(i, newNodesPerEdge, perm[edge_ord]));
@@ -461,9 +461,9 @@ add_face_nodes_to_elem_connectivity(
       faceConnectivity.at(face_rels[face_ord]);
     const std::vector<int>& ords = desc.face_node_connectivities(face_index);
 
-    ThrowAssert(nodeIds.size() == ords.size());
-    ThrowAssert(desc.newNodesPerFace == static_cast<int>(ords.size()));
-    ThrowAssert(desc.newNodesPerFace == newNodesPerEdge * newNodesPerEdge);
+    STK_ThrowAssert(nodeIds.size() == ords.size());
+    STK_ThrowAssert(desc.newNodesPerFace == static_cast<int>(ords.size()));
+    STK_ThrowAssert(desc.newNodesPerFace == newNodesPerEdge * newNodesPerEdge);
 
     for (int j = 0; j < newNodesPerEdge; ++j) {
       for (int i = 0; i < newNodesPerEdge; ++i) {
@@ -505,7 +505,7 @@ destroy_entity(stk::mesh::BulkData& bulk, stk::mesh::Entity entity)
     for (unsigned irel = 0; irel < relatives.size(); ++irel) {
       bool del =
         bulk.destroy_relation(relatives[irel], entity, relative_ordinals[irel]);
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         del, "Failed to disconnect entity: " +
                std::to_string(bulk.identifier(entity)));
     }
@@ -524,10 +524,10 @@ destroy_entities(
     selector, bulk.get_buckets(rank, selector), entities);
 
   for (stk::mesh::Entity entity : entities) {
-    ThrowRequire(bulk.is_valid(entity));
+    STK_ThrowRequire(bulk.is_valid(entity));
     if (bulk.bucket(entity).owned()) {
       bool destroyed = destroy_entity(bulk, entity);
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         destroyed,
         "Failed to destroy entity: " + std::to_string(bulk.identifier(entity)));
     }
@@ -585,7 +585,8 @@ set_coordinates_hex(
   std::array<double, 3> physCoords;
 
   bucket_loop(elem_buckets, [&](const stk::mesh::Entity elem) {
-    ThrowAssert(desc.nodesPerElement == static_cast<int>(bulk.num_nodes(elem)));
+    STK_ThrowAssert(
+      desc.nodesPerElement == static_cast<int>(bulk.num_nodes(elem)));
     const stk::mesh::Entity* node_rels = bulk.begin_nodes(elem);
 
     for (int ord = 0; ord < 8; ++ord) {
@@ -626,13 +627,13 @@ make_base_nodes_to_elem_map_at_boundary(
   NodesElemMap nodesToElemMap;
   stk::mesh::EntityIdVector parents(baseNumNodes);
   bucket_loop(baseElemSideBuckets, [&](stk::mesh::Entity side) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       mesh.num_elements(side) == 1u,
       "Multiple elements attached to boundary side");
     const stk::mesh::Entity parent_elem = mesh.begin_elements(side)[0];
 
     const auto* node_rels = mesh.begin_nodes(parent_elem);
-    ThrowAssert(mesh.num_nodes(parent_elem) == parents.size());
+    STK_ThrowAssert(mesh.num_nodes(parent_elem) == parents.size());
 
     for (int j = 0; j < desc.nodesInBaseElement; ++j) {
       parents.at(j) = mesh.identifier(node_rels[j]);
@@ -653,7 +654,7 @@ exposed_side_to_super_elem_map(
    * Generates a map between each exposed face and the super-element
    * notionally attached to that exposed face.
    */
-  ThrowRequire(part_vector_is_valid_and_nonempty(
+  STK_ThrowRequire(part_vector_is_valid_and_nonempty(
     super_elem_part_vector(base_elem_mesh_parts)));
 
   const auto& superElemBuckets = bulk.get_buckets(
@@ -670,7 +671,7 @@ exposed_side_to_super_elem_map(
 
   bucket_loop(superElemBuckets, [&](stk::mesh::Entity superElem) {
     const auto* node_rels = bulk.begin_nodes(superElem);
-    ThrowAssert(
+    STK_ThrowAssert(
       static_cast<int>(bulk.num_nodes(superElem)) > baseNumNodes ||
       desc.polyOrder == 1);
 
@@ -685,7 +686,7 @@ exposed_side_to_super_elem_map(
     if (it != nodesToElemMap.end()) {
       const stk::mesh::Entity baseElem = it->second;
       auto result = elemToSuperElemMap.insert({baseElem, superElem});
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         result.second, "Multiple superElems with same parent nodes as the base "
                        "elements found");
     }
@@ -702,11 +703,11 @@ exposed_side_to_super_elem_map(
 
   bucket_loop(boundary_buckets, [&](stk::mesh::Entity side) {
     ;
-    ThrowAssert(bulk.num_elements(side) == 1u);
+    STK_ThrowAssert(bulk.num_elements(side) == 1u);
     const stk::mesh::Entity baseElem = bulk.begin_elements(side)[0];
     const stk::mesh::Entity superElem = elemToSuperElemMap.at(baseElem);
     auto result = exposedSideToSuperElemMap.insert({side, superElem});
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       result.second, "Multiple super elements associated with the same face");
   });
 
diff --git a/src/element_promotion/PromotedElementIO.C b/src/element_promotion/PromotedElementIO.C
index 301ea039a..c108d5b54 100644
--- a/src/element_promotion/PromotedElementIO.C
+++ b/src/element_promotion/PromotedElementIO.C
@@ -49,7 +49,7 @@ size_t
 num_sub_elements(
   const int dim, const stk::mesh::BucketVector& buckets, unsigned polyOrder)
 {
-  ThrowRequire(dim == 2 || dim == 3);
+  STK_ThrowRequire(dim == 2 || dim == 3);
   int faceSize = std::pow(polyOrder, dim - 1);
   int elemSize = std::pow(polyOrder, dim);
   auto side_rank =
@@ -90,11 +90,11 @@ PromotedElementIO::PromotedElementIO(
 
   databaseIO = Ioss::IOFactory::create(
     "exodus", fileName_, Ioss::WRITE_RESULTS, bulkData_.parallel(), properties);
-  ThrowRequire(databaseIO != nullptr && databaseIO->ok(true));
+  STK_ThrowRequire(databaseIO != nullptr && databaseIO->ok(true));
 
   output_ = std::make_unique<Ioss::Region>(
     databaseIO, "HighOrderOutput"); // sink for databaseIO
-  ThrowRequire(output_ != nullptr);
+  STK_ThrowRequire(output_ != nullptr);
 
   const stk::mesh::BucketVector& elem_buckets = bulkData_.get_buckets(
     stk::topology::ELEM_RANK, stk::mesh::selectUnion(baseParts));
@@ -103,10 +103,10 @@ PromotedElementIO::PromotedElementIO(
   std::vector<stk::mesh::EntityId> subElemIds;
 
   bulkData.generate_new_ids(stk::topology::ELEM_RANK, numSubElems, subElemIds);
-  ThrowRequire(subElemIds.size() == numSubElems);
+  STK_ThrowRequire(subElemIds.size() == numSubElems);
 
   superElemParts_ = super_elem_part_vector(baseParts);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     part_vector_is_valid_and_nonempty(superElemParts_),
     "Not all element parts have a super-element mirror");
 
@@ -135,7 +135,7 @@ PromotedElementIO::write_database_data(double currentTime)
   nodeBlock_->get_field_data("ids", ids);
 
   for (const auto& pair : fields_) {
-    ThrowRequire(pair.second != nullptr);
+    STK_ThrowRequire(pair.second != nullptr);
     const stk::mesh::FieldBase& field = *pair.second;
     if (field.type_is<int>()) {
       put_data_on_node_block<int32_t>(*nodeBlock_, ids, field, nodeBuckets);
@@ -148,7 +148,7 @@ PromotedElementIO::write_database_data(double currentTime)
     } else if (field.type_is<double>()) {
       put_data_on_node_block<double>(*nodeBlock_, ids, field, nodeBuckets);
     } else {
-      ThrowRequireMsg(false, "Unsupported type for output");
+      STK_ThrowRequireMsg(false, "Unsupported type for output");
     }
   }
 
@@ -164,8 +164,8 @@ PromotedElementIO::put_data_on_node_block(
   const stk::mesh::FieldBase& field,
   const stk::mesh::BucketVector& buckets) const
 {
-  ThrowRequire(field.type_is<T>());
-  int fieldLength = field.max_size(stk::topology::NODE_RANK);
+  STK_ThrowRequire(field.type_is<T>());
+  int fieldLength = field.max_size();
   std::vector<T> flatArray(count_entities(buckets) * fieldLength);
 
   size_t index = 0;
@@ -198,10 +198,10 @@ PromotedElementIO::write_elem_block_definitions(
       auto block = std::make_unique<Ioss::ElementBlock>(
         databaseIO, baseElemPart->name(), baseElemPart->topology().name(),
         numSubElems);
-      ThrowRequireMsg(block != nullptr, "Element block creation failed");
+      STK_ThrowRequireMsg(block != nullptr, "Element block creation failed");
 
       auto result = elementBlockPointers_.insert({ip, block.get()});
-      ThrowRequireMsg(result.second, "Attempted to add redundant part");
+      STK_ThrowRequireMsg(result.second, "Attempted to add redundant part");
 
       output_->add(block.release());
     }
@@ -217,7 +217,7 @@ PromotedElementIO::write_node_block_definitions(
   auto nodeCount = count_entities(nodeBuckets);
   auto nodeBlock = std::make_unique<Ioss::NodeBlock>(
     databaseIO, "nodeblock", nodeCount, nDim_);
-  ThrowRequireMsg(nodeBlock != nullptr, "Node block creation failed");
+  STK_ThrowRequireMsg(nodeBlock != nullptr, "Node block creation failed");
   nodeBlock_ = nodeBlock.get();
   output_->add(nodeBlock.release());
 }
@@ -328,7 +328,7 @@ PromotedElementIO::add_fields(const std::vector<stk::mesh::FieldBase*>& fields)
       } else if (field.type_is<uint64_t>() || field.type_is<int64_t>()) {
         iossType = Ioss::Field::INT64;
       } else {
-        ThrowRequireMsg(
+        STK_ThrowRequireMsg(
           field.type_is<double>(),
           "Only (u)int32, (u)int64, and double fields supported");
       }
@@ -345,7 +345,7 @@ std::string
 PromotedElementIO::storage_name(const stk::mesh::FieldBase& field) const
 {
   std::string storageType;
-  switch (field.max_size(stk::topology::NODE_RANK)) {
+  switch (field.max_size()) {
   case 1: {
     storageType = "scalar";
     break;
diff --git a/src/element_promotion/PromotedPartHelper.C b/src/element_promotion/PromotedPartHelper.C
index 3643144c8..daf8762ee 100644
--- a/src/element_promotion/PromotedPartHelper.C
+++ b/src/element_promotion/PromotedPartHelper.C
@@ -63,11 +63,11 @@ bool
 check_parts_for_promotion(const stk::mesh::PartVector& parts)
 {
   for (const auto* ip : parts) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       ip != nullptr, "An invalid part was designated for promotion.");
 
     if (ip->topology().rank() == stk::topology::ELEM_RANK) {
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         super_elem_part(*ip) != nullptr,
         "No super element mirror for part requesting promotion");
     }
@@ -78,10 +78,10 @@ check_parts_for_promotion(const stk::mesh::PartVector& parts)
       }
     } else {
       for (const auto* is : part.subsets()) {
-        ThrowRequireMsg(
+        STK_ThrowRequireMsg(
           is != nullptr,
           "An invalid subset part was designated for promotion.");
-        ThrowRequireMsg(
+        STK_ThrowRequireMsg(
           super_subset_part(*is) != nullptr,
           "A subset part lacks its super side mirror.");
         const stk::mesh::Part& subsetPart = *is;
@@ -103,7 +103,7 @@ super_element_suffix()
 std::string
 super_element_part_name(std::string base_name)
 {
-  ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem part");
+  STK_ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem part");
   return (base_name + super_element_suffix());
 }
 //--------------------------------------------------------------------------
@@ -116,7 +116,7 @@ super_subset_part_name(std::string base_name)
   // Note: there's a 32 character limit on the maximum length of the
   // part name (Ioss throws a warning message), so economy on characters is good
 
-  ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem part");
+  STK_ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem part");
   auto first_token = base_name.substr(0, base_name.find_first_of('_'));
   auto last_token =
     base_name.substr(base_name.find_last_of("_"), base_name.length());
@@ -135,7 +135,7 @@ super_subset_part_name(
   // Note: there's a 32 character limit on the maximum length of the
   // part name (Ioss throws a warning message), so economy on characters is good
 
-  //    ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem
+  //    STK_ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem
   //    part"); auto first_token = base_name.substr(0,
   //    base_name.find_first_of('_')); auto last_token =
   //    base_name.substr(base_name.find_last_of("_"),base_name.length());
@@ -144,7 +144,7 @@ super_subset_part_name(
   //        + "_superside" + std::to_string(numSideNodes)
   //        + last_token;
 
-  ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem part");
+  STK_ThrowAssertMsg(!base_name.empty(), "Empty base name for super elem part");
   auto first_token = base_name.substr(0, base_name.find_first_of('_'));
   auto last_token =
     base_name.substr(base_name.find_last_of("_"), base_name.length());
@@ -157,7 +157,7 @@ super_subset_part_name(
 std::string
 base_element_part_name(std::string super_name)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     super_name.find(super_element_suffix()) != std::string::npos,
     "Not a super-element part name!");
   return (super_name.erase(
@@ -194,7 +194,7 @@ base_elem_part_from_super_elem_part(const stk::mesh::Part& super_elem_part)
 stk::mesh::Part*
 super_elem_part(const stk::mesh::Part* part)
 {
-  ThrowAssert(part != nullptr);
+  STK_ThrowAssert(part != nullptr);
   return (
     part->mesh_meta_data().get_part(super_element_part_name(part->name())));
 }
@@ -202,7 +202,7 @@ super_elem_part(const stk::mesh::Part* part)
 void
 transform_to_super_elem_part_vector(stk::mesh::PartVector& parts)
 {
-  ThrowAssert(part_vector_is_valid_and_nonempty(parts));
+  STK_ThrowAssert(part_vector_is_valid_and_nonempty(parts));
   std::transform(
     parts.begin(), parts.end(), parts.begin(),
     [](stk::mesh::Part* part) { return super_elem_part(part); });
@@ -234,7 +234,7 @@ base_ranked_parts(
   stk::topology::rank_t rank,
   bool with_subsets)
 {
-  ThrowAssert(part_vector_is_valid_and_nonempty(parts));
+  STK_ThrowAssert(part_vector_is_valid_and_nonempty(parts));
 
   stk::mesh::PartVector elemParts;
   std::copy_if(
@@ -276,7 +276,7 @@ base_face_parts(const stk::mesh::PartVector& parts)
 stk::mesh::PartVector
 only_super_parts(const stk::mesh::PartVector& parts)
 {
-  ThrowAssert(part_vector_is_valid_and_nonempty(parts));
+  STK_ThrowAssert(part_vector_is_valid_and_nonempty(parts));
 
   stk::mesh::PartVector elemParts;
   std::copy_if(
@@ -298,7 +298,7 @@ only_super_parts(const stk::mesh::PartVector& parts)
 stk::mesh::PartVector
 only_super_elem_parts(const stk::mesh::PartVector& parts)
 {
-  ThrowAssert(part_vector_is_valid_and_nonempty(parts));
+  STK_ThrowAssert(part_vector_is_valid_and_nonempty(parts));
 
   stk::mesh::PartVector elemParts;
   std::copy_if(
@@ -310,7 +310,7 @@ only_super_elem_parts(const stk::mesh::PartVector& parts)
 stk::mesh::PartVector
 only_super_side_parts(const stk::mesh::PartVector& parts)
 {
-  ThrowAssert(part_vector_is_valid_and_nonempty(parts));
+  STK_ThrowAssert(part_vector_is_valid_and_nonempty(parts));
 
   stk::mesh::PartVector faceParts;
   std::copy_if(
diff --git a/src/kernel/MomentumOpenAdvDiffElemKernel.C b/src/kernel/MomentumOpenAdvDiffElemKernel.C
index 21f1a437a..a95fcb430 100644
--- a/src/kernel/MomentumOpenAdvDiffElemKernel.C
+++ b/src/kernel/MomentumOpenAdvDiffElemKernel.C
@@ -470,7 +470,7 @@ MomentumOpenAdvDiffElemKernel<BcAlgTraits>::execute(
         break;
       }
       default:
-        NGP_ThrowErrorMsg("invalid entrainment method");
+        STK_NGP_ThrowErrorMsg("invalid entrainment method");
       }
     }
 
diff --git a/src/master_element/Hex8CVFEM.C b/src/master_element/Hex8CVFEM.C
index 697bc82be..a738805e9 100644
--- a/src/master_element/Hex8CVFEM.C
+++ b/src/master_element/Hex8CVFEM.C
@@ -51,13 +51,13 @@ hex8_derivative(
       "hex8_derivative: Error in deriv or par_coord array index 0");
 #else
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     8 == deriv.extent(1), "hex8_derivative: Error in derivative array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == deriv.extent(2), "hex8_derivative: Error in derivative array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == par_coord.extent(1), "hex8_derivative: Error in derivative array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     deriv.extent(0) == par_coord.extent(0),
     "hex8_derivative: Error in derivative array");
 
@@ -141,34 +141,34 @@ hex_gradient_operator(
   //
   const unsigned nint = deriv.extent(0);
   const unsigned npe = deriv.extent(1);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == deriv.extent(2), "hex_gradient_operator: Error in derivative array");
 
   const unsigned nelem = cordel.extent(0);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == cordel.extent(1),
     "hex_gradient_operator: Error in coorindate array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == cordel.extent(2), "hex_gradient_operator: Error in coorindate array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == gradop.extent(0), "hex_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == gradop.extent(1),
     "hex_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == gradop.extent(2), "hex_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == gradop.extent(3), "hex_gradient_operator: Error in gradient array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == det_j.extent(0),
     "hex_gradient_operator: Error in determinent array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == det_j.extent(1),
     "hex_gradient_operator: Error in determinent array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == err.extent(0), "hex_gradient_operator: Error in error array");
 
   const double realmin = std::numeric_limits<double>::min();
diff --git a/src/master_element/Hex8FEM.C b/src/master_element/Hex8FEM.C
index 090a33d76..a39477352 100644
--- a/src/master_element/Hex8FEM.C
+++ b/src/master_element/Hex8FEM.C
@@ -60,11 +60,12 @@ generic_grad_op_3d(
   static_assert(OutputViewType::rank == 3, "Weight view assumed to be 3D");
   static_assert(AlgTraits::nDim_ == 3, "3D method");
 
-  ThrowAssert(AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
-  ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
-  ThrowAssert(weights.extent(0) == referenceGradWeights.extent(0));
-  ThrowAssert(weights.extent(1) == referenceGradWeights.extent(1));
-  ThrowAssert(weights.extent(2) == referenceGradWeights.extent(2));
+  STK_ThrowAssert(
+    AlgTraits::nodesPerElement_ == referenceGradWeights.extent(1));
+  STK_ThrowAssert(AlgTraits::nDim_ == referenceGradWeights.extent(2));
+  STK_ThrowAssert(weights.extent(0) == referenceGradWeights.extent(0));
+  STK_ThrowAssert(weights.extent(1) == referenceGradWeights.extent(1));
+  STK_ThrowAssert(weights.extent(2) == referenceGradWeights.extent(2));
 
   for (unsigned ip = 0; ip < AlgTraits::numGp_; ++ip) {
     ftype jact[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
@@ -95,7 +96,7 @@ generic_grad_op_3d(
 
     detj(ip) = jact[0][0] * invJac[0][0] + jact[1][0] * invJac[1][0] +
                jact[2][0] * invJac[2][0];
-    ThrowAssertMsg(
+    STK_ThrowAssertMsg(
       stk::simd::are_any(detj(ip) > +tiny_positive_value()),
       "Problem with determinant");
 
diff --git a/src/master_element/MasterElementRepo.C b/src/master_element/MasterElementRepo.C
index 308193eb9..2a4bda8ef 100644
--- a/src/master_element/MasterElementRepo.C
+++ b/src/master_element/MasterElementRepo.C
@@ -194,7 +194,7 @@ find_topo(
          "pointer."
       << " The pointer was not found in the master element database:"
       << host_ptr << std::endl;
-    ThrowRequire(host_map.find(host_ptr) != dev_map.end());
+    STK_ThrowRequire(host_map.find(host_ptr) != dev_map.end());
   }
   return theTopo;
 }
@@ -215,7 +215,7 @@ find_me(
       NaluEnv::self().naluOutputP0() << v.first.name() << std::endl;
     NaluEnv::self().naluOutputP0()
       << " Add topology to MasterElementRepo::find_me()" << std::endl;
-    ThrowRequire(it != me_map.end());
+    STK_ThrowRequire(it != me_map.end());
   }
   MasterElement* theElem = it->second;
   return theElem;
diff --git a/src/master_element/Pyr5CVFEM.C b/src/master_element/Pyr5CVFEM.C
index bf59d1c9b..00803cc1b 100644
--- a/src/master_element/Pyr5CVFEM.C
+++ b/src/master_element/Pyr5CVFEM.C
@@ -65,34 +65,34 @@ pyr_gradient_operator(
   //
   const unsigned nint = deriv.extent(0);
   const unsigned npe = deriv.extent(1);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == deriv.extent(2), "pyr_gradient_operator: Error in derivative array");
 
   const unsigned nelem = cordel.extent(0);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == cordel.extent(1),
     "pyr_gradient_operator: Error in coorindate array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == cordel.extent(2), "pyr_gradient_operator: Error in coorindate array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == gradop.extent(0), "pyr_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == gradop.extent(1),
     "pyr_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == gradop.extent(2), "pyr_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == gradop.extent(3), "pyr_gradient_operator: Error in gradient array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == det_j.extent(0),
     "pyr_gradient_operator: Error in determinent array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == det_j.extent(1),
     "pyr_gradient_operator: Error in determinent array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == err.extent(0), "pyr_gradient_operator: Error in error array");
 
   const double realmin = std::numeric_limits<double>::min();
@@ -1429,7 +1429,7 @@ PyrSCS::sidePcoords_to_elemPcoords(
   const double* side_pcoords,
   double* elem_pcoords)
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     side_ordinal >= 0 && side_ordinal <= 4,
     "Invalid pyramid side ordinal " + std::to_string(side_ordinal));
 
diff --git a/src/master_element/Tet4CVFEM.C b/src/master_element/Tet4CVFEM.C
index fd1be623a..07f05f4d8 100644
--- a/src/master_element/Tet4CVFEM.C
+++ b/src/master_element/Tet4CVFEM.C
@@ -89,34 +89,34 @@ tet_gradient_operator(
   //
   const unsigned nint = deriv.extent(0);
   const unsigned npe = deriv.extent(1);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == deriv.extent(2), "tet_gradient_operator: Error in derivative array");
 
   const unsigned nelem = cordel.extent(0);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == cordel.extent(1),
     "tet_gradient_operator: Error in coorindate array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == cordel.extent(2), "tet_gradient_operator: Error in coorindate array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == gradop.extent(0), "tet_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == gradop.extent(1),
     "tet_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == gradop.extent(2), "tet_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == gradop.extent(3), "tet_gradient_operator: Error in gradient array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == det_j.extent(0),
     "tet_gradient_operator: Error in determinent array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == det_j.extent(1),
     "tet_gradient_operator: Error in determinent array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == err.extent(0), "tet_gradient_operator: Error in error array");
 
   const double realmin = std::numeric_limits<double>::min();
@@ -336,7 +336,7 @@ TetSCV::determinant_scv(
     // compute volume using an equivalent polyhedron
     volume(icv) = hex_volume_grandy(ehexcoords);
     // check for negative volume
-    // ThrowAssertMsg( volume(icv) < 0.0, "ERROR in TetSCV::determinant,
+    // STK_ThrowAssertMsg( volume(icv) < 0.0, "ERROR in TetSCV::determinant,
     // negative volume.");
   }
 }
diff --git a/src/master_element/Wed6CVFEM.C b/src/master_element/Wed6CVFEM.C
index 96d1f9bc0..fea1e7ee8 100644
--- a/src/master_element/Wed6CVFEM.C
+++ b/src/master_element/Wed6CVFEM.C
@@ -52,34 +52,34 @@ wed_gradient_operator(
   //
   const unsigned nint = deriv.extent(0);
   const unsigned npe = deriv.extent(1);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == deriv.extent(2), "wed_gradient_operator: Error in derivative array");
 
   const unsigned nelem = cordel.extent(0);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == cordel.extent(1),
     "wed_gradient_operator: Error in coorindate array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == cordel.extent(2), "wed_gradient_operator: Error in coorindate array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == gradop.extent(0), "wed_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == gradop.extent(1),
     "wed_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     npe == gradop.extent(2), "wed_gradient_operator: Error in gradient array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     3 == gradop.extent(3), "wed_gradient_operator: Error in gradient array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nint == det_j.extent(0),
     "wed_gradient_operator: Error in determinent array");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == det_j.extent(1),
     "wed_gradient_operator: Error in determinent array");
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     nelem == err.extent(0), "wed_gradient_operator: Error in error array");
 
   const double realmin = std::numeric_limits<double>::min();
diff --git a/src/matrix_free/ConductionFields.C b/src/matrix_free/ConductionFields.C
index bc6426142..2bc4e4363 100644
--- a/src/matrix_free/ConductionFields.C
+++ b/src/matrix_free/ConductionFields.C
@@ -28,8 +28,8 @@ get_ngp_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
   auto field = stk::mesh::get_updated_ngp_field<double>(
     *meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
diff --git a/src/matrix_free/ConductionGatheredFieldManager.C b/src/matrix_free/ConductionGatheredFieldManager.C
index 6878662d0..7466cfaae 100644
--- a/src/matrix_free/ConductionGatheredFieldManager.C
+++ b/src/matrix_free/ConductionGatheredFieldManager.C
@@ -33,8 +33,8 @@ get_ngp_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
   return stk::mesh::get_updated_ngp_field<double>(
     *meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
diff --git a/src/matrix_free/ConductionOperator.C b/src/matrix_free/ConductionOperator.C
index 79836f46c..071a73628 100644
--- a/src/matrix_free/ConductionOperator.C
+++ b/src/matrix_free/ConductionOperator.C
@@ -103,9 +103,9 @@ ConductionLinearizedResidualOperator<p>::apply(
   double beta) const
 {
   stk::mesh::ProfilingBlock pf("LinearizedResidualOperator<p>::apply");
-  ThrowRequire(trans == Teuchos::NO_TRANS);
-  ThrowRequire(alpha == 1.0);
-  ThrowRequire(beta == 0.0);
+  STK_ThrowRequire(trans == Teuchos::NO_TRANS);
+  STK_ThrowRequire(alpha == 1.0);
+  STK_ThrowRequire(beta == 0.0);
   if (exporter_.getTargetMap()->isDistributed()) {
     cached_sln_.doImport(owned_sln, exporter_, Tpetra::INSERT);
     cached_rhs_.putScalar(0.);
diff --git a/src/matrix_free/ContinuityOperator.C b/src/matrix_free/ContinuityOperator.C
index cdaa750d4..98614237d 100644
--- a/src/matrix_free/ContinuityOperator.C
+++ b/src/matrix_free/ContinuityOperator.C
@@ -109,9 +109,9 @@ ContinuityLinearizedResidualOperator<p>::apply(
 {
   stk::mesh::ProfilingBlock pf(
     "ContinuityLinearizedResidualOperator<p>::apply");
-  ThrowRequire(trans == Teuchos::NO_TRANS);
-  ThrowRequire(alpha == 1.0);
-  ThrowRequire(beta == 0.0);
+  STK_ThrowRequire(trans == Teuchos::NO_TRANS);
+  STK_ThrowRequire(alpha == 1.0);
+  STK_ThrowRequire(beta == 0.0);
 
   if (exporter_.getTargetMap()->isDistributed()) {
     {
diff --git a/src/matrix_free/GreenGaussGradientOperator.C b/src/matrix_free/GreenGaussGradientOperator.C
index e78e026a4..7f4f7a1a3 100644
--- a/src/matrix_free/GreenGaussGradientOperator.C
+++ b/src/matrix_free/GreenGaussGradientOperator.C
@@ -90,7 +90,7 @@ GradientLinearizedResidualOperator<p>::apply(
 {
   stk::mesh::ProfilingBlock pf("GradientLinearizedResidualOperator<p>::apply");
 
-  ThrowRequire(owned_sln.getNumVectors() == 3);
+  STK_ThrowRequire(owned_sln.getNumVectors() == 3);
 
   if (exporter_.getTargetMap()->isDistributed()) {
     {
diff --git a/src/matrix_free/LowMachFields.C b/src/matrix_free/LowMachFields.C
index b424aa244..05aebc338 100644
--- a/src/matrix_free/LowMachFields.C
+++ b/src/matrix_free/LowMachFields.C
@@ -33,8 +33,8 @@ get_and_sync_ngp_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
 
   auto field = stk::mesh::get_updated_ngp_field<T>(
diff --git a/src/matrix_free/LowMachGatheredFieldManager.C b/src/matrix_free/LowMachGatheredFieldManager.C
index 8ce6a6749..60afffd21 100644
--- a/src/matrix_free/LowMachGatheredFieldManager.C
+++ b/src/matrix_free/LowMachGatheredFieldManager.C
@@ -48,8 +48,8 @@ get_synced_ngp_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
 
   auto field = stk::mesh::get_updated_ngp_field<T>(
diff --git a/src/matrix_free/LowMachUpdate.C b/src/matrix_free/LowMachUpdate.C
index a9f408eff..bbec658bb 100644
--- a/src/matrix_free/LowMachUpdate.C
+++ b/src/matrix_free/LowMachUpdate.C
@@ -33,8 +33,8 @@ get_ngp_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
   return stk::mesh::get_updated_ngp_field<T>(
     *meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
diff --git a/src/matrix_free/MomentumOperator.C b/src/matrix_free/MomentumOperator.C
index ca7eee861..6d628c796 100644
--- a/src/matrix_free/MomentumOperator.C
+++ b/src/matrix_free/MomentumOperator.C
@@ -59,7 +59,7 @@ MomentumResidualOperator<p>::compute(mv_type& owned_rhs)
 {
   stk::mesh::ProfilingBlock pf("MomentumResidualOperator<p>::apply");
   if (exporter_.getTargetMap()->isDistributed()) {
-    ThrowRequire(owned_rhs.getLocalLength() == size_t(max_owned_row_id_));
+    STK_ThrowRequire(owned_rhs.getLocalLength() == size_t(max_owned_row_id_));
     local_compute(cached_rhs_.getLocalViewDevice(Tpetra::Access::ReadWrite));
     {
       stk::mesh::ProfilingBlock pfinner("export from owned-shared to owned");
@@ -117,16 +117,16 @@ MomentumLinearizedResidualOperator<p>::apply(
   double beta) const
 {
   stk::mesh::ProfilingBlock pf("MomentumLinearizedResidualOperator<p>::apply");
-  ThrowRequire(trans == Teuchos::NO_TRANS);
-  ThrowRequire(alpha == 1.0);
-  ThrowRequire(beta == 0.0);
+  STK_ThrowRequire(trans == Teuchos::NO_TRANS);
+  STK_ThrowRequire(alpha == 1.0);
+  STK_ThrowRequire(beta == 0.0);
   if (exporter_.getTargetMap()->isDistributed()) {
     {
       stk::mesh::ProfilingBlock pfinner("import into owned-shared from owned");
       cached_sln_.doImport(owned_sln, exporter_, Tpetra::INSERT);
     }
 
-    ThrowRequire(owned_rhs.getLocalLength() == size_t(max_owned_row_id_));
+    STK_ThrowRequire(owned_rhs.getLocalLength() == size_t(max_owned_row_id_));
     local_apply(
       cached_sln_.getLocalViewDevice(Tpetra::Access::ReadWrite),
       cached_rhs_.getLocalViewDevice(Tpetra::Access::ReadWrite));
diff --git a/src/matrix_free/StkToTpetraMap.C b/src/matrix_free/StkToTpetraMap.C
index a0f379f7e..57c6b0fa4 100644
--- a/src/matrix_free/StkToTpetraMap.C
+++ b/src/matrix_free/StkToTpetraMap.C
@@ -195,7 +195,7 @@ make_owned_and_shared_row_map(
     (meta.globally_shared_part() & active_linsys) - meta.locally_owned_part();
   const auto num_shared = count_local_nodes(mesh, shared_selector);
 
-  ThrowRequire(
+  STK_ThrowRequire(
     count_local_nodes(mesh, active_linsys) == num_owned + num_shared);
   Kokkos::View<global_ordinal_type*> row_ids(
     "oas_row_ids", count_local_nodes(mesh, active_linsys));
diff --git a/src/matrix_free/StrongDirichletBC.C b/src/matrix_free/StrongDirichletBC.C
index b81838ae6..19d93d058 100644
--- a/src/matrix_free/StrongDirichletBC.C
+++ b/src/matrix_free/StrongDirichletBC.C
@@ -55,7 +55,8 @@ dirichlet_residual(
   tpetra_view_type yout)
 {
   stk::mesh::ProfilingBlock pf("vector_dirichlet_residual");
-  ThrowRequireMsg(yout.extent_int(1) == 3, "length is " << yout.extent_int(1));
+  STK_ThrowRequireMsg(
+    yout.extent_int(1) == 3, "length is " << yout.extent_int(1));
 
 #if defined(KOKKOS_ENABLE_HIP)
   using policy_type = Kokkos::MDRangePolicy<
@@ -85,8 +86,8 @@ dirichlet_linearized(
   tpetra_view_type yout)
 {
   stk::mesh::ProfilingBlock pf("dirichlet_linearized");
-  ThrowRequire(yout.extent_int(0) == xin.extent_int(0));
-  ThrowRequire(yout.extent_int(1) == xin.extent_int(1));
+  STK_ThrowRequire(yout.extent_int(0) == xin.extent_int(0));
+  STK_ThrowRequire(yout.extent_int(1) == xin.extent_int(1));
 
 #if defined(KOKKOS_ENABLE_HIP)
   using policy_type = Kokkos::MDRangePolicy<
diff --git a/src/mesh_motion/FrameBase.C b/src/mesh_motion/FrameBase.C
index 00b05d1e6..36ae39be3 100644
--- a/src/mesh_motion/FrameBase.C
+++ b/src/mesh_motion/FrameBase.C
@@ -110,7 +110,7 @@ FrameBase::populate_part_vec(const YAML::Node& node)
     partNamesVec.end()) {
     partNamesVec.clear();
     for (const auto* part : meta_.get_mesh_parts()) {
-      ThrowRequire(part);
+      STK_ThrowRequire(part);
       if (part->topology().rank() == stk::topology::ELEMENT_RANK) {
         partNamesVec.push_back(part->name());
       }
diff --git a/src/ngp_algorithms/NodalGradBndryElemAlg.C b/src/ngp_algorithms/NodalGradBndryElemAlg.C
index 60207bdad..c877cfec8 100644
--- a/src/ngp_algorithms/NodalGradBndryElemAlg.C
+++ b/src/ngp_algorithms/NodalGradBndryElemAlg.C
@@ -53,21 +53,21 @@ NodalGradBndryElemAlg<AlgTraits, PhiType, GradPhiType, ViewHelperType>::
       MasterElementRepo::get_surface_master_element_on_dev(AlgTraits::topo_))
 {
   if (phiSize_ == 1u) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       gradPhiSize_ == AlgTraits::nDim_,
       "NodalGradBndryElemAlg called with scalar input field '"
         << phi->name() << "' but with non-vector output field '"
         << gradPhi->name() << "' of length " << gradPhiSize_ << " (should be "
         << AlgTraits::nDim_ << ")");
   } else if (phiSize_ == AlgTraits::nDim_) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       gradPhiSize_ == AlgTraits::nDim_ * AlgTraits::nDim_,
       "NodalGradBndryElemAlg called with vector input field '"
         << phi->name() << "' but with non-tensor output field '"
         << gradPhi->name() << "' of length " << gradPhiSize_ << " (should be "
         << AlgTraits::nDim_ * AlgTraits::nDim_ << ")");
   } else {
-    ThrowErrorMsg(
+    STK_ThrowErrorMsg(
       "NodalGradBndryElemAlg called with an input field '"
       << phi->name()
       << "' that is not a scalar or a vector.  "
diff --git a/src/ngp_algorithms/NodalGradEdgeAlg.C b/src/ngp_algorithms/NodalGradEdgeAlg.C
index 7a53cfcf5..e1257f26b 100644
--- a/src/ngp_algorithms/NodalGradEdgeAlg.C
+++ b/src/ngp_algorithms/NodalGradEdgeAlg.C
@@ -32,21 +32,21 @@ NodalGradEdgeAlg<PhiType, GradPhiType>::NodalGradEdgeAlg(
 {
   const int gradPhiSize = max_extent(*gradPhi, 0);
   if (dim1_ == 1) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       gradPhiSize == dim2_, "NodalGradEdgeAlg called with scalar input field '"
                               << phi->name()
                               << "' but with non-vector output field '"
                               << gradPhi->name() << "' of length "
                               << gradPhiSize << " (should be " << dim2_ << ")");
   } else if (dim1_ == dim2_) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       gradPhiSize == dim2_ * dim2_,
       "NodalGradBndryElemAlg called with vector input field '"
         << phi->name() << "' but with non-tensor output field '"
         << gradPhi->name() << "' of length " << gradPhiSize << " (should be "
         << dim2_ * dim2_ << ")");
   } else {
-    ThrowErrorMsg(
+    STK_ThrowErrorMsg(
       "NodalGradBndryElemAlg called with an input field '"
       << phi->name()
       << "' that is not a scalar or a vector.  "
diff --git a/src/ngp_algorithms/NodalGradElemAlg.C b/src/ngp_algorithms/NodalGradElemAlg.C
index bdcd702d4..27882deb1 100644
--- a/src/ngp_algorithms/NodalGradElemAlg.C
+++ b/src/ngp_algorithms/NodalGradElemAlg.C
@@ -49,21 +49,21 @@ NodalGradElemAlg<AlgTraits, PhiType, GradPhiType, ViewHelperType>::
       MasterElementRepo::get_surface_master_element_on_dev(AlgTraits::topo_))
 {
   if (phiSize_ == 1u) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       gradPhiSize_ == AlgTraits::nDim_,
       "NodalGradElemAlg called with scalar input field '"
         << phi->name() << "' but with non-vector output field '"
         << gradPhi->name() << "' of length " << gradPhiSize_ << " (should be "
         << AlgTraits::nDim_ << ")");
   } else if (phiSize_ == AlgTraits::nDim_) {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       gradPhiSize_ == AlgTraits::nDim_ * AlgTraits::nDim_,
       "NodalGradElemAlg called with vector input field '"
         << phi->name() << "' but with non-tensor output field '"
         << gradPhi->name() << "' of length " << gradPhiSize_ << " (should be "
         << AlgTraits::nDim_ * AlgTraits::nDim_ << ")");
   } else {
-    ThrowErrorMsg(
+    STK_ThrowErrorMsg(
       "NodalGradBndryElemAlg called with an input field '"
       << phi->name()
       << "' that is not a scalar or a vector.  "
diff --git a/src/overset/ExtOverset.C b/src/overset/ExtOverset.C
index d6274d2a1..b4665a41f 100644
--- a/src/overset/ExtOverset.C
+++ b/src/overset/ExtOverset.C
@@ -181,7 +181,7 @@ ExtOverset::register_solution(const std::vector<std::string>& fnames)
   if (!hasOverset_)
     return ncomp;
 
-  ThrowAssert(fnames.size() > 0u);
+  STK_ThrowAssert(fnames.size() > 0u);
   // Store field names for update solution phase
   slnFieldNames_ = fnames;
 
@@ -206,7 +206,7 @@ ExtOverset::update_solution()
   if (!hasOverset_)
     return;
 
-  ThrowAssert(slnFieldNames_.size() > 0u);
+  STK_ThrowAssert(slnFieldNames_.size() > 0u);
 
 #ifdef NALU_USES_TIOGA
   for (auto* realm : time_.realmVec_) {
diff --git a/src/overset/OversetManager.C b/src/overset/OversetManager.C
index a720c18ac..e5e1f518e 100644
--- a/src/overset/OversetManager.C
+++ b/src/overset/OversetManager.C
@@ -97,7 +97,7 @@ OversetManager::overset_orphan_node_field_update(
     const int num_nodes = bulkData_->num_nodes(owningElement);
 
     // sanity check on num nodes
-    ThrowAssert(num_nodes == nodesPerElement);
+    STK_ThrowAssert(num_nodes == nodesPerElement);
 
     for (int ni = 0; ni < num_nodes; ++ni) {
       stk::mesh::Entity node = elem_node_rels[ni];
diff --git a/src/overset/OversetManagerTIOGA.C b/src/overset/OversetManagerTIOGA.C
index e67e058bc..a7c04b5fa 100644
--- a/src/overset/OversetManagerTIOGA.C
+++ b/src/overset/OversetManagerTIOGA.C
@@ -35,7 +35,7 @@ OversetManagerTIOGA::OversetManagerTIOGA(
     tiogaIface_(
       *this, oversetUserData.oversetBlocks_, realm.get_coordinates_name())
 {
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     metaData_->spatial_dimension() == 3u, "TIOGA only supports 3-D meshes.");
 }
 
diff --git a/src/overset/TiogaSTKIface.C b/src/overset/TiogaSTKIface.C
index 7667ab211..935c5d85c 100644
--- a/src/overset/TiogaSTKIface.C
+++ b/src/overset/TiogaSTKIface.C
@@ -395,7 +395,7 @@ TiogaSTKIface::populate_overset_info()
   std::unordered_set<stk::mesh::EntityId> seenIDs;
 
   // Ensure that the oversetInfoVec has been cleared out
-  ThrowAssert(osetInfo.size() == 0);
+  STK_ThrowAssert(osetInfo.size() == 0);
 
   sierra::nalu::VectorFieldType* coords =
     meta_.get_field<double>(stk::topology::NODE_RANK, coordsName_);
diff --git a/src/overset/overset_utils.C b/src/overset/overset_utils.C
index c4a88090f..4b8a6dffd 100644
--- a/src/overset/overset_utils.C
+++ b/src/overset/overset_utils.C
@@ -23,9 +23,9 @@ get_overset_field_data(Realm& realm, std::vector<std::string> fnames)
 
   for (const auto& ff : fnames) {
     auto* fld = meta.get_field(stk::topology::NODE_RANK, ff);
-    ThrowAssert(fld != nullptr);
+    STK_ThrowAssert(fld != nullptr);
 
-    const int col = fld->max_size(fld->entity_rank());
+    const int col = fld->max_size();
     fields.emplace_back(fld, row, col);
   }
 
diff --git a/src/property_evaluator/GenericPropAlgorithm.C b/src/property_evaluator/GenericPropAlgorithm.C
index 3e4136dbd..8bbb9a84e 100644
--- a/src/property_evaluator/GenericPropAlgorithm.C
+++ b/src/property_evaluator/GenericPropAlgorithm.C
@@ -48,7 +48,7 @@ GenericPropAlgorithm::execute()
   }
 
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   // empty independet variable list; hence "Generic"
   std::vector<double> indVarList(1, 0.0);
diff --git a/src/property_evaluator/InverseDualVolumePropAlgorithm.C b/src/property_evaluator/InverseDualVolumePropAlgorithm.C
index 8dc1f8e68..029223c6d 100644
--- a/src/property_evaluator/InverseDualVolumePropAlgorithm.C
+++ b/src/property_evaluator/InverseDualVolumePropAlgorithm.C
@@ -38,7 +38,7 @@ InverseDualVolumePropAlgorithm::execute()
 {
 
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   stk::mesh::Selector selector = stk::mesh::selectUnion(partVec_);
 
diff --git a/src/property_evaluator/InversePropAlgorithm.C b/src/property_evaluator/InversePropAlgorithm.C
index 8a9e17752..0554f23aa 100644
--- a/src/property_evaluator/InversePropAlgorithm.C
+++ b/src/property_evaluator/InversePropAlgorithm.C
@@ -44,7 +44,7 @@ InversePropAlgorithm::execute()
 {
 
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   stk::mesh::Selector selector = stk::mesh::selectUnion(partVec_);
 
diff --git a/src/property_evaluator/LinearPropAlgorithm.C b/src/property_evaluator/LinearPropAlgorithm.C
index e97d3990d..0fb53028a 100644
--- a/src/property_evaluator/LinearPropAlgorithm.C
+++ b/src/property_evaluator/LinearPropAlgorithm.C
@@ -42,7 +42,7 @@ LinearPropAlgorithm::execute()
 {
 
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   stk::mesh::Selector selector = stk::mesh::selectUnion(partVec_);
 
diff --git a/src/property_evaluator/TemperaturePropAlgorithm.C b/src/property_evaluator/TemperaturePropAlgorithm.C
index d68e04b16..4ed3e36fe 100644
--- a/src/property_evaluator/TemperaturePropAlgorithm.C
+++ b/src/property_evaluator/TemperaturePropAlgorithm.C
@@ -48,7 +48,7 @@ TemperaturePropAlgorithm::execute()
 {
 
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   std::vector<double> indVarList(1);
 
diff --git a/src/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.C b/src/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.C
index 280fb4c2c..64bb4f657 100644
--- a/src/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.C
+++ b/src/property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.C
@@ -9,6 +9,7 @@
 
 #include <Algorithm.h>
 #include <property_evaluator/ThermalConductivityFromPrandtlPropAlgorithm.h>
+#include <FieldManager.h>
 #include <FieldTypeDef.h>
 #include <Realm.h>
 
@@ -31,19 +32,12 @@ namespace nalu {
 //--------------------------------------------------------------------------
 ThermalConductivityFromPrandtlPropAlgorithm::
   ThermalConductivityFromPrandtlPropAlgorithm(
-    Realm& realm,
-    const stk::mesh::PartVector& part_vec,
-    ScalarFieldType* thermalCond,
-    ScalarFieldType* specHeat,
-    ScalarFieldType* viscosity,
-    const double Pr)
-  : Algorithm(realm, part_vec),
-    thermalCond_(thermalCond),
-    specHeat_(specHeat),
-    viscosity_(viscosity),
-    Pr_(Pr)
+    Realm& realm, const stk::mesh::PartVector& part_vec, const double Pr)
+  : Algorithm(realm, part_vec), Pr_(Pr)
 {
-  // does nothing
+  fieldManager_.register_field<double>("thermal_conductivity", part_vec);
+  fieldManager_.register_field<double>("specific_heat", part_vec);
+  fieldManager_.register_field<double>("viscosity", part_vec);
 }
 
 //--------------------------------------------------------------------------
@@ -53,32 +47,30 @@ void
 ThermalConductivityFromPrandtlPropAlgorithm::execute()
 {
   // make sure that partVec_ is size one
-  ThrowAssert(partVec_.size() == 1);
+  STK_ThrowAssert(partVec_.size() == 1);
 
   stk::mesh::Selector selector = stk::mesh::selectUnion(partVec_);
 
   stk::mesh::BucketVector const& node_buckets =
     realm_.get_buckets(stk::topology::NODE_RANK, selector);
 
-  thermalCond_->sync_to_host();
-  specHeat_->sync_to_host();
-  viscosity_->sync_to_host();
+  auto thermalCond =
+    fieldManager_.get_legacy_smart_field<double, tags::READ_WRITE>(
+      "thermal_conductivity");
+  const auto specHeat =
+    fieldManager_.get_legacy_smart_field<double, tags::READ>("specific_heat");
+  const auto viscosity =
+    fieldManager_.get_legacy_smart_field<double, tags::READ>("viscosity");
 
   for (stk::mesh::BucketVector::const_iterator ib = node_buckets.begin();
        ib != node_buckets.end(); ++ib) {
     stk::mesh::Bucket& b = **ib;
     const stk::mesh::Bucket::size_type length = b.size();
 
-    double* thermalCond = stk::mesh::field_data(*thermalCond_, b);
-    const double* specHeat = stk::mesh::field_data(*specHeat_, b);
-    const double* viscosity = stk::mesh::field_data(*viscosity_, b);
-
     for (stk::mesh::Bucket::size_type k = 0; k < length; ++k) {
-      thermalCond[k] = specHeat[k] * viscosity[k] / Pr_;
+      thermalCond(b)[k] = specHeat(b)[k] * viscosity(b)[k] / Pr_;
     }
   }
-  thermalCond_->modify_on_host();
-  thermalCond_->sync_to_device();
 }
 
 } // namespace nalu
diff --git a/src/user_functions/KovasznayVelocityAuxFunction.C b/src/user_functions/KovasznayVelocityAuxFunction.C
index 354bae2ed..85fad22e5 100644
--- a/src/user_functions/KovasznayVelocityAuxFunction.C
+++ b/src/user_functions/KovasznayVelocityAuxFunction.C
@@ -26,7 +26,7 @@ KovasznayVelocityAuxFunction::KovasznayVelocityAuxFunction(
     kx_(2 * std::acos(-1.)),
     ky_(2 * std::acos(-1.))
 {
-  ThrowRequireMsg(endPos == 2, "Only 2D for Kovasznay flow");
+  STK_ThrowRequireMsg(endPos == 2, "Only 2D for Kovasznay flow");
 }
 
 void
diff --git a/src/utils/ComputeVectorDivergence.C b/src/utils/ComputeVectorDivergence.C
index 0478792d5..c087a23cf 100644
--- a/src/utils/ComputeVectorDivergence.C
+++ b/src/utils/ComputeVectorDivergence.C
@@ -196,7 +196,7 @@ compute_vector_divergence(
       int num_nodes = b->num_nodes(k);
 
       // sanity check on num nodes
-      ThrowAssert(num_nodes == nodesPerFace);
+      STK_ThrowAssert(num_nodes == nodesPerFace);
 
       for (int ni = 0; ni < num_nodes; ++ni) {
         stk::mesh::Entity node = face_node_rels[ni];
diff --git a/src/utils/StkHelpers.C b/src/utils/StkHelpers.C
index 379bf9350..ae95f84cb 100644
--- a/src/utils/StkHelpers.C
+++ b/src/utils/StkHelpers.C
@@ -61,14 +61,14 @@ get_elem_topo(const Realm& realm, const stk::mesh::Part& surfacePart)
   std::vector<const stk::mesh::Part*> blockParts =
     realm.meta_data().get_blocks_touching_surface(&surfacePart);
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     blockParts.size() >= 1,
     "Error, expected at least 1 block for surface " << surfacePart.name());
 
   stk::topology elemTopo = blockParts[0]->topology();
   if (blockParts.size() > 1) {
     for (size_t i = 1; i < blockParts.size(); ++i) {
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         blockParts[i]->topology() == elemTopo,
         "Error, found blocks of different topology connected to surface '"
           << surfacePart.name() << "', " << elemTopo << " and "
@@ -76,7 +76,7 @@ get_elem_topo(const Realm& realm, const stk::mesh::Part& surfacePart)
     }
   }
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     elemTopo != stk::topology::INVALID_TOPOLOGY,
     "Error, didn't find valid topology block for surface "
       << surfacePart.name());
diff --git a/src/wind_energy/ABLForcingAlgorithm.C b/src/wind_energy/ABLForcingAlgorithm.C
index 668000f57..8c885300d 100644
--- a/src/wind_energy/ABLForcingAlgorithm.C
+++ b/src/wind_energy/ABLForcingAlgorithm.C
@@ -165,7 +165,7 @@ ABLForcingAlgorithm::create_interp_arrays(
 
   // Check that all timesteps contain values for all the heights
   for (auto vx : inpArr) {
-    ThrowAssert((vx.size() == (nHeights + 1)));
+    STK_ThrowAssert((vx.size() == (nHeights + 1)));
   }
 
   auto nTimes = inpArr.size();
diff --git a/src/wind_energy/LidarPatterns.C b/src/wind_energy/LidarPatterns.C
index c29d1abb2..188e226f0 100644
--- a/src/wind_energy/LidarPatterns.C
+++ b/src/wind_energy/LidarPatterns.C
@@ -176,7 +176,8 @@ ScanningLidarSegmentGenerator::load(const YAML::Node& node)
   double sweep_angle_in_degrees = 20;
   get_if_present(
     node, "sweep_angle", sweep_angle_in_degrees, sweep_angle_in_degrees);
-  ThrowRequireMsg(sweep_angle_in_degrees > 0, "Sweep angle must be positive");
+  STK_ThrowRequireMsg(
+    sweep_angle_in_degrees > 0, "Sweep angle must be positive");
   sweep_angle_ = convert::degrees_to_radians(sweep_angle_in_degrees);
 
   double step_delta_angle_in_degrees = 1;
@@ -185,17 +186,18 @@ ScanningLidarSegmentGenerator::load(const YAML::Node& node)
     step_delta_angle_in_degrees);
   step_delta_angle_ = convert::degrees_to_radians(step_delta_angle_in_degrees);
 
-  ThrowRequireMsg(step_delta_angle_ > 0, "step delta angle must be positive");
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
+    step_delta_angle_ > 0, "step delta angle must be positive");
+  STK_ThrowRequireMsg(
     step_delta_angle_ <= sweep_angle_,
     "step delta angle must be less than full sweep");
 
   get_if_present(node, "stare_time", stare_time_, stare_time_);
-  ThrowRequireMsg(stare_time_ > 0, "stare time must be positive");
+  STK_ThrowRequireMsg(stare_time_ > 0, "stare time must be positive");
 
   get_if_present(
     node, "reset_time_delta", reset_time_delta_, reset_time_delta_);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     reset_time_delta_ >= 0, "reset time delta must be semi-positive");
 
   get_required(node, "beam_length", beam_length_);
@@ -320,8 +322,8 @@ SpinnerLidarSegmentGenerator::load(const YAML::Node& node)
   NaluEnv::self().naluOutputP0()
     << "LidarLineOfSite::SpinnerLidarSegmentGenerator::load" << std::endl;
 
-  ThrowRequireMsg(node["center"], "Lidar center must be provided");
-  ThrowRequireMsg(node["axis"], "Lidar axis must be provided");
+  STK_ThrowRequireMsg(node["center"], "Lidar center must be provided");
+  STK_ThrowRequireMsg(node["axis"], "Lidar axis must be provided");
 
   lidarCenter_ = to_array3(node["center"].as<Coordinates>());
 
@@ -365,7 +367,7 @@ SpinnerLidarSegmentGenerator::load(const YAML::Node& node)
     normalize_vec3(groundNormal_.data());
   }
 
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     std::abs(ddot(groundNormal_.data(), laserAxis_.data(), 3)) <
       small_positive_value(),
     "Ground and laser axes must be orthogonal");
@@ -561,7 +563,7 @@ RadarSegmentGenerator::load(const YAML::Node& node)
   double sweep_angle_in_degrees = 20;
   get_if_present(
     node, "sweep_angle", sweep_angle_in_degrees, sweep_angle_in_degrees);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     sweep_angle_in_degrees >= 0, "Sweep angle must be semipositive");
   sweep_angle_ = convert::degrees_to_radians(sweep_angle_in_degrees);
 
@@ -591,7 +593,7 @@ RadarSegmentGenerator::load(const YAML::Node& node)
 
   get_if_present(
     node, "reset_time_delta", reset_time_delta_, reset_time_delta_);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     reset_time_delta_ >= 0, "reset time delta must be semi-positive");
 
   if (specified_count > 0 && specified_count != 8) {
diff --git a/src/wind_energy/SyntheticLidar.C b/src/wind_energy/SyntheticLidar.C
index 7c2b25cad..73fe14877 100644
--- a/src/wind_energy/SyntheticLidar.C
+++ b/src/wind_energy/SyntheticLidar.C
@@ -517,7 +517,7 @@ LidarLineOfSite::output_cone_filtered(
   double dtratio)
 {
   const auto* radar = dynamic_cast<RadarSegmentGenerator*>(segGen.get());
-  ThrowRequire(radar);
+  STK_ThrowRequire(radar);
 
   auto seg = radar->generate(time());
   if (!(seg.valid || always_output_)) {
@@ -531,7 +531,7 @@ LidarLineOfSite::output_cone_filtered(
   const auto& rays = radar_data_.rays;
 
   const auto dn = (npoints_ - 1);
-  ThrowRequireMsg(dn >= 0, "At least two points required");
+  STK_ThrowRequireMsg(dn >= 0, "At least two points required");
 
   const vs::Vector dx(
     (seg.tip_[0] - seg.tail_[0]) / dn, (seg.tip_[1] - seg.tail_[1]) / dn,
@@ -613,7 +613,7 @@ LidarLineOfSite::output_cone_filtered(
         Ioss::FileInfo::create_path(name_);
       }
       // only text for now, check at parse
-      ThrowRequire(output_type_ == Output::TEXT);
+      STK_ThrowRequire(output_type_ == Output::TEXT);
       output_txt_los(time(), points, avg_line_velocity, npoints_, *file_);
     }
   }
@@ -903,7 +903,7 @@ LidarLOS::load(const YAML::Node& node, DataProbePostProcessing* probes)
     if (output_type == "dataprobes") {
       LidarLineOfSite lidarLOS;
       auto lidarDBSpec = lidarLOS.determine_line_of_site_info(node);
-      ThrowRequireMsg(
+      STK_ThrowRequireMsg(
         probes, "Lidar processing with dataprobe output "
                 "requires valid data probe object");
       probes->add_external_data_probe_spec_info(lidarDBSpec.release());
@@ -917,7 +917,7 @@ LidarLOS::load(const YAML::Node& node, DataProbePostProcessing* probes)
         double phi = 0;
         get_required(cone_grid_spec, "cone_angle", phi);
         phi = convert::degrees_to_radians(phi);
-        ThrowRequire(phi > 0);
+        STK_ThrowRequire(phi > 0);
 
         int nphi = 0;
         get_required(cone_grid_spec, "num_circles", nphi);
diff --git a/src/xfer/LocalVolumeSearch.C b/src/xfer/LocalVolumeSearch.C
index a84c07f73..8c9baf4a8 100644
--- a/src/xfer/LocalVolumeSearch.C
+++ b/src/xfer/LocalVolumeSearch.C
@@ -206,7 +206,7 @@ compute_local_coordinates(
     return compute_local_coordinates_t<AlgTraitsHex8>(
       bulk, coord_field, elem, point);
   default: {
-    ThrowRequire(topo == stk::topology::PYRAMID_5);
+    STK_ThrowRequire(topo == stk::topology::PYRAMID_5);
     return compute_local_coordinates_t<AlgTraitsPyr5>(
       bulk, coord_field, elem, point);
   }
@@ -271,7 +271,7 @@ interpolate_field(
     return interpolate_field_t<AlgTraitsHex8>(
       bulk, elem, field_prev, field, x, dtratio);
   default: {
-    ThrowRequire(topo == stk::topology::PYRAMID_5);
+    STK_ThrowRequire(topo == stk::topology::PYRAMID_5);
     return interpolate_field_t<AlgTraitsPyr5>(
       bulk, elem, field_prev, field, x, dtratio);
   }
diff --git a/unit_tests/UnitTestHelperObjects.h b/unit_tests/UnitTestHelperObjects.h
index dbc4acc55..5c6f35c60 100644
--- a/unit_tests/UnitTestHelperObjects.h
+++ b/unit_tests/UnitTestHelperObjects.h
@@ -205,13 +205,13 @@ struct EdgeHelperObjects : public HelperObjectsBase
   template <typename T, class... Args>
   void create(stk::mesh::Part* part, Args&&... args)
   {
-    ThrowRequire(edgeAlg == nullptr);
+    STK_ThrowRequire(edgeAlg == nullptr);
     edgeAlg = new T(realm, part, &eqSystem, std::forward<Args>(args)...);
   }
 
   virtual void execute() override
   {
-    ThrowRequire(edgeAlg != nullptr);
+    STK_ThrowRequire(edgeAlg != nullptr);
     edgeAlg->execute();
 
     Kokkos::deep_copy(linsys->hostNumSumIntoCalls_, linsys->numSumIntoCalls_);
diff --git a/unit_tests/UnitTestHexElementPromotion.C b/unit_tests/UnitTestHexElementPromotion.C
index 8ecebbbca..d35a02e91 100644
--- a/unit_tests/UnitTestHexElementPromotion.C
+++ b/unit_tests/UnitTestHexElementPromotion.C
@@ -61,7 +61,7 @@ protected:
     surfSubPart = nullptr;
     topo = stk::topology::HEX_8;
     hexPart = fixture->m_elem_parts[0];
-    ThrowRequire(hexPart != nullptr);
+    STK_ThrowRequire(hexPart != nullptr);
     coordField =
       &meta->declare_field<double>(stk::topology::NODE_RANK, "coords");
     intField =
@@ -201,7 +201,7 @@ TEST_F(PromoteElementHexTest, node_sharing)
   init(2, 1, 1, polyOrder);
 
   promote_mesh();
-  ThrowRequire(!bulk->in_modifiable_state());
+  STK_ThrowRequire(!bulk->in_modifiable_state());
 
   stk::mesh::EntityIdVector sharedNodeIds = {2, 5, 8, 11, 21, 22, 23, 24, 33};
 
diff --git a/unit_tests/UnitTestKokkosME.h b/unit_tests/UnitTestKokkosME.h
index 4e418f033..640c47a88 100644
--- a/unit_tests/UnitTestKokkosME.h
+++ b/unit_tests/UnitTestKokkosME.h
@@ -120,7 +120,7 @@ class KokkosMEViews
 #endif
   )
   {
-    ThrowAssertMsg(
+    STK_ThrowAssertMsg(
       partVec_.size() == 1, "KokkosMEViews unit-test assumes partVec_.size==1");
 
 #if !defined(KOKKOS_ENABLE_GPU)
diff --git a/unit_tests/UnitTestKokkosMEBC.h b/unit_tests/UnitTestKokkosMEBC.h
index efd51262f..2baea957d 100644
--- a/unit_tests/UnitTestKokkosMEBC.h
+++ b/unit_tests/UnitTestKokkosMEBC.h
@@ -97,9 +97,9 @@ class KokkosMEBC
 #endif
   )
   {
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       partVec_.size() == 1, "KokkosMEViews unit-test assumes partVec_.size==1");
-    ThrowRequireMsg(
+    STK_ThrowRequireMsg(
       !bulk_->get_buckets(meta_->side_rank(), *partVec_[0]).empty(),
       "part does not contain side-ranked elements");
 
diff --git a/unit_tests/UnitTestLinearSystem.h b/unit_tests/UnitTestLinearSystem.h
index 8b7f856c2..9548f8029 100644
--- a/unit_tests/UnitTestLinearSystem.h
+++ b/unit_tests/UnitTestLinearSystem.h
@@ -313,7 +313,7 @@ class TestLinearSystem : public sierra::nalu::LinearSystem
         rhs_(i) = rhs[i];
       }
       const size_t numRows = rhs.size();
-      ThrowAssert(numRows * numRows == lhs.size());
+      STK_ThrowAssert(numRows * numRows == lhs.size());
       for (size_t i = 0; i < numRows; ++i) {
         for (size_t j = 0; j < numRows; ++j) {
           lhs_(i, j) = lhs[numRows * i + j];
diff --git a/unit_tests/UnitTestScratchViews.C b/unit_tests/UnitTestScratchViews.C
index 70832f1fa..fabf64e93 100644
--- a/unit_tests/UnitTestScratchViews.C
+++ b/unit_tests/UnitTestScratchViews.C
@@ -134,7 +134,7 @@ do_the_test(
         sierra::nalu::get_shmem_view_1D<DoubleType, TeamType, ShmemType>(
           team, rhsSize);
 
-      NGP_ThrowAssert(scrviews.total_bytes() != 0);
+      STK_NGP_ThrowAssert(scrviews.total_bytes() != 0);
       const size_t bucketLen = b.size();
 
       Kokkos::parallel_for(
@@ -309,7 +309,7 @@ do_the_smdata_test(
       sierra::nalu::SharedMemData<TeamType, ShmemType> smdata(
         team, ngpMesh.get_spatial_dimension(), dataNGP, numNodes, rhsSize);
 
-      NGP_ThrowAssert(smdata.simdPrereqData.total_bytes() != 0);
+      STK_NGP_ThrowAssert(smdata.simdPrereqData.total_bytes() != 0);
       const size_t bucketLen = b.size();
 
       Kokkos::parallel_for(
diff --git a/unit_tests/UnitTestSideIsInElement.C b/unit_tests/UnitTestSideIsInElement.C
index 45a2b2ca4..3bbfc0e82 100644
--- a/unit_tests/UnitTestSideIsInElement.C
+++ b/unit_tests/UnitTestSideIsInElement.C
@@ -32,7 +32,7 @@ randomly_perturb_element_coords(
   const stk::mesh::Entity* node_rels,
   const VectorFieldType& coordField)
 {
-  int dim = coordField.max_size(stk::topology::NODE_RANK);
+  int dim = coordField.max_size();
   auto rot = unit_test_utils::random_rotation_matrix(dim, rng);
 
   std::uniform_real_distribution<double> random_perturb(0.125, 0.25);
diff --git a/unit_tests/UnitTestSideWriter.C b/unit_tests/UnitTestSideWriter.C
index eee3d6651..ec058fc1e 100644
--- a/unit_tests/UnitTestSideWriter.C
+++ b/unit_tests/UnitTestSideWriter.C
@@ -31,8 +31,9 @@ public:
     stk::mesh::put_field_on_mesh(
       *test_field, meta->universal_part(), 1, &minus_one);
 
+    double minus_one_vec[3] = {-1.0, -1.0, -1.0};
     stk::mesh::put_field_on_mesh(
-      *test_vector_field, meta->universal_part(), 3, &minus_one);
+      *test_vector_field, meta->universal_part(), 3, minus_one_vec);
     stk::io::set_field_output_type(
       *test_vector_field, stk::io::FieldOutputType::VECTOR_3D);
 
diff --git a/unit_tests/UnitTestTpetra.C b/unit_tests/UnitTestTpetra.C
index 365a632c8..d74146441 100644
--- a/unit_tests/UnitTestTpetra.C
+++ b/unit_tests/UnitTestTpetra.C
@@ -45,7 +45,7 @@ get_TpetraLinearSystem(unit_test_utils::NaluTest& naluObj)
 
   sierra::nalu::TpetraLinearSystem* tpetraLinsys =
     dynamic_cast<sierra::nalu::TpetraLinearSystem*>(linsys);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     tpetraLinsys != nullptr, "Expected TpetraLinearSystem to be non-null");
 
   return tpetraLinsys;
@@ -63,7 +63,7 @@ create_algorithm(sierra::nalu::Realm& realm, stk::mesh::Part& part)
       *eqsys, part, eqsys->solverAlgDriver_->solverAlgorithmMap_);
 
   EXPECT_TRUE(solverAlgResult.second);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     solverAlgResult.first != nullptr,
     "Error, failed to obtain non-null solver-algorithm object.");
 
@@ -88,11 +88,11 @@ get_AssembleElemSolverAlgorithm(unit_test_utils::NaluTest& naluObj)
   auto solverAlgMap = eqsys->solverAlgDriver_->solverAlgorithmMap_;
   EXPECT_EQ(1u, solverAlgMap.size());
   sierra::nalu::SolverAlgorithm* solverAlg = solverAlgMap.begin()->second;
-  ThrowRequireMsg(solverAlg != nullptr, "Error, null solver-algorithm");
+  STK_ThrowRequireMsg(solverAlg != nullptr, "Error, null solver-algorithm");
 
   sierra::nalu::AssembleElemSolverAlgorithm* assembleElemSolverAlgorithm =
     dynamic_cast<sierra::nalu::AssembleElemSolverAlgorithm*>(solverAlg);
-  ThrowRequireMsg(
+  STK_ThrowRequireMsg(
     assembleElemSolverAlgorithm != nullptr,
     "Error, failed to dynamic_cast to AssembleElemSolverAlgorithm.");
   return assembleElemSolverAlgorithm;
diff --git a/unit_tests/UnitTestTpetraHelperObjects.h b/unit_tests/UnitTestTpetraHelperObjects.h
index 663b6ce9d..d9e1cbdbf 100644
--- a/unit_tests/UnitTestTpetraHelperObjects.h
+++ b/unit_tests/UnitTestTpetraHelperObjects.h
@@ -254,13 +254,13 @@ struct TpetraHelperObjectsEdge : public TpetraHelperObjectsBase
   template <typename T, class... Args>
   void create(stk::mesh::Part* part, Args&&... args)
   {
-    ThrowRequire(edgeAlg == nullptr);
+    STK_ThrowRequire(edgeAlg == nullptr);
     edgeAlg = new T(realm, part, &eqSystem, std::forward<Args>(args)...);
   }
 
   virtual void execute() override
   {
-    ThrowRequire(edgeAlg != nullptr);
+    STK_ThrowRequire(edgeAlg != nullptr);
     linsys->buildEdgeToNodeGraph({&realm.meta_data().universal_part()});
     linsys->finalizeLinearSystem();
 
diff --git a/unit_tests/UnitTestUtils.C b/unit_tests/UnitTestUtils.C
index d70b28a45..595933aca 100644
--- a/unit_tests/UnitTestUtils.C
+++ b/unit_tests/UnitTestUtils.C
@@ -59,7 +59,7 @@ perturb_coord_hex_8(stk::mesh::BulkData& bulk, double perturbSize)
   const auto& meta = bulk.mesh_meta_data();
   const sierra::nalu::VectorFieldType* coordField =
     dynamic_cast<const sierra::nalu::VectorFieldType*>(meta.coordinate_field());
-  ThrowRequire(coordField != nullptr);
+  STK_ThrowRequire(coordField != nullptr);
 
   for (const auto* ib :
        bulk.get_buckets(stk::topology::NODE_RANK, meta.locally_owned_part())) {
@@ -166,7 +166,7 @@ create_one_element(
   stk::mesh::get_selected_entities(
     surfaceSelector, bulk.get_buckets(meta.side_rank(), surfaceSelector),
     all_faces);
-  ThrowRequire(all_faces.size() == topo.num_sides());
+  STK_ThrowRequire(all_faces.size() == topo.num_sides());
 
   bulk.modification_begin();
   for (unsigned k = 0u; k < all_faces.size(); ++k) {
diff --git a/unit_tests/UnitTestUtils.h b/unit_tests/UnitTestUtils.h
index d2cd03ac2..62aece25c 100644
--- a/unit_tests/UnitTestUtils.h
+++ b/unit_tests/UnitTestUtils.h
@@ -96,9 +96,10 @@ class Hex8Mesh : public ::testing::Test
 
     double one = 1.0;
     double zero = 0.0;
+    double zeroVec[3] = {0.0, 0.0, 0.0};
     const stk::mesh::PartVector parts(1, &meta->universal_part());
     elemCentroidField =
-      fieldManager->register_field<double>("elemCentroid", parts, &zero);
+      fieldManager->register_field<double>("elemCentroid", parts, zeroVec);
     nodalPressureField =
       fieldManager->register_field<double>("nodalPressure", parts, &one);
     discreteLaplacianOfPressure =
diff --git a/unit_tests/kernels/UnitTestKernelUtils.C b/unit_tests/kernels/UnitTestKernelUtils.C
index 0b96c4dcf..1dbc8c54a 100644
--- a/unit_tests/kernels/UnitTestKernelUtils.C
+++ b/unit_tests/kernels/UnitTestKernelUtils.C
@@ -876,7 +876,7 @@ calc_edge_area_vec(
   const auto& buckets = bulk.get_buckets(stk::topology::ELEMENT_RANK, sel);
 
   for (auto b : buckets) {
-    ThrowRequire(b->topology() == topo);
+    STK_ThrowRequire(b->topology() == topo);
 
     const auto bktlen = b->size();
     for (size_t ie = 0; ie < bktlen; ++ie) {
diff --git a/unit_tests/kernels/UnitTestKernelUtils.h b/unit_tests/kernels/UnitTestKernelUtils.h
index 6ad10f321..2f134e893 100644
--- a/unit_tests/kernels/UnitTestKernelUtils.h
+++ b/unit_tests/kernels/UnitTestKernelUtils.h
@@ -1454,7 +1454,7 @@ class VOFKernelHex8Mesh : public TestKernelHex8Mesh
     stk::mesh::put_field_on_mesh(
       *volumeOfFluid_, meta_->universal_part(), nullptr);
     stk::mesh::put_field_on_mesh(
-      *dvolumeOfFluidDx_, meta_->universal_part(), nullptr);
+      *dvolumeOfFluidDx_, meta_->universal_part(), spatialDim_, nullptr);
     stk::mesh::put_field_on_mesh(
       *velocity_, meta_->universal_part(), spatialDim_, nullptr);
     stk::io::set_field_output_type(
@@ -1486,7 +1486,7 @@ class VOFKernelHex8Mesh : public TestKernelHex8Mesh
   }
 
   sierra::nalu::ScalarFieldType* volumeOfFluid_{nullptr};
-  sierra::nalu::ScalarFieldType* dvolumeOfFluidDx_{nullptr};
+  sierra::nalu::VectorFieldType* dvolumeOfFluidDx_{nullptr};
   sierra::nalu::VectorFieldType* velocity_{nullptr};
   sierra::nalu::ScalarFieldType* density_{nullptr};
   sierra::nalu::ScalarFieldType* viscosity_{nullptr};
diff --git a/unit_tests/matrix_free/UnitTestConductionFields.C b/unit_tests/matrix_free/UnitTestConductionFields.C
index cc478d653..13130e81c 100644
--- a/unit_tests/matrix_free/UnitTestConductionFields.C
+++ b/unit_tests/matrix_free/UnitTestConductionFields.C
@@ -100,7 +100,7 @@ protected:
     stk::mesh::get_selected_entities(
       surfaceSelector, bulk->get_buckets(meta->side_rank(), surfaceSelector),
       all_faces);
-    ThrowRequire(all_faces.size() == topo.num_sides());
+    STK_ThrowRequire(all_faces.size() == topo.num_sides());
 
     bulk->modification_begin();
     for (unsigned k = 0u; k < all_faces.size(); ++k) {
diff --git a/unit_tests/matrix_free/UnitTestConductionUpdate.C b/unit_tests/matrix_free/UnitTestConductionUpdate.C
index 4f4054f8e..4bcab7ced 100644
--- a/unit_tests/matrix_free/UnitTestConductionUpdate.C
+++ b/unit_tests/matrix_free/UnitTestConductionUpdate.C
@@ -77,8 +77,8 @@ get_ngp_field(
   std::string name,
   stk::mesh::FieldState state = stk::mesh::StateNP1)
 {
-  ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
-  ThrowAssert(
+  STK_ThrowAssert(meta.get_field(stk::topology::NODE_RANK, name));
+  STK_ThrowAssert(
     meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
   return stk::mesh::get_updated_ngp_field<double>(
     *meta.get_field(stk::topology::NODE_RANK, name)->field_state(state));
diff --git a/unit_tests/matrix_free/UnitTestMomentumOperator.C b/unit_tests/matrix_free/UnitTestMomentumOperator.C
index ec7b2d297..ffa4e7f83 100644
--- a/unit_tests/matrix_free/UnitTestMomentumOperator.C
+++ b/unit_tests/matrix_free/UnitTestMomentumOperator.C
@@ -150,7 +150,7 @@ TEST_F(
 {
   {
     auto host_lhs = lhs.getLocalViewHost(Tpetra::Access::ReadWrite);
-    ThrowRequire(host_lhs.extent(1) == 3u);
+    STK_ThrowRequire(host_lhs.extent(1) == 3u);
     for (const auto* ib :
          bulk.get_buckets(stk::topology::NODE_RANK, active())) {
       for (auto node : *ib) {
diff --git a/unit_tests/matrix_free/UnitTestStkSimdConnectivityMap.C b/unit_tests/matrix_free/UnitTestStkSimdConnectivityMap.C
index 3939b8754..b7040e3ff 100644
--- a/unit_tests/matrix_free/UnitTestStkSimdConnectivityMap.C
+++ b/unit_tests/matrix_free/UnitTestStkSimdConnectivityMap.C
@@ -96,7 +96,7 @@ protected:
     stk::mesh::get_selected_entities(
       surfaceSelector, bulk.get_buckets(meta.side_rank(), surfaceSelector),
       all_faces);
-    ThrowRequire(all_faces.size() == topo.num_sides());
+    STK_ThrowRequire(all_faces.size() == topo.num_sides());
 
     std::vector<std::vector<double>> nodeLocations = {
       {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
diff --git a/unit_tests/matrix_free/UnitTestStkSimdFaceConnectivityMap.C b/unit_tests/matrix_free/UnitTestStkSimdFaceConnectivityMap.C
index c79c020f2..bbd35e0f1 100644
--- a/unit_tests/matrix_free/UnitTestStkSimdFaceConnectivityMap.C
+++ b/unit_tests/matrix_free/UnitTestStkSimdFaceConnectivityMap.C
@@ -107,7 +107,7 @@ protected:
     stk::mesh::get_selected_entities(
       surfaceSelector, bulk.get_buckets(meta.side_rank(), surfaceSelector),
       all_faces);
-    ThrowRequire(all_faces.size() == topo.num_sides());
+    STK_ThrowRequire(all_faces.size() == topo.num_sides());
 
     std::vector<std::vector<double>> nodeLocations = {
       {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
diff --git a/unit_tests/matrix_free/UnitTestStkSimdNodeConnectivityMap.C b/unit_tests/matrix_free/UnitTestStkSimdNodeConnectivityMap.C
index 0f7a4d1f3..510694c1b 100644
--- a/unit_tests/matrix_free/UnitTestStkSimdNodeConnectivityMap.C
+++ b/unit_tests/matrix_free/UnitTestStkSimdNodeConnectivityMap.C
@@ -91,7 +91,7 @@ protected:
     stk::mesh::get_selected_entities(
       surfaceSelector, bulk.get_buckets(meta.side_rank(), surfaceSelector),
       all_faces);
-    ThrowRequire(all_faces.size() == topo.num_sides());
+    STK_ThrowRequire(all_faces.size() == topo.num_sides());
 
     std::vector<std::vector<double>> nodeLocations = {
       {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
diff --git a/unit_tests/ngp_algorithms/UnitTestNgpAlgUtils.C b/unit_tests/ngp_algorithms/UnitTestNgpAlgUtils.C
index bfeaf94f1..f28c1fad9 100644
--- a/unit_tests/ngp_algorithms/UnitTestNgpAlgUtils.C
+++ b/unit_tests/ngp_algorithms/UnitTestNgpAlgUtils.C
@@ -48,7 +48,7 @@ linear_scalar_field(
       fieldData[2] = coordsData[2] * zCoeff;
     }
   } else {
-    ThrowErrorMsg(
+    STK_ThrowErrorMsg(
       "linear_scalar_field(): Field has unhandled length " << fieldLength);
   }
 
diff --git a/unit_tests/ngp_kernels/UnitTestNgpLoops.C b/unit_tests/ngp_kernels/UnitTestNgpLoops.C
index d649cceb4..f29fba07d 100644
--- a/unit_tests/ngp_kernels/UnitTestNgpLoops.C
+++ b/unit_tests/ngp_kernels/UnitTestNgpLoops.C
@@ -55,7 +55,7 @@ public:
       *velocity, stk::io::FieldOutputType::VECTOR_3D);
     stk::mesh::put_field_on_mesh(
       *massFlowRate, meta->universal_part(), hex8SCS.num_integration_points(),
-      &zero);
+      nullptr);
     stk::mesh::put_field_on_mesh(*mdotEdge, meta->universal_part(), &zero);
   }
 
diff --git a/wind-utils b/wind-utils
index fc9710c41..a055d9438 160000
--- a/wind-utils
+++ b/wind-utils
@@ -1 +1 @@
-Subproject commit fc9710c419111f5476d50586fd1a8e12db52b713
+Subproject commit a055d9438f1348f8bdeb5d320343c83c469de220