diff --git a/include/ArrayND.h b/include/ArrayND.h
new file mode 100644
index 0000000000..75ba6c3606
--- /dev/null
+++ b/include/ArrayND.h
@@ -0,0 +1,180 @@
+#ifndef ARRAY_ND_H
+#define ARRAY_ND_H
+
+#include "Kokkos_Macros.hpp"
+#include <type_traits>
+
+namespace sierra::nalu {
+
+// stack array set to interoperate with Kokkos views
+template <typename ArrayType, typename = void>
+struct ArrayND
+{
+};
+
+template <typename ArrayType, decltype(std::rank_v<ArrayType>) r>
+using enable_if_rank = std::enable_if_t<std::rank_v<ArrayType> == r>;
+
+template <typename ArrayType>
+struct ArrayND<ArrayType, enable_if_rank<ArrayType, 1>>
+{
+  KOKKOS_DEFAULTED_FUNCTION constexpr ArrayND() = default;
+
+  static constexpr int rank = 1;
+  using value_type = std::remove_all_extents_t<ArrayType>;
+  value_type internal_data_[std::extent_v<ArrayType, 0>];
+
+  [[nodiscard]] static constexpr int extent_int(int /*unused*/)
+  {
+    return int(std::extent_v<ArrayType, 0>);
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
+  operator()(int i) const noexcept
+  {
+    return internal_data_[i];
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
+  operator[](int i) const noexcept
+  {
+    return internal_data_[i];
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type&
+  operator()(int i) noexcept
+  {
+    return internal_data_[i];
+  }
+};
+
+template <typename ArrayType>
+struct ArrayND<ArrayType, enable_if_rank<ArrayType, 2>>
+{
+  KOKKOS_DEFAULTED_FUNCTION constexpr ArrayND() = default;
+
+  static constexpr int rank = 2;
+
+  using value_type = std::remove_all_extents_t<ArrayType>;
+  value_type internal_data_[std::extent_v<ArrayType, 0>]
+                           [std::extent_v<ArrayType, 1>];
+
+  [[nodiscard]] static constexpr int extent_int(int n)
+  {
+    return (n == 0) ? int(std::extent_v<ArrayType, 0>)
+                    : int(std::extent_v<ArrayType, 1>);
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
+  operator()(int j, int i) const noexcept
+  {
+    return internal_data_[j][i];
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type&
+  operator()(int j, int i) noexcept
+  {
+    return internal_data_[j][i];
+  }
+};
+
+template <typename ArrayType>
+struct ArrayND<ArrayType, enable_if_rank<ArrayType, 3>>
+{
+  KOKKOS_DEFAULTED_FUNCTION constexpr ArrayND() = default;
+
+  static constexpr int rank = 3;
+  [[nodiscard]] static constexpr int extent_int(int n)
+  {
+    return (n == 0)   ? int(std::extent_v<ArrayType, 0>)
+           : (n == 1) ? int(std::extent_v<ArrayType, 1>)
+                      : int(std::extent_v<ArrayType, 2>);
+  }
+
+  using value_type = std::remove_all_extents_t<ArrayType>;
+  value_type internal_data_[std::extent_v<ArrayType, 0>]
+                           [std::extent_v<ArrayType, 1>]
+                           [std::extent_v<ArrayType, 2>];
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr auto
+  operator()(int k, int j, int i) const noexcept
+  {
+    return internal_data_[k][j][i];
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr auto&
+  operator()(int k, int j, int i) noexcept
+  {
+    return internal_data_[k][j][i];
+  }
+};
+
+template <typename ArrayType>
+struct ArrayND<ArrayType, enable_if_rank<ArrayType, 4>>
+{
+  KOKKOS_DEFAULTED_FUNCTION constexpr ArrayND() = default;
+
+  static constexpr int rank = 4;
+  [[nodiscard]] static constexpr int extent_int(int n)
+  {
+    return (n == 0)   ? int(std::extent_v<ArrayType, 0>)
+           : (n == 1) ? int(std::extent_v<ArrayType, 1>)
+           : (n == 2) ? int(std::extent_v<ArrayType, 2>)
+                      : int(std::extent_v<ArrayType, 3>);
+  }
+
+  using value_type = std::remove_all_extents_t<ArrayType>;
+  value_type
+    internal_data_[std::extent_v<ArrayType, 0>][std::extent_v<ArrayType, 1>]
+                  [std::extent_v<ArrayType, 2>][std::extent_v<ArrayType, 3>];
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
+  operator()(int l, int k, int j, int i) const noexcept
+  {
+    return internal_data_[l][k][j][i];
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type&
+  operator()(int l, int k, int j, int i) noexcept
+  {
+    return internal_data_[l][k][j][i];
+  }
+};
+
+template <typename ArrayType>
+struct ArrayND<ArrayType, enable_if_rank<ArrayType, 5>>
+{
+  KOKKOS_DEFAULTED_FUNCTION constexpr ArrayND() = default;
+
+  static constexpr int rank = 5;
+  [[nodiscard]] static constexpr int extent_int(int n)
+  {
+    return (n == 0)   ? int(std::extent_v<ArrayType, 0>)
+           : (n == 1) ? int(std::extent_v<ArrayType, 1>)
+           : (n == 2) ? int(std::extent_v<ArrayType, 2>)
+           : (n == 3) ? int(std::extent_v<ArrayType, 3>)
+                      : int(std::extent_v<ArrayType, 4>);
+  }
+
+  using value_type = std::remove_all_extents_t<ArrayType>;
+  value_type
+    internal_data_[std::extent_v<ArrayType, 0>][std::extent_v<ArrayType, 1>]
+                  [std::extent_v<ArrayType, 2>][std::extent_v<ArrayType, 3>]
+                  [std::extent_v<ArrayType, 4>];
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
+  operator()(int m, int l, int k, int j, int i) const noexcept
+  {
+    return internal_data_[m][l][k][j][i];
+  }
+
+  [[nodiscard]] KOKKOS_FORCEINLINE_FUNCTION constexpr value_type&
+  operator()(int m, int l, int k, int j, int i) noexcept
+  {
+    return internal_data_[m][l][k][j][i];
+  }
+};
+
+} // namespace sierra::nalu
+
+#endif
\ No newline at end of file
diff --git a/include/matrix_free/Coefficients.h b/include/matrix_free/Coefficients.h
index 3ff31c5572..d51d900c05 100644
--- a/include/matrix_free/Coefficients.h
+++ b/include/matrix_free/Coefficients.h
@@ -10,7 +10,7 @@
 #ifndef CVFEM_COEFFICIENTS_H
 #define CVFEM_COEFFICIENTS_H
 
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 // An option to use a symmetric version of
 // the P=2 CVFEM operator.  This loses one to two orders
@@ -32,10 +32,10 @@ struct Coeffs
 template <>
 struct Coeffs<1>
 {
-  using nodal_matrix_type = LocalArray<double[2][2]>;
-  using scs_matrix_type = LocalArray<double[1][2]>;
-  using linear_nodal_matrix_type = LocalArray<double[2][2]>;
-  using linear_scs_matrix_type = LocalArray<double[2][1]>;
+  using nodal_matrix_type = ArrayND<double[2][2]>;
+  using scs_matrix_type = ArrayND<double[1][2]>;
+  using linear_nodal_matrix_type = ArrayND<double[2][2]>;
+  using linear_scs_matrix_type = ArrayND<double[2][1]>;
   static constexpr nodal_matrix_type W = {{{0.75, 0.25}, {0.25, 0.75}}};
 
   static constexpr nodal_matrix_type D = {{{-0.5, +0.5}, {-0.5, 0.5}}};
@@ -46,17 +46,17 @@ struct Coeffs<1>
   static constexpr linear_nodal_matrix_type Nlin = {{{1, 0}, {0, 1}}};
   static constexpr linear_scs_matrix_type Ntlin = {{{0.5}, {0.5}}};
 
-  static constexpr LocalArray<double[2]> Wl = {{1, 1}};
+  static constexpr ArrayND<double[2]> Wl = {{1, 1}};
 };
 
 #if USE_SYMMETRIC
 template <>
 struct Coeffs<2>
 {
-  using nodal_matrix_type = LocalArray<double[3][3]>;
-  using scs_matrix_type = LocalArray<double[2][3]>;
-  using linear_nodal_matrix_type = LocalArray<double[2][3]>;
-  using linear_scs_matrix_type = LocalArray<double[2][2]>;
+  using nodal_matrix_type = ArrayND<double[3][3]>;
+  using scs_matrix_type = ArrayND<double[2][3]>;
+  using linear_nodal_matrix_type = ArrayND<double[2][3]>;
+  using linear_scs_matrix_type = ArrayND<double[2][2]>;
 
   static constexpr nodal_matrix_type W = {
     {{0.2561728395061728395062, 0.09876543209876543209877,
@@ -87,7 +87,7 @@ struct Coeffs<2>
     {{0.8333333333333333333333, 0.1666666666666666666667},
      {0.1666666666666666666667, 0.8333333333333333333333}}};
 
-  static constexpr LocalArray<double[3]> Wl = {
+  static constexpr ArrayND<double[3]> Wl = {
     {0.3333333333333333333333, 1.333333333333333333333,
      0.3333333333333333333333}};
 };
@@ -95,10 +95,10 @@ struct Coeffs<2>
 template <>
 struct Coeffs<2>
 {
-  using nodal_matrix_type = LocalArray<double[3][3]>;
-  using scs_matrix_type = LocalArray<double[2][3]>;
-  using linear_nodal_matrix_type = LocalArray<double[2][3]>;
-  using linear_scs_matrix_type = LocalArray<double[2][2]>;
+  using nodal_matrix_type = ArrayND<double[3][3]>;
+  using scs_matrix_type = ArrayND<double[2][3]>;
+  using linear_nodal_matrix_type = ArrayND<double[2][3]>;
+  using linear_scs_matrix_type = ArrayND<double[2][2]>;
 
   static constexpr nodal_matrix_type W = {{
     {+0.301258318378354124194, +0.15346642738699932044,
@@ -129,7 +129,7 @@ struct Coeffs<2>
     {{+0.7886751345948129, +0.21132486540518708},
      {+0.21132486540518708, +0.7886751345948129}}};
 
-  static constexpr LocalArray<double[3]> Wl = {
+  static constexpr ArrayND<double[3]> Wl = {
     {0.422649730810374235491, 1.154700538379251529018,
      0.422649730810374235491}};
 };
@@ -137,10 +137,10 @@ struct Coeffs<2>
 template <>
 struct Coeffs<3>
 {
-  using nodal_matrix_type = LocalArray<double[4][4]>;
-  using scs_matrix_type = LocalArray<double[3][4]>;
-  using linear_nodal_matrix_type = LocalArray<double[2][4]>;
-  using linear_scs_matrix_type = LocalArray<double[2][3]>;
+  using nodal_matrix_type = ArrayND<double[4][4]>;
+  using scs_matrix_type = ArrayND<double[3][4]>;
+  using linear_nodal_matrix_type = ArrayND<double[2][4]>;
+  using linear_scs_matrix_type = ArrayND<double[2][3]>;
 
   static constexpr nodal_matrix_type W = {
     {{0.1583333333333333333333, 0.08527003148341972055897,
@@ -190,7 +190,7 @@ struct Coeffs<3>
      {0.1127016653792582978610, 0.5000000000000000000000,
       0.8872983346207417021390}}};
 
-  static constexpr LocalArray<double[4]> Wl = {
+  static constexpr ArrayND<double[4]> Wl = {
     {0.225403330758516622964, 0.774596669241483377036, 0.774596669241483377036,
      0.225403330758516622964}};
 };
@@ -198,10 +198,10 @@ struct Coeffs<3>
 template <>
 struct Coeffs<4>
 {
-  using nodal_matrix_type = LocalArray<double[5][5]>;
-  using scs_matrix_type = LocalArray<double[4][5]>;
-  using linear_nodal_matrix_type = LocalArray<double[2][5]>;
-  using linear_scs_matrix_type = LocalArray<double[2][4]>;
+  using nodal_matrix_type = ArrayND<double[5][5]>;
+  using scs_matrix_type = ArrayND<double[4][5]>;
+  using linear_nodal_matrix_type = ArrayND<double[2][5]>;
+  using linear_scs_matrix_type = ArrayND<double[2][4]>;
 
   static constexpr nodal_matrix_type W = {
     {{0.09695253015044793809126, 0.05331486033726615325831,
@@ -268,7 +268,7 @@ struct Coeffs<4>
      {0.06943184420297371373110, 0.3300094782075718713443,
       0.6699905217924281286557, 0.9305681557970262307578}}};
 
-  static constexpr LocalArray<double[5]> Wl = {
+  static constexpr ArrayND<double[5]> Wl = {
     {0.13886368840594742478, 0.52115526800919631042, 0.679962087169712529605,
      0.52115526800919631042, 0.13886368840594742478}};
 };
diff --git a/include/matrix_free/ConductionInterior.h b/include/matrix_free/ConductionInterior.h
index 445140eeeb..8a48c91804 100644
--- a/include/matrix_free/ConductionInterior.h
+++ b/include/matrix_free/ConductionInterior.h
@@ -12,7 +12,7 @@
 
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include "Kokkos_Array.hpp"
 #include "Tpetra_MultiVector.hpp"
@@ -29,7 +29,7 @@ namespace impl {
 template <int p>
 struct conduction_residual_t
 {
-  using narray = LocalArray<ftype[p + 1][p + 1][p + 1]>;
+  using narray = ArrayND<ftype[p + 1][p + 1][p + 1]>;
 
   static void invoke(
     Kokkos::Array<double, 3> gammas,
@@ -47,7 +47,7 @@ namespace impl {
 template <int p>
 struct conduction_linearized_residual_t
 {
-  using narray = LocalArray<ftype[p + 1][p + 1][p + 1]>;
+  using narray = ArrayND<ftype[p + 1][p + 1][p + 1]>;
 
   static void invoke(
     double gamma,
diff --git a/include/matrix_free/ContinuityInterior.h b/include/matrix_free/ContinuityInterior.h
index e93efa4e47..44a695af3f 100644
--- a/include/matrix_free/ContinuityInterior.h
+++ b/include/matrix_free/ContinuityInterior.h
@@ -11,7 +11,7 @@
 #define CONTINUITY_INTERIOR_H
 
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 
 #include "Teuchos_RCP.hpp"
@@ -29,7 +29,7 @@ namespace impl {
 template <int p>
 struct continuity_residual_t
 {
-  using narray = LocalArray<ftype[p + 1][p + 1][p + 1]>;
+  using narray = ArrayND<ftype[p + 1][p + 1][p + 1]>;
 
   static void invoke(
     double scaling,
@@ -44,7 +44,7 @@ namespace impl {
 template <int p>
 struct continuity_linearized_residual_t
 {
-  using narray = LocalArray<ftype[p + 1][p + 1][p + 1]>;
+  using narray = ArrayND<ftype[p + 1][p + 1][p + 1]>;
 
   static void invoke(
     const_elem_offset_view<p> offsets,
diff --git a/include/matrix_free/ElementFluxIntegral.h b/include/matrix_free/ElementFluxIntegral.h
index 6713d0381e..888a3ed53e 100644
--- a/include/matrix_free/ElementFluxIntegral.h
+++ b/include/matrix_free/ElementFluxIntegral.h
@@ -13,7 +13,7 @@
 #include <cmath>
 
 #include "matrix_free/KokkosFramework.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/ShuffledAccess.h"
 #include "matrix_free/ElementSCSInterpolate.h"
 
@@ -88,7 +88,7 @@ KOKKOS_FORCEINLINE_FUNCTION void
 scalar_flux_divergence(int index, const InArray& in, OutArray& out)
 {
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[p + 1][p + 1]> scratch;
+    ArrayND<ftype[p + 1][p + 1]> scratch;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         ftype acc = 0;
@@ -131,7 +131,7 @@ advdiff_flux(
 {
   for (int l = 0; l < p; ++l) {
     enum { LEVEL_0 = 0, LEVEL_1 = 1 };
-    LocalArray<ftype[2][p + 1][p + 1]> scratch;
+    ArrayND<ftype[2][p + 1][p + 1]> scratch;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         ftype acc = 0;
@@ -222,7 +222,7 @@ diffusive_flux(
 {
   enum { LEVEL_0 = 0, LEVEL_1 = 1 };
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[2][p + 1][p + 1]> scratch;
+    ArrayND<ftype[2][p + 1][p + 1]> scratch;
 
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
@@ -313,7 +313,7 @@ scalar_flux_vector(
 {
   enum { LEVEL_0 = 0, LEVEL_1 = 1 };
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[2][p + 1][p + 1]> scratch;
+    ArrayND<ftype[2][p + 1][p + 1]> scratch;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         ftype in_scs(0);
diff --git a/include/matrix_free/ElementGradient.h b/include/matrix_free/ElementGradient.h
index ac974e6111..aa02df52e2 100644
--- a/include/matrix_free/ElementGradient.h
+++ b/include/matrix_free/ElementGradient.h
@@ -15,7 +15,7 @@
 #include "matrix_free/ElementSCSInterpolate.h"
 #include "matrix_free/GeometricFunctions.h"
 #include "matrix_free/KokkosFramework.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/ShuffledAccess.h"
 #include "matrix_free/TensorOperations.h"
 
@@ -25,10 +25,10 @@ namespace matrix_free {
 
 template <int p, typename BoxArray, typename InpArray>
 KOKKOS_FUNCTION
-  typename std::enable_if<InpArray::Rank == 4, LocalArray<ftype[3][3]>>::type
+  typename std::enable_if<InpArray::rank == 4, ArrayND<ftype[3][3]>>::type
   gradient_nodal(const BoxArray& box, const InpArray& u, int k, int j, int i)
 {
-  LocalArray<ftype[3][3]> gu_hat;
+  ArrayND<ftype[3][3]> gu_hat;
   for (int dj = 0; dj < 3; ++dj) {
     gu_hat(dj, 0) = 0;
     gu_hat(dj, 1) = 0;
@@ -41,17 +41,17 @@ KOKKOS_FUNCTION
     }
   }
 
-  LocalArray<ftype[3][3]> gu;
+  ArrayND<ftype[3][3]> gu;
   inv_transform_t(geom::linear_hex_jacobian<p>(box, k, j, i), gu_hat, gu);
   return gu;
 }
 
 template <int p, typename BoxArray, typename InpArray>
 KOKKOS_FUNCTION
-  typename std::enable_if<InpArray::Rank == 3, LocalArray<ftype[3]>>::type
+  typename std::enable_if<InpArray::rank == 3, ArrayND<ftype[3]>>::type
   gradient_nodal(const BoxArray& box, const InpArray& u, int k, int j, int i)
 {
-  LocalArray<ftype[3]> gu_hat;
+  ArrayND<ftype[3]> gu_hat;
   gu_hat(0) = 0;
   gu_hat(1) = 0;
   gu_hat(2) = 0;
@@ -61,7 +61,7 @@ KOKKOS_FUNCTION
     gu_hat(1) += D(j, q) * u(k, q, i);
     gu_hat(2) += D(k, q) * u(q, j, i);
   }
-  LocalArray<ftype[3]> gu;
+  ArrayND<ftype[3]> gu;
   inv_transform_t(geom::linear_hex_jacobian<p>(box, k, j, i), gu_hat, gu);
   return gu;
 }
@@ -73,7 +73,7 @@ template <
   typename UArrayType,
   typename UHatArrayType>
 KOKKOS_FORCEINLINE_FUNCTION
-  typename std::enable_if<UArrayType::Rank == 4, LocalArray<ftype[3][3]>>::type
+  typename std::enable_if<UArrayType::rank == 4, ArrayND<ftype[3][3]>>::type
   gradient_scs(
     const BoxArrayType& box,
     const UArrayType& u,
@@ -87,7 +87,7 @@ KOKKOS_FORCEINLINE_FUNCTION
   constexpr int dir_1 = (dir == XH) ? YH : XH;
   constexpr int dir_2 = (dir == ZH) ? YH : ZH;
 
-  LocalArray<ftype[3][3]> gu_hat;
+  ArrayND<ftype[3][3]> gu_hat;
   for (int d = 0; d < 3; ++d) {
     ftype acc = 0;
     for (int q = 0; q < p + 1; ++q) {
@@ -110,7 +110,7 @@ KOKKOS_FORCEINLINE_FUNCTION
     }
     gu_hat(d, dir_2) = acc;
   }
-  LocalArray<ftype[3][3]> gu;
+  ArrayND<ftype[3][3]> gu;
   inv_transform_t(
     geom::linear_hex_jacobian_scs<p, dir>(box, l, s, r), gu_hat, gu);
   return gu;
@@ -123,7 +123,7 @@ template <
   typename UArrayType,
   typename UHatArrayType>
 KOKKOS_FORCEINLINE_FUNCTION
-  typename std::enable_if<UArrayType::Rank == 3, LocalArray<ftype[3]>>::type
+  typename std::enable_if<UArrayType::rank == 3, ArrayND<ftype[3]>>::type
   gradient_scs(
     const BoxArrayType& box,
     const UArrayType& u,
@@ -137,7 +137,7 @@ KOKKOS_FORCEINLINE_FUNCTION
   constexpr int dir_1 = (dir == XH) ? YH : XH;
   constexpr int dir_2 = (dir == ZH) ? YH : ZH;
 
-  LocalArray<ftype[3]> gu_hat;
+  ArrayND<ftype[3]> gu_hat;
   for (int d = 0; d < 3; ++d) {
     ftype acc = 0;
     for (int q = 0; q < p + 1; ++q) {
@@ -160,7 +160,7 @@ KOKKOS_FORCEINLINE_FUNCTION
     }
     gu_hat(dir_2) = acc;
   }
-  LocalArray<ftype[3]> gu;
+  ArrayND<ftype[3]> gu;
   inv_transform_t(
     geom::linear_hex_jacobian_scs<p, dir>(box, l, s, r), gu_hat, gu);
   return gu;
diff --git a/include/matrix_free/ElementSCSInterpolate.h b/include/matrix_free/ElementSCSInterpolate.h
index 906e0f9167..2bf9c46a32 100644
--- a/include/matrix_free/ElementSCSInterpolate.h
+++ b/include/matrix_free/ElementSCSInterpolate.h
@@ -14,7 +14,7 @@
 
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/ShuffledAccess.h"
 
 namespace sierra {
@@ -62,7 +62,7 @@ grad_scs(const ScalarView& phi, int l, int s, int r)
     return val;
   }
   case dir_1: {
-    LocalArray<ftype[p + 1]> val_array;
+    ArrayND<ftype[p + 1]> val_array;
     for (int q = 0; q < p + 1; ++q) {
       val_array(q) = interp_scs<p, dir>(phi, l, s, q);
     }
@@ -75,7 +75,7 @@ grad_scs(const ScalarView& phi, int l, int s, int r)
     return val;
   }
   default: {
-    LocalArray<ftype[p + 1]> val_array;
+    ArrayND<ftype[p + 1]> val_array;
     for (int q = 0; q < p + 1; ++q) {
       val_array(q) = interp_scs<p, dir>(phi, l, q, r);
     }
@@ -108,7 +108,7 @@ grad_scs(const ScalarView& phi, int l, int s, int r, int d)
   }
 
   if (dk == dir_1) {
-    LocalArray<ftype[p + 1]> val_array;
+    ArrayND<ftype[p + 1]> val_array;
     for (int q = 0; q < p + 1; ++q) {
       val_array(q) = interp_scs<p, dir>(phi, l, s, q, d);
     }
@@ -122,7 +122,7 @@ grad_scs(const ScalarView& phi, int l, int s, int r, int d)
   }
 
   if (dk == dir_2) {
-    LocalArray<ftype[p + 1]> val_array;
+    ArrayND<ftype[p + 1]> val_array;
     for (int q = 0; q < p + 1; ++q) {
       val_array(q) = interp_scs<p, dir>(phi, l, q, r, d);
     }
diff --git a/include/matrix_free/ElementVolumeIntegral.h b/include/matrix_free/ElementVolumeIntegral.h
index e52e5cc4f9..1a9ab9be50 100644
--- a/include/matrix_free/ElementVolumeIntegral.h
+++ b/include/matrix_free/ElementVolumeIntegral.h
@@ -14,7 +14,7 @@
 
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/ShuffledAccess.h"
 
 namespace sierra {
@@ -35,7 +35,7 @@ apply_mass(
   static constexpr auto vandermonde = Coeffs<p>::W;
   for (int k = 0; k < p + 1; ++k) {
     for (int j = 0; j < p + 1; ++j) {
-      LocalArray<ftype[p + 1]> scratch;
+      ArrayND<ftype[p + 1]> scratch;
       for (int i = 0; i < p + 1; ++i) {
         scratch(i) = -gamma * vol(index, k, j, i) * delta(k, j, i);
       }
@@ -50,7 +50,7 @@ apply_mass(
   }
 
   for (int i = 0; i < p + 1; ++i) {
-    LocalArray<ftype[p + 1][p + 1]> scratch;
+    ArrayND<ftype[p + 1][p + 1]> scratch;
     for (int k = 0; k < p + 1; ++k) {
       for (int j = 0; j < p + 1; ++j) {
         ftype acc(0);
@@ -118,7 +118,7 @@ consistent_mass_time_derivative(
 
   for (int k = 0; k < p + 1; ++k) {
     for (int j = 0; j < p + 1; ++j) {
-      LocalArray<ftype[p + 1]> scratch_1d;
+      ArrayND<ftype[p + 1]> scratch_1d;
       for (int i = 0; i < p + 1; ++i) {
         scratch_1d(i) =
           -vol(index, k, j, i) *
@@ -137,7 +137,7 @@ consistent_mass_time_derivative(
   }
 
   for (int i = 0; i < p + 1; ++i) {
-    LocalArray<ftype[p + 1][p + 1]> scratch;
+    ArrayND<ftype[p + 1][p + 1]> scratch;
     for (int k = 0; k < p + 1; ++k) {
       for (int j = 0; j < p + 1; ++j) {
         scratch(k, j) = 0;
diff --git a/include/matrix_free/FilterDiagonal.h b/include/matrix_free/FilterDiagonal.h
index c8c3f59cf7..a3cfb5271a 100644
--- a/include/matrix_free/FilterDiagonal.h
+++ b/include/matrix_free/FilterDiagonal.h
@@ -13,7 +13,7 @@
 #include "Tpetra_MultiVector.hpp"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
diff --git a/include/matrix_free/GeometricFunctions.h b/include/matrix_free/GeometricFunctions.h
index b72899b936..3b7aa07e62 100644
--- a/include/matrix_free/GeometricFunctions.h
+++ b/include/matrix_free/GeometricFunctions.h
@@ -81,7 +81,7 @@ hex_jacobian_component_scs(const BoxArray& box, int l, int s, int r)
 }
 
 template <int p, int dk, typename BoxArray>
-KOKKOS_FUNCTION LocalArray<typename BoxArray::value_type[3][3]>
+KOKKOS_FUNCTION ArrayND<typename BoxArray::value_type[3][3]>
 linear_hex_jacobian_scs(const BoxArray& box, int k, int j, int i)
 {
   enum { XH = 0, YH = 1, ZH = 2 };
@@ -103,14 +103,14 @@ linear_hex_jacobian_scs(const BoxArray& box, int k, int j, int i)
 }
 
 template <int p, int dk, typename BoxArray>
-KOKKOS_FUNCTION LocalArray<typename BoxArray::value_type[3][3]>
+KOKKOS_FUNCTION ArrayND<typename BoxArray::value_type[3][3]>
 linear_hex_invjact_scs(const BoxArray& box, int k, int j, int i)
 {
   return invert_transpose_matrix(linear_hex_jacobian_scs<p, dk>(box, k, j, i));
 }
 
 template <int p, int dk, typename BoxArrayType>
-KOKKOS_FUNCTION LocalArray<ftype[3]>
+KOKKOS_FUNCTION ArrayND<ftype[3]>
 linear_area(const BoxArrayType& box, int k, int j, int i)
 {
   enum { XH = 0, YH = 1, ZH = 2 };
@@ -122,13 +122,13 @@ linear_area(const BoxArrayType& box, int k, int j, int i)
   const auto dy_ds2 = hex_jacobian_component_scs<p, dk, ds2, YH>(box, k, j, i);
   const auto dz_ds1 = hex_jacobian_component_scs<p, dk, ds1, ZH>(box, k, j, i);
   const auto dz_ds2 = hex_jacobian_component_scs<p, dk, ds2, ZH>(box, k, j, i);
-  return LocalArray<ftype[3]>{
+  return ArrayND<ftype[3]>{
     {dy_ds1 * dz_ds2 - dz_ds1 * dy_ds2, dz_ds1 * dx_ds2 - dx_ds1 * dz_ds2,
      dx_ds1 * dy_ds2 - dy_ds1 * dx_ds2}};
 }
 
 template <int p, int dk, typename BoxArray>
-KOKKOS_FUNCTION LocalArray<ftype[3]>
+KOKKOS_FUNCTION ArrayND<ftype[3]>
 laplacian_metric(const BoxArray& box, int k, int j, int i)
 {
   enum { XH = 0, YH = 1, ZH = 2 };
@@ -141,7 +141,7 @@ laplacian_metric(const BoxArray& box, int k, int j, int i)
 
   switch (dk) {
   case XH: {
-    return LocalArray<ftype[3]>{
+    return ArrayND<ftype[3]>{
       {-inv_detj * (adj_jac(XH, XH) * adj_jac(XH, XH) +
                     adj_jac(XH, YH) * adj_jac(XH, YH) +
                     adj_jac(XH, ZH) * adj_jac(XH, ZH)),
@@ -153,7 +153,7 @@ laplacian_metric(const BoxArray& box, int k, int j, int i)
                     adj_jac(XH, ZH) * adj_jac(ZH, ZH))}};
   }
   case YH: {
-    return LocalArray<ftype[3]>{
+    return ArrayND<ftype[3]>{
       {-inv_detj * (adj_jac(YH, XH) * adj_jac(YH, XH) +
                     adj_jac(YH, YH) * adj_jac(YH, YH) +
                     adj_jac(YH, ZH) * adj_jac(YH, ZH)),
@@ -165,7 +165,7 @@ laplacian_metric(const BoxArray& box, int k, int j, int i)
                     adj_jac(YH, ZH) * adj_jac(ZH, ZH))}};
   }
   default:
-    return LocalArray<ftype[3]>{
+    return ArrayND<ftype[3]>{
       {-inv_detj * (adj_jac(ZH, XH) * adj_jac(ZH, XH) +
                     adj_jac(ZH, YH) * adj_jac(ZH, YH) +
                     adj_jac(ZH, ZH) * adj_jac(ZH, ZH)),
@@ -219,12 +219,12 @@ hex_jacobian_component(
 }
 
 template <int p, typename BoxArray>
-KOKKOS_FUNCTION LocalArray<typename BoxArray::value_type[3][3]>
+KOKKOS_FUNCTION ArrayND<typename BoxArray::value_type[3][3]>
 linear_hex_jacobian(const BoxArray& box, int k, int j, int i)
 {
   static constexpr auto coeff = Coeffs<p>::Nlin;
   enum { XH = 0, YH = 1, ZH = 2 };
-  LocalArray<typename BoxArray::value_type[3][3]> jac;
+  ArrayND<typename BoxArray::value_type[3][3]> jac;
   jac(0, 0) = hex_jacobian_component<p, XH, XH>(coeff, box, k, j, i);
   jac(0, 1) = hex_jacobian_component<p, XH, YH>(coeff, box, k, j, i);
   jac(0, 2) = hex_jacobian_component<p, XH, ZH>(coeff, box, k, j, i);
diff --git a/include/matrix_free/HexVertexCoordinates.h b/include/matrix_free/HexVertexCoordinates.h
index 0a1a254bca..2d154c4175 100644
--- a/include/matrix_free/HexVertexCoordinates.h
+++ b/include/matrix_free/HexVertexCoordinates.h
@@ -10,7 +10,7 @@
 #ifndef HEX_VERTEX_COORDINATES_H
 #define HEX_VERTEX_COORDINATES_H
 
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/KokkosViewTypes.h"
 
@@ -19,11 +19,11 @@ namespace nalu {
 namespace matrix_free {
 
 template <int p, typename ElemCoordsArray>
-KOKKOS_FUNCTION LocalArray<ftype[3][8]>
+KOKKOS_FUNCTION ArrayND<ftype[3][8]>
 hex_vertex_coordinates(const ElemCoordsArray& xc)
 {
-  static_assert(ElemCoordsArray::Rank == 4, "");
-  LocalArray<ftype[3][8]> box;
+  static_assert(ElemCoordsArray::rank == 4, "");
+  ArrayND<ftype[3][8]> box;
   for (int d = 0; d < 3; ++d) {
     box(d, 0) = xc(0, 0, 0, d);
     box(d, 1) = xc(0, 0, p, d);
@@ -38,11 +38,11 @@ hex_vertex_coordinates(const ElemCoordsArray& xc)
 }
 
 template <int p, typename ElemCoordsArray>
-KOKKOS_FUNCTION LocalArray<ftype[3][8]>
+KOKKOS_FUNCTION ArrayND<ftype[3][8]>
 hex_vertex_coordinates(int index, const ElemCoordsArray& xc)
 {
-  static_assert(ElemCoordsArray::Rank == 5, "");
-  LocalArray<ftype[3][8]> box;
+  static_assert(ElemCoordsArray::rank == 5, "");
+  ArrayND<ftype[3][8]> box;
   for (int d = 0; d < 3; ++d) {
     box(d, 0) = xc(index, 0, 0, 0, d);
     box(d, 1) = xc(index, 0, 0, p, d);
@@ -57,11 +57,11 @@ hex_vertex_coordinates(int index, const ElemCoordsArray& xc)
 }
 
 template <typename ElemCoordsArray>
-KOKKOS_FUNCTION LocalArray<ftype[3][8]>
+KOKKOS_FUNCTION ArrayND<ftype[3][8]>
 hex_vertex_coordinates(int n, int m, int l, const ElemCoordsArray& xc)
 {
-  static_assert(ElemCoordsArray::Rank == 4, "");
-  LocalArray<ftype[3][8]> box;
+  static_assert(ElemCoordsArray::rank == 4, "");
+  ArrayND<ftype[3][8]> box;
   for (int d = 0; d < 3; ++d) {
     box(d, 0) = xc(n + 0, m + 0, l + 0, d);
     box(d, 1) = xc(n + 0, m + 0, l + 1, d);
@@ -76,10 +76,10 @@ hex_vertex_coordinates(int n, int m, int l, const ElemCoordsArray& xc)
 }
 
 template <int p>
-KOKKOS_FUNCTION LocalArray<ftype[3][4]>
+KOKKOS_FUNCTION ArrayND<ftype[3][4]>
 face_vertex_coordinates(int index, const const_face_vector_view<p>& xc)
 {
-  LocalArray<ftype[3][4]> base_box;
+  ArrayND<ftype[3][4]> base_box;
   for (int d = 0; d < 3; ++d) {
     base_box(d, 0) = xc(index, 0, 0, d);
     base_box(d, 1) = xc(index, 0, p, d);
diff --git a/include/matrix_free/LocalArray.h b/include/matrix_free/LocalArray.h
deleted file mode 100644
index bd81bb038f..0000000000
--- a/include/matrix_free/LocalArray.h
+++ /dev/null
@@ -1,182 +0,0 @@
-// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
-// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
-// Northwest Research Associates. Under the terms of Contract DE-NA0003525
-// with NTESS, the U.S. Government retains certain rights in this software.
-//
-// This software is released under the BSD 3-clause license. See LICENSE file
-// for more details.
-//
-
-#ifndef LOCALARRAY_H
-#define LOCALARRAY_H
-
-#include "Kokkos_Macros.hpp"
-#include "matrix_free/KokkosFramework.h"
-#include <type_traits>
-
-namespace sierra {
-namespace nalu {
-namespace matrix_free {
-
-template <typename ArrayType, typename = void>
-struct LocalArray
-{
-};
-
-template <typename ArrayType>
-struct alignas(alignment) LocalArray<
-  ArrayType,
-  typename std::enable_if<std::rank<ArrayType>::value == 1>::type>
-{
-  static constexpr int Rank = 1;
-  using value_type = typename std::remove_all_extents<ArrayType>::type;
-  static constexpr int extent_0 = std::extent<ArrayType>::value;
-  value_type internal_data_[extent_0];
-
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator()(int i) const noexcept
-  {
-    return internal_data_[i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator[](int i) const noexcept
-  {
-    return internal_data_[i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION value_type& operator()(int i) noexcept
-  {
-    return internal_data_[i];
-  }
-};
-
-template <typename ArrayType>
-struct alignas(alignment) LocalArray<
-  ArrayType,
-  typename std::enable_if<std::rank<ArrayType>::value == 2>::type>
-{
-  static constexpr int Rank = 2;
-  static constexpr int extent_0 = std::extent<ArrayType, 0>::value;
-  static constexpr int extent_1 = std::extent<ArrayType, 1>::value;
-  using value_type = typename std::remove_all_extents<ArrayType>::type;
-  value_type internal_data_[extent_0][extent_1];
-
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator()(int j, int i) const noexcept
-  {
-    return internal_data_[j][i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION value_type& operator()(int j, int i) noexcept
-  {
-    return internal_data_[j][i];
-  }
-};
-
-template <typename ArrayType>
-struct alignas(alignment) LocalArray<
-  ArrayType,
-  typename std::enable_if<std::rank<ArrayType>::value == 3>::type>
-{
-  static constexpr int Rank = 3;
-  static constexpr int extent_0 = std::extent<ArrayType, 0>::value;
-  static constexpr int extent_1 = std::extent<ArrayType, 1>::value;
-  static constexpr int extent_2 = std::extent<ArrayType, 2>::value;
-  using value_type = typename std::remove_all_extents<ArrayType>::type;
-  value_type internal_data_[extent_0][extent_1][extent_2];
-
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator()(int k, int j, int i) const noexcept
-  {
-    return internal_data_[k][j][i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION value_type&
-  operator()(int k, int j, int i) noexcept
-  {
-    return internal_data_[k][j][i];
-  }
-};
-
-template <typename ArrayType>
-struct alignas(alignment) LocalArray<
-  ArrayType,
-  typename std::enable_if<std::rank<ArrayType>::value == 4>::type>
-{
-  static constexpr int Rank = 4;
-  static constexpr int extent_0 = std::extent<ArrayType, 0>::value;
-  static constexpr int extent_1 = std::extent<ArrayType, 1>::value;
-  static constexpr int extent_2 = std::extent<ArrayType, 2>::value;
-  static constexpr int extent_3 = std::extent<ArrayType, 3>::value;
-  using value_type = typename std::remove_all_extents<ArrayType>::type;
-  value_type internal_data_[extent_0][extent_1][extent_2][extent_3];
-
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator()(int l, int k, int j, int i) const noexcept
-  {
-    return internal_data_[l][k][j][i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION value_type&
-  operator()(int l, int k, int j, int i) noexcept
-  {
-    return internal_data_[l][k][j][i];
-  }
-};
-
-template <typename ArrayType>
-struct alignas(alignment) LocalArray<
-  ArrayType,
-  typename std::enable_if<std::rank<ArrayType>::value == 5>::type>
-{
-  static constexpr int Rank = 5;
-  static constexpr int extent_0 = std::extent<ArrayType, 0>::value;
-  static constexpr int extent_1 = std::extent<ArrayType, 1>::value;
-  static constexpr int extent_2 = std::extent<ArrayType, 2>::value;
-  static constexpr int extent_3 = std::extent<ArrayType, 3>::value;
-  static constexpr int extent_4 = std::extent<ArrayType, 4>::value;
-  using value_type = typename std::remove_all_extents<ArrayType>::type;
-  value_type internal_data_[extent_0][extent_1][extent_2][extent_3][extent_4];
-
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator()(int m, int l, int k, int j, int i) const noexcept
-  {
-    return internal_data_[m][l][k][j][i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION value_type&
-  operator()(int m, int l, int k, int j, int i) noexcept
-  {
-    return internal_data_[m][l][k][j][i];
-  }
-};
-
-template <typename ArrayType>
-struct alignas(alignment) LocalArray<
-  ArrayType,
-  typename std::enable_if<std::rank<ArrayType>::value == 6>::type>
-{
-  static constexpr int Rank = 6;
-  static constexpr int extent_0 = std::extent<ArrayType, 0>::value;
-  static constexpr int extent_1 = std::extent<ArrayType, 1>::value;
-  static constexpr int extent_2 = std::extent<ArrayType, 2>::value;
-  static constexpr int extent_3 = std::extent<ArrayType, 3>::value;
-  static constexpr int extent_4 = std::extent<ArrayType, 4>::value;
-  static constexpr int extent_5 = std::extent<ArrayType, 5>::value;
-
-  using value_type = typename std::remove_all_extents<ArrayType>::type;
-  value_type internal_data_[extent_0][extent_1][extent_2][extent_3][extent_4]
-                           [extent_5];
-
-  KOKKOS_FORCEINLINE_FUNCTION constexpr value_type
-  operator()(int n, int m, int l, int k, int j, int i) const noexcept
-  {
-    return internal_data_[n][m][l][k][j][i];
-  }
-  KOKKOS_FORCEINLINE_FUNCTION value_type&
-  operator()(int n, int m, int l, int k, int j, int i) noexcept
-  {
-    return internal_data_[n][m][l][k][j][i];
-  }
-};
-
-} // namespace matrix_free
-} // namespace nalu
-} // namespace sierra
-
-#endif
diff --git a/include/matrix_free/MaxCourantReynolds.h b/include/matrix_free/MaxCourantReynolds.h
index 10a6735167..1dd16c923a 100644
--- a/include/matrix_free/MaxCourantReynolds.h
+++ b/include/matrix_free/MaxCourantReynolds.h
@@ -12,7 +12,7 @@
 
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
diff --git a/include/matrix_free/MomentumInterior.h b/include/matrix_free/MomentumInterior.h
index 2b69460a44..431fea93c6 100644
--- a/include/matrix_free/MomentumInterior.h
+++ b/include/matrix_free/MomentumInterior.h
@@ -11,7 +11,7 @@
 #define MOMENTUM_INTERIOR_H
 
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 
 #include "Kokkos_Array.hpp"
@@ -32,7 +32,7 @@ namespace impl {
 template <int p>
 struct momentum_residual_t
 {
-  using narray = LocalArray<ftype[p + 1][p + 1][p + 1]>;
+  using narray = ArrayND<ftype[p + 1][p + 1][p + 1]>;
 
   static void invoke(
     Kokkos::Array<double, 3> gammas,
@@ -57,7 +57,7 @@ namespace impl {
 template <int p>
 struct momentum_linearized_residual_t
 {
-  using narray = LocalArray<ftype[p + 1][p + 1][p + 1]>;
+  using narray = ArrayND<ftype[p + 1][p + 1][p + 1]>;
 
   static void invoke(
     double proj_time_scale,
diff --git a/include/matrix_free/NodeOrderMap.h b/include/matrix_free/NodeOrderMap.h
index c69636c436..5051ae3bc3 100644
--- a/include/matrix_free/NodeOrderMap.h
+++ b/include/matrix_free/NodeOrderMap.h
@@ -10,7 +10,7 @@
 #ifndef NODEORDERMAP_H
 #define NODEORDERMAP_H
 
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
@@ -24,14 +24,14 @@ struct StkNodeOrderMapping
 template <>
 struct StkNodeOrderMapping<1>
 {
-  using node_map_type = LocalArray<int[2][2][2]>;
+  using node_map_type = ArrayND<int[2][2][2]>;
   static constexpr node_map_type map = {{{{0, 1}, {3, 2}}, {{4, 5}, {7, 6}}}};
 };
 
 template <>
 struct StkNodeOrderMapping<2>
 {
-  using node_map_type = LocalArray<int[3][3][3]>;
+  using node_map_type = ArrayND<int[3][3][3]>;
   static constexpr node_map_type map = {
     {{{0, 8, 1}, {11, 21, 9}, {3, 10, 2}},
      {{12, 25, 13}, {23, 20, 24}, {15, 26, 14}},
@@ -41,7 +41,7 @@ struct StkNodeOrderMapping<2>
 template <>
 struct StkNodeOrderMapping<3>
 {
-  using node_map_type = LocalArray<int[4][4][4]>;
+  using node_map_type = ArrayND<int[4][4][4]>;
   static constexpr node_map_type map = {
     {{{0, 8, 9, 1}, {15, 32, 34, 10}, {14, 33, 35, 11}, {3, 13, 12, 2}},
      {{16, 48, 49, 18}, {40, 56, 57, 44}, {42, 58, 59, 45}, {22, 53, 52, 20}},
@@ -52,7 +52,7 @@ struct StkNodeOrderMapping<3>
 template <>
 struct StkNodeOrderMapping<4>
 {
-  using node_map_type = LocalArray<int[5][5][5]>;
+  using node_map_type = ArrayND<int[5][5][5]>;
   static constexpr node_map_type map = {
     {{{0, 8, 9, 10, 1},
       {19, 44, 47, 50, 11},
@@ -89,21 +89,21 @@ struct StkFaceNodeMapping
 template <>
 struct StkFaceNodeMapping<1>
 {
-  using node_map_type = LocalArray<int[2][2]>;
+  using node_map_type = ArrayND<int[2][2]>;
   static constexpr node_map_type map = {{{0, 1}, {3, 2}}};
 };
 
 template <>
 struct StkFaceNodeMapping<2>
 {
-  using node_map_type = LocalArray<int[3][3]>;
+  using node_map_type = ArrayND<int[3][3]>;
   static constexpr node_map_type map = {{{0, 4, 1}, {7, 8, 5}, {3, 6, 2}}};
 };
 
 template <>
 struct StkFaceNodeMapping<3>
 {
-  using node_map_type = LocalArray<int[4][4]>;
+  using node_map_type = ArrayND<int[4][4]>;
   static constexpr node_map_type map = {
     {{0, 4, 5, 1}, {11, 12, 13, 6}, {10, 14, 15, 7}, {3, 9, 8, 2}}};
 };
@@ -111,7 +111,7 @@ struct StkFaceNodeMapping<3>
 template <>
 struct StkFaceNodeMapping<4>
 {
-  using node_map_type = LocalArray<int[5][5]>;
+  using node_map_type = ArrayND<int[5][5]>;
   static constexpr node_map_type map = {
     {{0, 4, 5, 6, 1},
      {15, 16, 17, 18, 7},
diff --git a/include/matrix_free/ScalarFluxBC.h b/include/matrix_free/ScalarFluxBC.h
index 21ff3b81cd..b974361624 100644
--- a/include/matrix_free/ScalarFluxBC.h
+++ b/include/matrix_free/ScalarFluxBC.h
@@ -14,7 +14,7 @@
 
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
diff --git a/include/matrix_free/ShuffledAccess.h b/include/matrix_free/ShuffledAccess.h
index 87de7e5be1..a21fd9d2a8 100644
--- a/include/matrix_free/ShuffledAccess.h
+++ b/include/matrix_free/ShuffledAccess.h
@@ -11,7 +11,7 @@
 #define SHUFFLED_ACCESS_H
 
 #include "matrix_free/KokkosFramework.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include "Kokkos_Macros.hpp"
 
diff --git a/include/matrix_free/StrongDirichletBC.h b/include/matrix_free/StrongDirichletBC.h
index 38b6f4e952..b0d54eee55 100644
--- a/include/matrix_free/StrongDirichletBC.h
+++ b/include/matrix_free/StrongDirichletBC.h
@@ -13,7 +13,7 @@
 #include <Tpetra_MultiVector.hpp>
 
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
diff --git a/include/matrix_free/TensorOperations.h b/include/matrix_free/TensorOperations.h
index bdf4ca6077..e457c6276d 100644
--- a/include/matrix_free/TensorOperations.h
+++ b/include/matrix_free/TensorOperations.h
@@ -10,7 +10,7 @@
 #ifndef TENSOR_OPERATIONS_H
 #define TENSOR_OPERATIONS_H
 
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include "Kokkos_Macros.hpp"
 
@@ -24,7 +24,7 @@ namespace matrix_free {
 
 template <typename Scalar>
 KOKKOS_FORCEINLINE_FUNCTION Scalar
-determinant(const LocalArray<Scalar[3][3]>& jac)
+determinant(const ArrayND<Scalar[3][3]>& jac)
 {
   return jac(XH, XH) * (jac(YH, YH) * jac(ZH, ZH) - jac(YH, ZH) * jac(ZH, YH)) -
          jac(XH, YH) * (jac(YH, XH) * jac(ZH, ZH) - jac(YH, ZH) * jac(ZH, XH)) +
@@ -33,8 +33,7 @@ determinant(const LocalArray<Scalar[3][3]>& jac)
 
 template <typename Scalar>
 KOKKOS_FORCEINLINE_FUNCTION void
-adjugate_matrix(
-  const LocalArray<Scalar[3][3]>& mat, LocalArray<Scalar[3][3]>& adj)
+adjugate_matrix(const ArrayND<Scalar[3][3]>& mat, ArrayND<Scalar[3][3]>& adj)
 {
   adj(XH, XH) = mat(YH, YH) * mat(ZH, ZH) - mat(ZH, YH) * mat(YH, ZH);
   adj(XH, YH) = mat(YH, ZH) * mat(ZH, XH) - mat(ZH, ZH) * mat(YH, XH);
@@ -48,10 +47,10 @@ adjugate_matrix(
 }
 
 template <typename Scalar>
-KOKKOS_FORCEINLINE_FUNCTION LocalArray<Scalar[3][3]>
-invert_matrix(const LocalArray<Scalar[3][3]>& mat)
+KOKKOS_FORCEINLINE_FUNCTION ArrayND<Scalar[3][3]>
+invert_matrix(const ArrayND<Scalar[3][3]>& mat)
 {
-  LocalArray<Scalar[3][3]> inv;
+  ArrayND<Scalar[3][3]> inv;
   auto inv_detj = 1. / determinant(mat);
   inv(XH, XH) =
     inv_detj * (mat(YH, YH) * mat(ZH, ZH) - mat(ZH, YH) * mat(YH, ZH));
@@ -76,10 +75,10 @@ invert_matrix(const LocalArray<Scalar[3][3]>& mat)
 }
 
 template <typename Scalar>
-KOKKOS_FORCEINLINE_FUNCTION LocalArray<Scalar[3][3]>
-adjugate_matrix(const LocalArray<Scalar[3][3]>& mat)
+KOKKOS_FORCEINLINE_FUNCTION ArrayND<Scalar[3][3]>
+adjugate_matrix(const ArrayND<Scalar[3][3]>& mat)
 {
-  return LocalArray<Scalar[3][3]>{
+  return ArrayND<Scalar[3][3]>{
     {{mat(YH, YH) * mat(ZH, ZH) - mat(ZH, YH) * mat(YH, ZH),
       mat(YH, ZH) * mat(ZH, XH) - mat(ZH, ZH) * mat(YH, XH),
       mat(YH, XH) * mat(ZH, YH) - mat(ZH, XH) * mat(YH, YH)},
@@ -94,7 +93,7 @@ adjugate_matrix(const LocalArray<Scalar[3][3]>& mat)
 template <typename InpScalar, typename OutScalar>
 KOKKOS_FORCEINLINE_FUNCTION void
 transform(
-  const LocalArray<InpScalar[3][3]>& A,
+  const ArrayND<InpScalar[3][3]>& A,
   const Kokkos::Array<OutScalar, 3>& x,
   Kokkos::Array<OutScalar, 3>& y)
 {
@@ -106,9 +105,9 @@ transform(
 template <typename InpScalar, typename OutScalar>
 KOKKOS_FORCEINLINE_FUNCTION void
 transform(
-  const LocalArray<InpScalar[3][3]>& A,
-  const LocalArray<OutScalar[3][3]>& x,
-  LocalArray<OutScalar[3][3]>& y)
+  const ArrayND<InpScalar[3][3]>& A,
+  const ArrayND<OutScalar[3][3]>& x,
+  ArrayND<OutScalar[3][3]>& y)
 {
   for (int d = 0; d < 3; ++d) {
     y(d, XH) =
@@ -123,9 +122,9 @@ transform(
 template <typename InpScalar, typename OutScalar>
 KOKKOS_FORCEINLINE_FUNCTION void
 inv_transform_t(
-  const LocalArray<InpScalar[3][3]>& A,
-  const LocalArray<OutScalar[3]>& b,
-  LocalArray<OutScalar[3]>& x)
+  const ArrayND<InpScalar[3][3]>& A,
+  const ArrayND<OutScalar[3]>& b,
+  ArrayND<OutScalar[3]>& x)
 {
   auto inv_det = 1. / determinant(A);
   x(0) = ((A(YH, YH) * A(ZH, ZH) - A(YH, ZH) * A(ZH, YH)) * b(0) +
@@ -147,9 +146,9 @@ inv_transform_t(
 template <typename InpScalar, typename OutScalar>
 KOKKOS_FORCEINLINE_FUNCTION void
 inv_transform_t(
-  const LocalArray<InpScalar[3][3]>& A,
-  const LocalArray<OutScalar[3][3]>& b,
-  LocalArray<OutScalar[3][3]>& x)
+  const ArrayND<InpScalar[3][3]>& A,
+  const ArrayND<OutScalar[3][3]>& b,
+  ArrayND<OutScalar[3][3]>& x)
 {
   auto inv_det = 1. / determinant(A);
   for (int d = 0; d < 3; ++d) {
@@ -171,10 +170,10 @@ inv_transform_t(
 }
 
 template <typename Scalar>
-KOKKOS_FORCEINLINE_FUNCTION LocalArray<Scalar[3][3]>
-square(const LocalArray<Scalar[3][3]>& a)
+KOKKOS_FORCEINLINE_FUNCTION ArrayND<Scalar[3][3]>
+square(const ArrayND<Scalar[3][3]>& a)
 {
-  LocalArray<ftype[3][3]> b;
+  ArrayND<ftype[3][3]> b;
   for (int dj = 0; dj < 3; ++dj) {
     for (int di = 0; di < 3; ++di) {
       b(dj, di) =
diff --git a/src/matrix_free/Coefficients.C b/src/matrix_free/Coefficients.C
index 969d947ccd..8ab09de197 100644
--- a/src/matrix_free/Coefficients.C
+++ b/src/matrix_free/Coefficients.C
@@ -8,7 +8,7 @@
 //
 
 #include "matrix_free/Coefficients.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
@@ -42,10 +42,10 @@ constexpr Coeffs<4>::scs_matrix_type Coeffs<4>::Dt;
 constexpr Coeffs<4>::linear_nodal_matrix_type Coeffs<4>::Nlin;
 constexpr Coeffs<4>::linear_scs_matrix_type Coeffs<4>::Ntlin;
 
-constexpr LocalArray<double[2]> Coeffs<1>::Wl;
-constexpr LocalArray<double[3]> Coeffs<2>::Wl;
-constexpr LocalArray<double[4]> Coeffs<3>::Wl;
-constexpr LocalArray<double[5]> Coeffs<4>::Wl;
+constexpr ArrayND<double[2]> Coeffs<1>::Wl;
+constexpr ArrayND<double[3]> Coeffs<2>::Wl;
+constexpr ArrayND<double[4]> Coeffs<3>::Wl;
+constexpr ArrayND<double[5]> Coeffs<4>::Wl;
 
 } // namespace matrix_free
 } // namespace nalu
diff --git a/src/matrix_free/ConductionDiagonal.C b/src/matrix_free/ConductionDiagonal.C
index 3ecb7ca15c..931e9fe761 100644
--- a/src/matrix_free/ConductionDiagonal.C
+++ b/src/matrix_free/ConductionDiagonal.C
@@ -14,7 +14,7 @@
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/ShuffledAccess.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include <KokkosInterface.h>
 #include <Kokkos_ScatterView.hpp>
@@ -82,7 +82,7 @@ conduction_diagonal_t<p>::invoke(
       constexpr auto vandermonde = Coeffs<p>::W;
       constexpr auto Wl = Coeffs<p>::Wl;
 
-      LocalArray<ftype[p + 1][p + 1][p + 1]> lhs;
+      ArrayND<ftype[p + 1][p + 1][p + 1]> lhs;
       for (int k = 0; k < p + 1; ++k) {
         const auto gammaWk = gamma * Wl(k);
         for (int j = 0; j < p + 1; ++j) {
diff --git a/src/matrix_free/ConductionInterior.C b/src/matrix_free/ConductionInterior.C
index 41f55bd003..488ea1502f 100644
--- a/src/matrix_free/ConductionInterior.C
+++ b/src/matrix_free/ConductionInterior.C
@@ -14,7 +14,7 @@
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include <KokkosInterface.h>
 #include <Kokkos_ScatterView.hpp>
@@ -94,7 +94,7 @@ conduction_linearized_residual_t<p>::invoke(
   Kokkos::parallel_for(
     "conduction_linop", offsets.extent_int(0), KOKKOS_LAMBDA(int index) {
       narray delta;
-      LocalArray<int[p + 1][p + 1][p + 1][simd_len]> idx;
+      ArrayND<int[p + 1][p + 1][p + 1][simd_len]> idx;
       const auto valid_length = valid_offset<p>(index, offsets);
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
diff --git a/src/matrix_free/ContinuityInterior.C b/src/matrix_free/ContinuityInterior.C
index 6bf35a1f7e..880472a85a 100644
--- a/src/matrix_free/ContinuityInterior.C
+++ b/src/matrix_free/ContinuityInterior.C
@@ -12,7 +12,7 @@
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/ElementFluxIntegral.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/ValidSimdLength.h"
 
@@ -41,7 +41,7 @@ continuity_residual_t<p>::invoke(
   auto yout_scatter = Kokkos::Experimental::create_scatter_view(yout);
   Kokkos::parallel_for(
     DeviceRangePolicy(0, offsets.extent_int(0)), KOKKOS_LAMBDA(int index) {
-      LocalArray<ftype[p + 1][p + 1][p + 1]> elem_rhs;
+      ArrayND<ftype[p + 1][p + 1][p + 1]> elem_rhs;
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
           for (int i = 0; i < p + 1; ++i) {
@@ -85,7 +85,7 @@ continuity_linearized_residual_t<p>::invoke(
   Kokkos::parallel_for(
     DeviceRangePolicy(0, offsets.extent_int(0)), KOKKOS_LAMBDA(int index) {
       narray delta;
-      LocalArray<int[p + 1][p + 1][p + 1][simd_len]> idx;
+      ArrayND<int[p + 1][p + 1][p + 1][simd_len]> idx;
       const auto valid_length = valid_offset<p>(index, offsets);
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
diff --git a/src/matrix_free/FilterDiagonal.C b/src/matrix_free/FilterDiagonal.C
index 722dcefc0d..a215171e49 100644
--- a/src/matrix_free/FilterDiagonal.C
+++ b/src/matrix_free/FilterDiagonal.C
@@ -13,7 +13,7 @@
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/LinSysInfo.h"
 #include <KokkosInterface.h>
 
@@ -40,7 +40,7 @@ filter_diagonal_t<p>::invoke(
   auto yout_scatter = Kokkos::Experimental::create_scatter_view(yout);
   Kokkos::parallel_for(
     DeviceRangePolicy(0, offsets.extent_int(0)), KOKKOS_LAMBDA(int index) {
-      LocalArray<ftype[p + 1][p + 1][p + 1]> lhs;
+      ArrayND<ftype[p + 1][p + 1][p + 1]> lhs;
       if (lumped) {
         for (int k = 0; k < p + 1; ++k) {
           static constexpr auto Wl = Coeffs<p>::Wl;
diff --git a/src/matrix_free/GreenGaussBoundaryClosure.C b/src/matrix_free/GreenGaussBoundaryClosure.C
index bf89307b1a..d98887c169 100644
--- a/src/matrix_free/GreenGaussBoundaryClosure.C
+++ b/src/matrix_free/GreenGaussBoundaryClosure.C
@@ -12,7 +12,7 @@
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/KokkosFramework.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/ValidSimdLength.h"
 
@@ -42,7 +42,7 @@ vector_component_flux(
   FaceRankOutput& out,
   int component)
 {
-  LocalArray<ftype[p + 1][p + 1]> scratch;
+  ArrayND<ftype[p + 1][p + 1]> scratch;
 
   for (int j = 0; j < p + 1; ++j) {
     for (int i = 0; i < p + 1; ++i) {
@@ -89,7 +89,7 @@ gradient_boundary_closure_t<p>::invoke(
   Kokkos::parallel_for(
     DeviceRangePolicy(0, offsets.extent_int(0)), KOKKOS_LAMBDA(int index) {
       for (int d = 0; d < 3; ++d) {
-        LocalArray<ftype[p + 1][p + 1]> rhs;
+        ArrayND<ftype[p + 1][p + 1]> rhs;
         vector_component_flux<p>(index, q, areav, rhs, d);
 
         auto accessor = yout_scatter.access();
diff --git a/src/matrix_free/GreenGaussGradientInterior.C b/src/matrix_free/GreenGaussGradientInterior.C
index 5958d8e84c..0a5df9db66 100644
--- a/src/matrix_free/GreenGaussGradientInterior.C
+++ b/src/matrix_free/GreenGaussGradientInterior.C
@@ -14,7 +14,7 @@
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include <KokkosInterface.h>
 #include "Kokkos_ScatterView.hpp"
@@ -43,7 +43,7 @@ gradient_residual_t<p>::invoke(
   auto range = range_type({0, 0}, {offsets.extent_int(0), 3});
   Kokkos::parallel_for(
     range, KOKKOS_LAMBDA(int index, int d) {
-      LocalArray<ftype[p + 1][p + 1][p + 1]> rhs;
+      ArrayND<ftype[p + 1][p + 1][p + 1]> rhs;
       if (lumped) {
         for (int k = 0; k < p + 1; ++k) {
           for (int j = 0; j < p + 1; ++j) {
@@ -55,7 +55,7 @@ gradient_residual_t<p>::invoke(
           }
         }
       } else {
-        LocalArray<ftype[p + 1][p + 1][p + 1]> scratch;
+        ArrayND<ftype[p + 1][p + 1][p + 1]> scratch;
         for (int k = 0; k < p + 1; ++k) {
           for (int j = 0; j < p + 1; ++j) {
             for (int i = 0; i < p + 1; ++i) {
@@ -116,7 +116,7 @@ filter_linearized_residual_t<p>::invoke(
   Kokkos::parallel_for(
     offsets.extent_int(0), KOKKOS_LAMBDA(int index) {
       const auto length = valid_offset<p>(index, offsets);
-      LocalArray<int[p + 1][p + 1][p + 1][simd_len]> idx;
+      ArrayND<int[p + 1][p + 1][p + 1][simd_len]> idx;
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
           for (int i = 0; i < p + 1; ++i) {
@@ -128,7 +128,7 @@ filter_linearized_residual_t<p>::invoke(
       }
 
       for (int d = 0; d < 3; ++d) {
-        LocalArray<ftype[p + 1][p + 1][p + 1]> delta;
+        ArrayND<ftype[p + 1][p + 1][p + 1]> delta;
         for (int k = 0; k < p + 1; ++k) {
           for (int j = 0; j < p + 1; ++j) {
             for (int i = 0; i < p + 1; ++i) {
@@ -139,7 +139,7 @@ filter_linearized_residual_t<p>::invoke(
           }
         }
 
-        LocalArray<ftype[p + 1][p + 1][p + 1]> rhs;
+        ArrayND<ftype[p + 1][p + 1][p + 1]> rhs;
         if (lumped) {
           for (int k = 0; k < p + 1; ++k) {
             for (int j = 0; j < p + 1; ++j) {
diff --git a/src/matrix_free/LinearAdvectionMetric.C b/src/matrix_free/LinearAdvectionMetric.C
index a998e66d31..9928aca0c6 100644
--- a/src/matrix_free/LinearAdvectionMetric.C
+++ b/src/matrix_free/LinearAdvectionMetric.C
@@ -12,7 +12,7 @@
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/KokkosFramework.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/ShuffledAccess.h"
 
@@ -60,7 +60,7 @@ corrected_momentum_flux_coefficient(
   }
 
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[p + 1][p + 1]> scratch;
+    ArrayND<ftype[p + 1][p + 1]> scratch;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         ftype acc = 0;
@@ -113,7 +113,7 @@ linear_advection_metric_t<p>::invoke(
   enum { XH = 0, YH = 1, ZH = 2 };
   Kokkos::parallel_for(
     DeviceRangePolicy(0, areas.extent_int(0)), KOKKOS_LAMBDA(int index) {
-      LocalArray<ftype[p + 1][p + 1][p + 1][3]> rhou_corr;
+      ArrayND<ftype[p + 1][p + 1][p + 1][3]> rhou_corr;
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
           for (int i = 0; i < p + 1; ++i) {
diff --git a/src/matrix_free/LinearDiffusionMetric.C b/src/matrix_free/LinearDiffusionMetric.C
index 3e0139ca55..1ce47ba7aa 100644
--- a/src/matrix_free/LinearDiffusionMetric.C
+++ b/src/matrix_free/LinearDiffusionMetric.C
@@ -14,7 +14,7 @@
 #include "matrix_free/HexVertexCoordinates.h"
 #include "matrix_free/GeometricFunctions.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 
 #include <KokkosInterface.h>
@@ -38,7 +38,7 @@ diffusion_metric_t<p>::invoke(
     KOKKOS_LAMBDA(int index) {
       static constexpr auto ntilde = Coeffs<p>::Nt;
 
-      LocalArray<ftype[3][p + 1][p + 1][p + 1]> interp;
+      ArrayND<ftype[3][p + 1][p + 1][p + 1]> interp;
       {
         const auto alpha_elem = Kokkos::subview(
           alpha, index, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
diff --git a/src/matrix_free/LinearExposedAreas.C b/src/matrix_free/LinearExposedAreas.C
index 3180c81384..0c9c1a8aec 100644
--- a/src/matrix_free/LinearExposedAreas.C
+++ b/src/matrix_free/LinearExposedAreas.C
@@ -16,7 +16,7 @@
 #include "matrix_free/Coefficients.h"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
@@ -42,9 +42,9 @@ jacobian_component(
 }
 
 template <typename CoeffArray>
-KOKKOS_FUNCTION LocalArray<ftype[3]>
+KOKKOS_FUNCTION ArrayND<ftype[3]>
 face_area(
-  const LocalArray<ftype[3][4]>& base_box, const CoeffArray& nlin, int j, int i)
+  const ArrayND<ftype[3][4]>& base_box, const CoeffArray& nlin, int j, int i)
 {
   enum { XH = 0, YH = 1, ZH = 2 };
   enum { DS1 = 0, DS2 = 1 };
@@ -55,7 +55,7 @@ face_area(
   const auto dy_ds2 = jacobian_component<YH, DS2>(base_box, nlin, j, i);
   const auto dz_ds1 = jacobian_component<ZH, DS1>(base_box, nlin, j, i);
   const auto dz_ds2 = jacobian_component<ZH, DS2>(base_box, nlin, j, i);
-  return LocalArray<ftype[3]>{
+  return ArrayND<ftype[3]>{
     {dy_ds1 * dz_ds2 - dz_ds1 * dy_ds2, dz_ds1 * dx_ds2 - dx_ds1 * dz_ds2,
      dx_ds1 * dy_ds2 - dy_ds1 * dx_ds2}};
 }
diff --git a/src/matrix_free/LinearVolume.C b/src/matrix_free/LinearVolume.C
index 6dd35cc2f3..0a9f1e4228 100644
--- a/src/matrix_free/LinearVolume.C
+++ b/src/matrix_free/LinearVolume.C
@@ -13,7 +13,7 @@
 #include "matrix_free/GeometricFunctions.h"
 #include "matrix_free/KokkosFramework.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/TensorOperations.h"
 
diff --git a/src/matrix_free/LocalDualNodalVolume.C b/src/matrix_free/LocalDualNodalVolume.C
index 44cfa53248..cda16ade83 100644
--- a/src/matrix_free/LocalDualNodalVolume.C
+++ b/src/matrix_free/LocalDualNodalVolume.C
@@ -34,7 +34,7 @@ compute_volumes(const BoxArray& box, OutArray& out)
 
   for (int k = 0; k < p + 1; ++k) {
     for (int j = 0; j < p + 1; ++j) {
-      LocalArray<ftype[p + 1]> scratch;
+      ArrayND<ftype[p + 1]> scratch;
       for (int i = 0; i < p + 1; ++i) {
         scratch(i) =
           determinant<ftype>(geom::linear_hex_jacobian<p>(box, k, j, i));
@@ -50,7 +50,7 @@ compute_volumes(const BoxArray& box, OutArray& out)
   }
 
   for (int i = 0; i < p + 1; ++i) {
-    LocalArray<ftype[p + 1][p + 1]> scratch;
+    ArrayND<ftype[p + 1][p + 1]> scratch;
     for (int k = 0; k < p + 1; ++k) {
       for (int j = 0; j < p + 1; ++j) {
         ftype acc(0);
@@ -91,7 +91,7 @@ local_dual_nodal_volume_t<p>::invoke(
       const auto box = hex_vertex_coordinates<p>(index, xc);
       const auto valid_length = valid_offset<p>(index, conn);
 
-      LocalArray<ftype[p + 1][p + 1][p + 1]> vols;
+      ArrayND<ftype[p + 1][p + 1][p + 1]> vols;
       compute_volumes<p>(box, vols);
 
       for (int k = 0; k < p + 1; ++k) {
diff --git a/src/matrix_free/MaxCourantReynolds.C b/src/matrix_free/MaxCourantReynolds.C
index 5615f28d3f..85453f4c88 100644
--- a/src/matrix_free/MaxCourantReynolds.C
+++ b/src/matrix_free/MaxCourantReynolds.C
@@ -14,7 +14,7 @@
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/ElementVolumeIntegral.h"
 #include "matrix_free/ElementSCSInterpolate.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include <KokkosInterface.h>
 #include <Kokkos_NumericTraits.hpp>
 #include <stk_simd/Simd.hpp>
diff --git a/src/matrix_free/MomentumDiagonal.C b/src/matrix_free/MomentumDiagonal.C
index 60622c41ac..26cc6691df 100644
--- a/src/matrix_free/MomentumDiagonal.C
+++ b/src/matrix_free/MomentumDiagonal.C
@@ -14,7 +14,7 @@
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/ShuffledAccess.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include <KokkosInterface.h>
 #include <Kokkos_ScatterView.hpp>
 
@@ -81,7 +81,7 @@ advdiff_diagonal_t<p>::invoke(
   auto yout_scatter = Kokkos::Experimental::create_scatter_view(yout);
   Kokkos::parallel_for(
     DeviceRangePolicy(0, offsets.extent_int(0)), KOKKOS_LAMBDA(int index) {
-      LocalArray<ftype[p + 1][p + 1][p + 1]> lhs;
+      ArrayND<ftype[p + 1][p + 1][p + 1]> lhs;
       // generally works better to lump the mass term to the diagonal here
       for (int k = 0; k < p + 1; ++k) {
         static constexpr auto Wl = Coeffs<p>::Wl;
diff --git a/src/matrix_free/MomentumInterior.C b/src/matrix_free/MomentumInterior.C
index 27cc034b17..2e8f10a13a 100644
--- a/src/matrix_free/MomentumInterior.C
+++ b/src/matrix_free/MomentumInterior.C
@@ -20,7 +20,7 @@
 #include "matrix_free/TensorOperations.h"
 #include "matrix_free/ValidSimdLength.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/ElementSCSInterpolate.h"
 #include <KokkosInterface.h>
 
@@ -41,13 +41,13 @@ template <
   typename AdvArrayType,
   typename ViscosityArrayType,
   typename UArrayType>
-KOKKOS_FORCEINLINE_FUNCTION LocalArray<ftype[3]>
+KOKKOS_FORCEINLINE_FUNCTION ArrayND<ftype[3]>
 momentum_flux(
   const BoxArrayType& box,
   const AdvArrayType& adv,
   const ViscosityArrayType& visc,
   const UArrayType& u,
-  const LocalArray<ftype[p + 1][p + 1][3]>& uhat,
+  const ArrayND<ftype[p + 1][p + 1][3]>& uhat,
   int l,
   int s,
   int r)
@@ -62,7 +62,7 @@ momentum_flux(
   const auto one_third_divu =
     (1. / 3.) * (gu(XH, XH) + gu(YH, YH) + gu(ZH, ZH));
 
-  LocalArray<ftype[3]> flux;
+  ArrayND<ftype[3]> flux;
   flux(0) = 2 * visc_ip *
               ((gu(XH, XH) - one_third_divu) * areav(XH) +
                0.5 * (gu(XH, YH) + gu(YH, XH)) * areav(YH) +
@@ -101,9 +101,9 @@ momentum_flux_difference(
   RHSArrayType& rhs)
 {
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[p + 1][p + 1][3]> flux;
+    ArrayND<ftype[p + 1][p + 1][3]> flux;
     {
-      LocalArray<ftype[p + 1][p + 1][3]> uhat;
+      ArrayND<ftype[p + 1][p + 1][3]> uhat;
       for (int s = 0; s < p + 1; ++s) {
         for (int r = 0; r < p + 1; ++r) {
           for (int d = 0; d < 3; ++d) {
@@ -124,7 +124,7 @@ momentum_flux_difference(
       }
     }
     for (int d = 0; d < 3; ++d) {
-      LocalArray<ftype[p + 1][p + 1]> fbar;
+      ArrayND<ftype[p + 1][p + 1]> fbar;
       for (int s = 0; s < p + 1; ++s) {
         for (int r = 0; r < p + 1; ++r) {
           ftype acc = 0;
@@ -169,7 +169,7 @@ momentum_mass(
   static constexpr auto vandermonde = Coeffs<p>::W;
   for (int k = 0; k < p + 1; ++k) {
     for (int j = 0; j < p + 1; ++j) {
-      LocalArray<ftype[p + 1][3]> scratch;
+      ArrayND<ftype[p + 1][3]> scratch;
       for (int i = 0; i < p + 1; ++i) {
         const auto vm1_val = vm1(index, k, j, i);
         const auto vp0_val = vp0(index, k, j, i);
@@ -198,7 +198,7 @@ momentum_mass(
 
   for (int d = 0; d < 3; ++d) {
     for (int i = 0; i < p + 1; ++i) {
-      LocalArray<ftype[p + 1][p + 1]> scratch;
+      ArrayND<ftype[p + 1][p + 1]> scratch;
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
           ftype acc(0);
@@ -246,7 +246,7 @@ momentum_residual_t<p>::invoke(
   auto yout_scatter = Kokkos::Experimental::create_scatter_view(yout);
   Kokkos::parallel_for(
     DeviceRangePolicy(0, offsets.extent_int(0)), KOKKOS_LAMBDA(int index) {
-      LocalArray<ftype[p + 1][p + 1][p + 1][3]> elem_rhs;
+      ArrayND<ftype[p + 1][p + 1][p + 1][3]> elem_rhs;
       if (p == 1) {
         static constexpr auto lumped = Coeffs<p>::Wl;
         for (int k = 0; k < p + 1; ++k) {
@@ -335,7 +335,7 @@ momentum_linearized_residual_t<p>::invoke(
   Kokkos::parallel_for(
     range, KOKKOS_LAMBDA(int index, int d) {
       const auto length = valid_offset<p>(index, offsets);
-      LocalArray<int[p + 1][p + 1][p + 1][simd_len]> idx;
+      ArrayND<int[p + 1][p + 1][p + 1][simd_len]> idx;
       narray delta;
       for (int k = 0; k < p + 1; ++k) {
         for (int j = 0; j < p + 1; ++j) {
diff --git a/src/matrix_free/NodeOrderMap.C b/src/matrix_free/NodeOrderMap.C
index 4c7a8bcc98..7f419e663a 100644
--- a/src/matrix_free/NodeOrderMap.C
+++ b/src/matrix_free/NodeOrderMap.C
@@ -8,7 +8,7 @@
 //
 
 #include "matrix_free/NodeOrderMap.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 namespace sierra {
 namespace nalu {
diff --git a/src/matrix_free/ScalarFluxBC.C b/src/matrix_free/ScalarFluxBC.C
index 9aaa3cae7a..6178ae077c 100644
--- a/src/matrix_free/ScalarFluxBC.C
+++ b/src/matrix_free/ScalarFluxBC.C
@@ -91,9 +91,9 @@ scalar_neumann_residual_t<p>::invoke(
   Kokkos::parallel_for(
     "flux_residual", DeviceRangePolicy(0, offsets.extent_int(0)),
     KOKKOS_LAMBDA(int index) {
-      LocalArray<ftype[p + 1][p + 1]> element_rhs;
+      ArrayND<ftype[p + 1][p + 1]> element_rhs;
       {
-        LocalArray<ftype[p + 1][p + 1]> scratch;
+        ArrayND<ftype[p + 1][p + 1]> scratch;
         scalar_flux<p>(index, dqdn, areav, scratch, element_rhs);
       }
 
diff --git a/src/matrix_free/SparsifiedEdgeLaplacian.C b/src/matrix_free/SparsifiedEdgeLaplacian.C
index 7bf4ae3ac5..2399a01694 100644
--- a/src/matrix_free/SparsifiedEdgeLaplacian.C
+++ b/src/matrix_free/SparsifiedEdgeLaplacian.C
@@ -135,7 +135,7 @@ assemble_sparsified_edge_laplacian_t<p>::invoke(
 {
   // map from the edge ordinals to an order
   // convenient for computing the edge laplacian
-  constexpr LocalArray<int[12][2][3]> edges = {{
+  constexpr ArrayND<int[12][2][3]> edges = {{
     {{0, 0, 0}, {0, 0, 1}}, // {0,1} .
     {{0, 1, 0}, {0, 1, 1}}, // {3,2} .
     {{1, 0, 0}, {1, 0, 1}}, // {4,5} .
@@ -163,7 +163,7 @@ assemble_sparsified_edge_laplacian_t<p>::invoke(
         for (int m = 0; m < p; ++m) {
           for (int l = 0; l < p; ++l) {
             const auto box = hex_vertex_coordinates(n, m, l, elem_coords);
-            LocalArray<ftype[12][2][2]> edge_lhs;
+            ArrayND<ftype[12][2][2]> edge_lhs;
             sparsified_laplacian_edge_lhs<p>(box, edge_lhs);
             for (int e = 0; e < 12; ++e) {
               const int ln = n + edges(e, 0, 0);
diff --git a/src/matrix_free/TransportCoefficients.C b/src/matrix_free/TransportCoefficients.C
index f4cf59f9bc..f80b6d429a 100644
--- a/src/matrix_free/TransportCoefficients.C
+++ b/src/matrix_free/TransportCoefficients.C
@@ -41,7 +41,7 @@ namespace matrix_free {
 namespace {
 
 KOKKOS_FUNCTION ftype
-wale_gradient_invariant(const LocalArray<ftype[3][3]>& dudx)
+wale_gradient_invariant(const ArrayND<ftype[3][3]>& dudx)
 {
   const auto dudx_sq = square(dudx);
   const auto one_third_trace =
@@ -67,7 +67,7 @@ wale_gradient_invariant(const LocalArray<ftype[3][3]>& dudx)
 }
 
 KOKKOS_FUNCTION ftype
-smag_gradient_invariant(const LocalArray<ftype[3][3]>& dudx)
+smag_gradient_invariant(const ArrayND<ftype[3][3]>& dudx)
 {
   ftype sij_sq = 0;
   for (int dj = 0; dj < 3; ++dj) {
diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt
index 2e4ccf8864..7179b5a864 100644
--- a/unit_tests/CMakeLists.txt
+++ b/unit_tests/CMakeLists.txt
@@ -1,5 +1,6 @@
 target_sources(${utest_ex_name} PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTest1ElemCoordCheck.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestArrayND.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestBasicKokkos.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestCopyAndInterleave.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestCreateOnDevice.C
diff --git a/unit_tests/matrix_free/StkSimdComparisons.h b/unit_tests/StkSimdComparisons.h
similarity index 84%
rename from unit_tests/matrix_free/StkSimdComparisons.h
rename to unit_tests/StkSimdComparisons.h
index 677018011d..a459baf657 100644
--- a/unit_tests/matrix_free/StkSimdComparisons.h
+++ b/unit_tests/StkSimdComparisons.h
@@ -11,45 +11,43 @@
 #define STK_SIMD_COMPARISONS_H
 
 #include "gtest/gtest.h"
-
-#include "matrix_free/KokkosViewTypes.h"
+#include "SimdInterface.h"
 
 // this allows calling the macros with a double instead of a DoubleType
 inline double
-get_simd_data_promote_double_to_doubletype(
-  sierra::nalu::matrix_free::ftype val, int ln)
+get_simd_data_promote_double_to_doubletype(sierra::nalu::DoubleType val, int ln)
 {
   return stk::simd::get_data(val, ln);
 }
 
 #define EXPECT_DOUBLETYPE_NEAR(val1, val2, abs_error)                          \
-  for (int ln = 0; ln < simd_len; ++ln)                                        \
+  for (int ln = 0; ln < simdLen; ++ln)                                         \
   EXPECT_PRED_FORMAT3(                                                         \
     ::testing::internal::DoubleNearPredFormat,                                 \
     get_simd_data_promote_double_to_doubletype(val1, ln),                      \
     get_simd_data_promote_double_to_doubletype(val2, ln), abs_error)
 
 #define ASSERT_DOUBLETYPE_NEAR(val1, val2, abs_error)                          \
-  for (int ln = 0; ln < simd_len; ++ln)                                        \
+  for (int ln = 0; ln < simdLen; ++ln)                                         \
   ASSERT_PRED_FORMAT3(                                                         \
     ::testing::internal::DoubleNearPredFormat,                                 \
     get_simd_data_promote_double_to_doubletype(val1, ln),                      \
     get_simd_data_promote_double_to_doubletype(val2, ln), abs_error)
 
 #define EXPECT_DOUBLETYPE_EQ(val1, val2)                                       \
-  for (int ln = 0; ln < simd_len; ++ln)                                        \
+  for (int ln = 0; ln < simdLen; ++ln)                                         \
   EXPECT_DOUBLE_EQ(                                                            \
     get_simd_data_promote_double_to_doubletype(val1, ln),                      \
     get_simd_data_promote_double_to_doubletype(val2, ln))
 
 #define EXPECT_DOUBLETYPE_GTEQ(val1, val2)                                     \
-  for (int ln = 0; ln < simd_len; ++ln)                                        \
+  for (int ln = 0; ln < simdLen; ++ln)                                         \
   EXPECT_PRED_FORMAT2(                                                         \
     ::testing::DoubleLE, get_simd_data_promote_double_to_doubletype(val2, ln), \
     get_simd_data_promote_double_to_doubletype(val1, ln))
 
 #define EXPECT_DOUBLETYPE_LTEQ(val1, val2)                                     \
-  for (int ln = 0; ln < simd_len; ++ln)                                        \
+  for (int ln = 0; ln < simdLen; ++ln)                                         \
   EXPECT_PRED_FORMAT2(                                                         \
     ::testing::DoubleLE, get_simd_data_promote_double_to_doubletype(val1, ln), \
     get_simd_data_promote_double_to_doubletype(val2, ln))
diff --git a/unit_tests/matrix_free/UnitTestLocalArray.C b/unit_tests/UnitTestArrayND.C
similarity index 51%
rename from unit_tests/matrix_free/UnitTestLocalArray.C
rename to unit_tests/UnitTestArrayND.C
index a82178847d..93dd401d67 100644
--- a/unit_tests/matrix_free/UnitTestLocalArray.C
+++ b/unit_tests/UnitTestArrayND.C
@@ -7,37 +7,35 @@
 // for more details.
 //
 
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/KokkosFramework.h"
 
 #include "StkSimdComparisons.h"
 #include <gtest/gtest.h>
 
-namespace sierra {
-namespace nalu {
-namespace matrix_free {
+namespace sierra::nalu {
 
-TEST(local_array, fill_double_1)
+TEST(array_nd, fill_double_1)
 {
-  LocalArray<double[4]> y = {{1, 1, 1, 1}};
+  ArrayND<double[4]> y = {{1, 1, 1, 1}};
   for (int i = 0; i < 4; ++i) {
     ASSERT_DOUBLE_EQ(y(i), 1.0);
     ASSERT_DOUBLE_EQ(y(i), y[i]);
   }
 }
 
-TEST(local_array, fill_ftypedouble_1)
+TEST(array_nd, fill_ftypedouble_1)
 {
-  LocalArray<ftype[4]> y = {{1, 1, 1, 1}};
+  ArrayND<DoubleType[4]> y = {{1, 1, 1, 1}};
   for (int i = 0; i < 4; ++i) {
     ASSERT_DOUBLETYPE_NEAR(y(i), 1.0, 1.0e-16);
     ASSERT_DOUBLETYPE_NEAR(y(i), y[i], 1.0e-16);
   }
 }
 
-TEST(local_array, fill_double_2)
+TEST(array_nd, fill_double_2)
 {
-  LocalArray<double[3][3]> y = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
+  ArrayND<double[3][3]> y = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
   for (int j = 0; j < 3; ++j) {
     for (int i = 0; i < 3; ++i) {
       ASSERT_DOUBLE_EQ(y(j, i), (i == j) ? 1.0 : 0.0);
@@ -45,9 +43,9 @@ TEST(local_array, fill_double_2)
   }
 }
 
-TEST(local_array, fill_double_3)
+TEST(array_nd, fill_double_3)
 {
-  LocalArray<double[3][3][3]> y;
+  ArrayND<double[3][3][3]> y;
   for (int k = 0; k < 3; ++k) {
     for (int j = 0; j < 3; ++j) {
       for (int i = 0; i < 3; ++i) {
@@ -58,9 +56,9 @@ TEST(local_array, fill_double_3)
   }
 }
 
-TEST(local_array, fill_double_4)
+TEST(array_nd, fill_double_4)
 {
-  LocalArray<double[3][3][4][3]> y;
+  ArrayND<double[3][3][4][3]> y;
   for (int l = 0; l < 3; ++l) {
     for (int k = 0; k < 3; ++k) {
       for (int j = 0; j < 4; ++j) {
@@ -73,9 +71,9 @@ TEST(local_array, fill_double_4)
   }
 }
 
-TEST(local_array, fill_double_5)
+TEST(array_nd, fill_double_5)
 {
-  LocalArray<double[7][3][3][4][3]> y;
+  ArrayND<double[7][3][3][4][3]> y;
   for (int m = 0; m < 7; ++m) {
     for (int l = 0; l < 3; ++l) {
       for (int k = 0; k < 3; ++k) {
@@ -90,6 +88,47 @@ TEST(local_array, fill_double_5)
   }
 }
 
-} // namespace matrix_free
-} // namespace nalu
-} // namespace sierra
+TEST(array_nd, static_rank)
+{
+  {
+    ArrayND<double[2]> s;
+    static_assert(decltype(s)::rank == 1);
+    static_assert(s.extent_int(0) == 2);
+  }
+
+  {
+    ArrayND<double[2][3]> s;
+    static_assert(decltype(s)::rank == 2);
+    static_assert(s.extent_int(0) == 2);
+    static_assert(s.extent_int(1) == 3);
+  }
+
+  {
+    ArrayND<double[2][3][4]> s;
+    static_assert(decltype(s)::rank == 3);
+    static_assert(s.extent_int(0) == 2);
+    static_assert(s.extent_int(1) == 3);
+    static_assert(s.extent_int(2) == 4);
+  }
+
+  {
+    ArrayND<double[2][3][4][5]> s;
+    static_assert(decltype(s)::rank == 4);
+    static_assert(s.extent_int(0) == 2);
+    static_assert(s.extent_int(1) == 3);
+    static_assert(s.extent_int(2) == 4);
+    static_assert(s.extent_int(3) == 5);
+  }
+
+  {
+    ArrayND<double[2][3][4][5][6]> s;
+    static_assert(decltype(s)::rank == 5);
+    static_assert(s.extent_int(0) == 2);
+    static_assert(s.extent_int(1) == 3);
+    static_assert(s.extent_int(2) == 4);
+    static_assert(s.extent_int(3) == 5);
+    static_assert(s.extent_int(4) == 6);
+  }
+}
+
+} // namespace sierra::nalu
diff --git a/unit_tests/matrix_free/CMakeLists.txt b/unit_tests/matrix_free/CMakeLists.txt
index 23e297d5ab..d5f0c8aa29 100644
--- a/unit_tests/matrix_free/CMakeLists.txt
+++ b/unit_tests/matrix_free/CMakeLists.txt
@@ -28,7 +28,6 @@ target_sources(${utest_ex_name} PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLinearExposedAreas.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLinearAdvectionMetric.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLinearVolume.C
-   ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLocalArray.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLocalDualNodalVolume.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLowMachUpdate.C
    ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestMatrixFreeSolver.C
diff --git a/unit_tests/matrix_free/UnitTestElementOperations.C b/unit_tests/matrix_free/UnitTestElementOperations.C
index 357ee146b1..2fb057049b 100644
--- a/unit_tests/matrix_free/UnitTestElementOperations.C
+++ b/unit_tests/matrix_free/UnitTestElementOperations.C
@@ -20,7 +20,7 @@
 #include "matrix_free/LobattoQuadratureRule.h"
 #include "matrix_free/PolynomialOrders.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "StkSimdComparisons.h"
 #include "gtest/gtest.h"
 
@@ -112,7 +112,7 @@ struct TensorPoly
 TEST(element_operations, scs_interp)
 {
   constexpr int p = inst::P2;
-  LocalArray<ftype[p + 1][p + 1][p + 1]> nodal_values;
+  ArrayND<ftype[p + 1][p + 1][p + 1]> nodal_values;
 
   auto poly = TensorPoly(p);
   for (int k = 0; k < p + 1; ++k) {
@@ -124,7 +124,7 @@ TEST(element_operations, scs_interp)
     }
   }
 
-  LocalArray<ftype[3][p + 1][p + 1][p + 1]> interp_values;
+  ArrayND<ftype[3][p + 1][p + 1][p + 1]> interp_values;
   interp_scs<p>(nodal_values, Coeffs<p>::Nt, interp_values);
 
   for (int k = 0; k < p + 1; ++k) {
@@ -161,7 +161,7 @@ TEST(element_operations, scs_interp)
 TEST(element_operations, integrate_volume)
 {
   constexpr int p = inst::P2;
-  LocalArray<ftype[p + 1][p + 1][p + 1]> nodal_values;
+  ArrayND<ftype[p + 1][p + 1][p + 1]> nodal_values;
 
   auto poly = TensorPoly(p);
   for (int k = 0; k < p + 1; ++k) {
@@ -180,7 +180,7 @@ TEST(element_operations, integrate_volume)
   }
   scs_end_loc[p + 1] = +1;
 
-  LocalArray<ftype[p + 1][p + 1][p + 1]> volumes;
+  ArrayND<ftype[p + 1][p + 1][p + 1]> volumes;
   for (int k = 0; k < p + 1; ++k) {
     for (int j = 0; j < p + 1; ++j) {
       for (int i = 0; i < p + 1; ++i) {
@@ -188,7 +188,7 @@ TEST(element_operations, integrate_volume)
       }
     }
   }
-  LocalArray<ftype[p + 1][p + 1][p + 1]> scratch;
+  ArrayND<ftype[p + 1][p + 1][p + 1]> scratch;
   volume<p>(nodal_values, scratch, volumes);
 
   for (int k = 0; k < p + 1; ++k) {
@@ -206,9 +206,9 @@ TEST(element_operations, integrate_volume)
 TEST(element_operations, nodal_grad)
 {
   constexpr int p = inst::P2;
-  LocalArray<ftype[p + 1][p + 1][p + 1]> nodal_values;
+  ArrayND<ftype[p + 1][p + 1][p + 1]> nodal_values;
 
-  LocalArray<ftype[p + 1][p + 1][p + 1][3]> xc;
+  ArrayND<ftype[p + 1][p + 1][p + 1][3]> xc;
 
   auto poly = TensorPoly(p);
   for (int k = 0; k < p + 1; ++k) {
@@ -246,9 +246,9 @@ TEST(element_operations, nodal_grad)
 TEST(element_operations, fp_grad)
 {
   constexpr int p = inst::P3;
-  LocalArray<ftype[p + 1][p + 1][p + 1]> nodal_values;
+  ArrayND<ftype[p + 1][p + 1][p + 1]> nodal_values;
 
-  LocalArray<ftype[p + 1][p + 1][p + 1][3]> xc;
+  ArrayND<ftype[p + 1][p + 1][p + 1][3]> xc;
 
   auto poly = TensorPoly(p);
   for (int k = 0; k < p + 1; ++k) {
@@ -267,7 +267,7 @@ TEST(element_operations, fp_grad)
 
   auto box = hex_vertex_coordinates<p>(xc);
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[p + 1][p + 1]> nhat;
+    ArrayND<ftype[p + 1][p + 1]> nhat;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         nhat(s, r) = interp_scs<p, 0>(nodal_values, l, s, r);
@@ -287,7 +287,7 @@ TEST(element_operations, fp_grad)
   }
 
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[p + 1][p + 1]> nhat;
+    ArrayND<ftype[p + 1][p + 1]> nhat;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         nhat(s, r) = interp_scs<p, 1>(nodal_values, l, s, r);
@@ -306,7 +306,7 @@ TEST(element_operations, fp_grad)
     }
   }
   for (int l = 0; l < p; ++l) {
-    LocalArray<ftype[p + 1][p + 1]> nhat;
+    ArrayND<ftype[p + 1][p + 1]> nhat;
     for (int s = 0; s < p + 1; ++s) {
       for (int r = 0; r < p + 1; ++r) {
         nhat(s, r) = interp_scs<p, 2>(nodal_values, l, s, r);
diff --git a/unit_tests/matrix_free/UnitTestLinearAdvectionMetric.C b/unit_tests/matrix_free/UnitTestLinearAdvectionMetric.C
index 3005b2cc8e..ab6ae9d103 100644
--- a/unit_tests/matrix_free/UnitTestLinearAdvectionMetric.C
+++ b/unit_tests/matrix_free/UnitTestLinearAdvectionMetric.C
@@ -14,8 +14,8 @@
 #include "matrix_free/LinearAreas.h"
 #include "matrix_free/LinearDiffusionMetric.h"
 #include "matrix_free/LobattoQuadratureRule.h"
-#include "matrix_free/LocalArray.h"
-#include "matrix_free/StkSimdComparisons.h"
+#include "ArrayND.h"
+#include "StkSimdComparisons.h"
 #include "matrix_free/TensorOperations.h"
 
 #include "Kokkos_Core.hpp"
@@ -31,7 +31,7 @@ void
 advection_single_cube_hex_p()
 {
   constexpr int poly = p;
-  LocalArray<double[3][3]> jac = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
+  ArrayND<double[3][3]> jac = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
 
   constexpr auto nodes = GLL<poly>::nodes;
   const int num_elems_1D = 32 / poly;
diff --git a/unit_tests/matrix_free/UnitTestLinearAreas.C b/unit_tests/matrix_free/UnitTestLinearAreas.C
index 7086db1571..ddfa5cd203 100644
--- a/unit_tests/matrix_free/UnitTestLinearAreas.C
+++ b/unit_tests/matrix_free/UnitTestLinearAreas.C
@@ -10,7 +10,7 @@
 #include "matrix_free/KokkosViewTypes.h"
 #include "matrix_free/LinearAreas.h"
 #include "matrix_free/LobattoQuadratureRule.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/TensorOperations.h"
 
 #include "StkSimdComparisons.h"
@@ -32,7 +32,7 @@ area_single_cube_hex_p()
 {
   constexpr int poly = p;
 
-  LocalArray<double[3][3]> jac = {
+  ArrayND<double[3][3]> jac = {
     {{+1.1, -2.6, 0}, {-1.2, .7, -0.2}, {10, -std::sqrt(3.), 12}}};
 
   constexpr auto nodes = GLL<poly>::nodes;
diff --git a/unit_tests/matrix_free/UnitTestLinearVolume.C b/unit_tests/matrix_free/UnitTestLinearVolume.C
index 6bbf9f8190..c4d2adf8dd 100644
--- a/unit_tests/matrix_free/UnitTestLinearVolume.C
+++ b/unit_tests/matrix_free/UnitTestLinearVolume.C
@@ -13,7 +13,7 @@
 #include "matrix_free/LobattoQuadratureRule.h"
 #include "matrix_free/TensorOperations.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 
 #include <Kokkos_Core.hpp>
 #include <stk_simd/Simd.hpp>
@@ -28,7 +28,7 @@ template <int poly>
 void
 single_affine_hex_p(bool cube)
 {
-  LocalArray<double[3][3]> jac = {{{0, 0, 1}, {1, 0, 0}, {0, 1, 0}}};
+  ArrayND<double[3][3]> jac = {{{0, 0, 1}, {1, 0, 0}, {0, 1, 0}}};
   if (!cube) {
     jac = {{{2, 1, 1.3333}, {0, 2, -1}, {1, 0, 2}}};
   }
diff --git a/unit_tests/matrix_free/UnitTestLocalDualNodalVolume.C b/unit_tests/matrix_free/UnitTestLocalDualNodalVolume.C
index 3ec4d7ed42..6022e1288d 100644
--- a/unit_tests/matrix_free/UnitTestLocalDualNodalVolume.C
+++ b/unit_tests/matrix_free/UnitTestLocalDualNodalVolume.C
@@ -14,7 +14,7 @@
 #include "matrix_free/LowMachFields.h"
 #include "matrix_free/StkLowMachFixture.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/StkSimdConnectivityMap.h"
 #include "matrix_free/TensorOperations.h"
 
@@ -34,7 +34,7 @@ skew_mesh(
   const stk::mesh::Selector& sel,
   stk::mesh::NgpField<double>& coords)
 {
-  LocalArray<double[3][3]> Q = {{{2, 1, 1.3333}, {0, 2, -1}, {1, 0, 2}}};
+  ArrayND<double[3][3]> Q = {{{2, 1, 1.3333}, {0, 2, -1}, {1, 0, 2}}};
   double detq = determinant(Q);
 
   stk::mesh::for_each_entity_run(
diff --git a/unit_tests/matrix_free/UnitTestTransportCoefficients.C b/unit_tests/matrix_free/UnitTestTransportCoefficients.C
index 489ce25a95..e31425462f 100644
--- a/unit_tests/matrix_free/UnitTestTransportCoefficients.C
+++ b/unit_tests/matrix_free/UnitTestTransportCoefficients.C
@@ -14,7 +14,7 @@
 #include "matrix_free/LowMachFields.h"
 #include "matrix_free/StkLowMachFixture.h"
 #include "matrix_free/KokkosViewTypes.h"
-#include "matrix_free/LocalArray.h"
+#include "ArrayND.h"
 #include "matrix_free/StkSimdConnectivityMap.h"
 #include "matrix_free/TensorOperations.h"
 
@@ -39,7 +39,7 @@ public:
 
   void skew_mesh()
   {
-    LocalArray<double[3][3]> Q = {{{2, 1, 1.3333}, {0, 2, -1}, {1, 0, 2}}};
+    ArrayND<double[3][3]> Q = {{{2, 1, 1.3333}, {0, 2, -1}, {1, 0, 2}}};
     auto coords = stk::mesh::get_updated_ngp_field<double>(coordinate_field());
 
     stk::mesh::for_each_entity_run(