diff --git a/src/DataStructures/CMakeLists.txt b/src/DataStructures/CMakeLists.txt
index 2dd61811c762..8abc889e82eb 100644
--- a/src/DataStructures/CMakeLists.txt
+++ b/src/DataStructures/CMakeLists.txt
@@ -8,6 +8,7 @@ set(LIBRARY_SOURCES
     Index.cpp
     IndexIterator.cpp
     Mesh.cpp
+    ModalVector.cpp
     SliceIterator.cpp
     StripeIterator.cpp
     VariablesHelpers.cpp
diff --git a/src/DataStructures/ModalVector.cpp b/src/DataStructures/ModalVector.cpp
new file mode 100644
index 000000000000..bf1e1d2be2d7
--- /dev/null
+++ b/src/DataStructures/ModalVector.cpp
@@ -0,0 +1,132 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "DataStructures/ModalVector.hpp"
+
+#include <algorithm>
+#include <pup.h>
+#include <pup_stl.h>
+
+#include "Utilities/StdHelpers.hpp"
+
+ModalVector::ModalVector(const size_t size, const double value)
+    : size_(size),
+      owned_data_(size_, value),
+      data_(owned_data_.data(), size_) {}
+
+template <class T, Requires<cpp17::is_same_v<T, double>>>
+ModalVector::ModalVector(std::initializer_list<T> list)
+    : size_(list.size()),
+      owned_data_(std::move(list)),
+      data_(owned_data_.data(), size_) {}
+
+ModalVector::ModalVector(double* start, size_t size)
+    : size_(size), owned_data_(0), data_(start, size_), owning_(false) {}
+
+/// \cond HIDDEN_SYMBOLS
+ModalVector::ModalVector(const ModalVector& rhs) {
+  size_ = rhs.size();
+  if (rhs.is_owning()) {
+    owned_data_ = rhs.owned_data_;
+  } else {
+    owned_data_ = InternalStorage_t(rhs.begin(), rhs.end());
+  }
+  data_ = decltype(data_){owned_data_.data(), size_};
+}
+
+ModalVector& ModalVector::operator=(const ModalVector& rhs) {
+  if (this == &rhs) {
+    return *this;
+  }
+  if (owning_) {
+    size_ = rhs.size();
+    if (rhs.is_owning()) {
+      owned_data_ = rhs.owned_data_;
+    } else {
+      owned_data_ = InternalStorage_t(rhs.begin(), rhs.end());
+    }
+    data_ = decltype(data_){owned_data_.data(), size_};
+  } else {
+    ASSERT(rhs.size() == size(), "Must copy into same size, not "
+                                     << rhs.size() << " into " << size());
+    std::copy(rhs.begin(), rhs.end(), begin());
+  }
+  return *this;
+}
+
+ModalVector::ModalVector(ModalVector&& rhs) noexcept {
+  size_ = rhs.size_;
+  owned_data_ = std::move(rhs.owned_data_);
+  // clang-tidy: move trivially copyable type, future proof in case impl
+  // changes
+  data_ = std::move(rhs.data_);  // NOLINT
+  owning_ = rhs.owning_;
+
+  rhs.owning_ = true;
+  rhs.size_ = 0;
+  rhs.data_ = decltype(rhs.data_){};
+}
+
+ModalVector&
+ModalVector::operator=(ModalVector&& rhs) noexcept {
+  if (this == &rhs) {
+    return *this;
+  }
+  if (owning_) {
+    size_ = rhs.size_;
+    owned_data_ = std::move(rhs.owned_data_);
+    // clang-tidy: move trivially copyable type, future proof in case impl
+    // changes
+    data_ = std::move(rhs.data_);  // NOLINT
+    owning_ = rhs.owning_;
+  } else {
+    ASSERT(rhs.size() == size(), "Must copy into same size, not "
+                                     << rhs.size() << " into " << size());
+    std::copy(rhs.begin(), rhs.end(), begin());
+  }
+  rhs.owning_ = true;
+  rhs.size_ = 0;
+  rhs.data_ = decltype(rhs.data_){};
+  return *this;
+}
+/// \endcond
+
+void ModalVector::pup(PUP::er& p) noexcept {  // NOLINT
+  p | size_;
+  if (p.isUnpacking()) {
+    owning_ = true;
+    p | owned_data_;
+    data_ = decltype(data_){owned_data_.data(), size_};
+  } else {
+    if (not owning_) {
+      owned_data_ =
+          InternalStorage_t(data_.data(), data_.data() + size_);  // NOLINT
+      p | owned_data_;
+      owned_data_.clear();
+    } else {
+      p | owned_data_;
+    }
+  }
+}
+
+std::ostream& operator<<(std::ostream& os, const ModalVector& d) {
+  // This function is inside the detail namespace StdHelpers.hpp
+  StdHelpers_detail::print_helper(os, d.begin(), d.end());
+  return os;
+}
+
+/// Equivalence operator for ModalVector
+bool operator==(const ModalVector& lhs, const ModalVector& rhs) {
+  return lhs.size() == rhs.size() and
+         std::equal(lhs.begin(), lhs.end(), rhs.begin());
+}
+
+/// Inequivalence operator for ModalVector
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs) {
+  return not(lhs == rhs);
+}
+
+/// \cond
+template
+ModalVector::ModalVector(std::initializer_list<double> list);
+/// \endcond
diff --git a/src/DataStructures/ModalVector.hpp b/src/DataStructures/ModalVector.hpp
new file mode 100644
index 000000000000..443d069fa091
--- /dev/null
+++ b/src/DataStructures/ModalVector.hpp
@@ -0,0 +1,490 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+/// \file
+/// Defines class ModalVector.
+
+#pragma once
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <initializer_list>
+#include <limits>
+#include <ostream>
+#include <type_traits>
+#include <vector>
+
+#include "DataStructures/DataVector.hpp"
+#include "ErrorHandling/Assert.hpp"
+#include "Utilities/ForceInline.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/PointerVector.hpp" // IWYU pragma: keep
+#include "Utilities/Requires.hpp"
+
+/// \cond HIDDEN_SYMBOLS
+namespace PUP {
+class er;
+}  // namespace PUP
+
+// clang-tidy: no using declarations in header files
+//             We want the std::abs to be used
+using std::abs;  // NOLINT
+/// \endcond
+
+/*!
+ * \ingroup DataStructuresGroup
+ * \brief A class for storing spectral coefficients on a mesh.
+ *
+ * A ModalVector holds an array of spectral coefficients, and can be
+ * either owning (the array is deleted when the ModalVector goes out of scope)
+ * or non-owning, meaning it just has a pointer to an array.
+ *
+ * Only basic mathematical operations are supported with ModalVectors. In
+ * addition to addition, subtraction, multiplication, division, there
+ * are the following element-wise operations:
+ *
+ * - abs
+ * - max
+ * - min
+ *
+ * In order to allow filtering, multiplication (*, *=) and division (/, /=)
+ * operations with a DataVectors (holding filters) is supported.
+ *
+ * Note : Because of no tagging feature yet, expressions containing purely
+ * DataVectors will evaluate and get copy-constructed to a ModalVector.
+ * Disallowing this would also mean disallowing expressions such as
+ * ModalVector coeffs_filtered(coeffs * filter) - where:
+ *  - coeffs is a ModalVector
+ *  - filter is a DataVector
+ */
+class ModalVector {
+  /// \cond HIDDEN_SYMBOLS
+  static constexpr int private_asserts() {
+    static_assert(std::is_nothrow_move_constructible<ModalVector>::value,
+                  "Missing move semantics");
+    return 0;
+  }
+  /// \endcond
+ public:
+  using value_type = double;
+  using allocator_type = std::allocator<value_type>;
+  using size_type = size_t;
+  using difference_type = std::ptrdiff_t;
+
+  // The type alias ElementType is needed because blaze::IsInvertible<T> is not
+  // SFINAE friendly to an ElementType type alias
+  using ElementType = double;
+
+ private:
+  /// The type of the "pointer" used internally
+  using InternalModalVector_t = PointerVector<double>;
+  /// The type used to store the data in
+  using InternalStorage_t = std::vector<double, allocator_type>;
+
+ public:
+  /// Create with the given size and value.
+  ///
+  /// \param size number of values
+  /// \param value the value to initialize each element.
+  explicit
+  ModalVector(size_t size,
+              double value = std::numeric_limits<double>::signaling_NaN());
+
+  /// Create a non-owning ModalVector that points to `start`
+  ModalVector(double* start, size_t size);
+
+  /// Create from an initializer list of doubles. All elements in the
+  /// `std::initializer_list` must have decimal points
+  // cppcheck-suppress syntaxError
+  template <class T, Requires<cpp17::is_same_v<T, double>> = nullptr>
+  ModalVector(std::initializer_list<T> list);
+
+  /// Empty ModalVector
+  ModalVector() = default;
+  /// \cond HIDDEN_SYMBOLS
+  ~ModalVector() = default;
+
+  ModalVector(const ModalVector& rhs);
+  ModalVector(ModalVector&& rhs) noexcept;
+  ModalVector& operator=(const ModalVector& rhs);
+  ModalVector& operator=(ModalVector&& rhs) noexcept;
+
+  // This is a converting constructor. clang-tidy complains that it's not
+  // explicit, but we want it to allow conversion.
+  // clang-tidy: mark as explicit (we want conversion to ModalVector)
+  template <typename VT, bool VF>
+  ModalVector(const blaze::Vector<VT, VF>& expression);  // NOLINT
+
+  template <typename VT, bool VF>
+  ModalVector& operator=(const blaze::Vector<VT, VF>& expression);
+  /// \endcond
+
+  /// Number of values stored
+  size_t size() const noexcept { return size_; }
+
+  // @{
+  /// Set the ModalVector to be a reference to another ModalVector
+  /// object
+  void set_data_ref(gsl::not_null<ModalVector*> rhs) noexcept {
+    set_data_ref(rhs->data(), rhs->size_);
+  }
+  void set_data_ref(double* start, size_t size) noexcept {
+    size_ = size;
+    owned_data_ = decltype(owned_data_){};
+    data_ = decltype(data_){start, size_};
+    owning_ = false;
+  }
+  // @}
+
+  /// Returns true if the class owns the data
+  bool is_owning() const noexcept { return owning_; }
+
+  // @{
+  /// Access ith element
+  double& operator[](const size_type i) noexcept {
+    ASSERT(i < size_, "i = " << i << ", size = " << size_);
+    // clang-tidy: do not use pointer arithmetic
+    return data_[i];  // NOLINT
+  }
+  const double& operator[](const size_type i) const noexcept {
+    ASSERT(i < size_, "i = " << i << ", size = " << size_);
+    // clang-tidy: do not use pointer arithmetic
+    return data_[i];  // NOLINT
+  }
+  // @}
+
+  // @{
+  /// Access to the pointer
+  double* data() noexcept { return data_.data(); }
+  const double* data() const noexcept { return data_.data(); }
+  // @}
+
+  // @{
+  /// Returns iterator to beginning of data
+  decltype(auto) begin() noexcept { return data_.begin(); }
+  decltype(auto) begin() const noexcept { return data_.begin(); }
+  // @}
+  // @{
+  /// Returns iterator to end of data
+  decltype(auto) end() noexcept { return data_.end(); }
+  decltype(auto) end() const noexcept { return data_.end(); }
+  // @}
+
+  /// Serialization for Charm++
+  // clang-tidy: google-runtime-references
+  void pup(PUP::er& p) noexcept;  // NOLINT
+
+  // @{
+  /// See the Blaze library documentation for details on these functions since
+  /// they merely forward to Blaze.
+  /// https://bitbucket.org/blaze-lib/blaze/overview
+  ModalVector& operator=(const double& rhs) noexcept {
+    data_ = rhs;
+    return *this;
+  }
+
+  ModalVector& operator+=(const ModalVector& rhs) noexcept {
+    data_ += rhs.data_;
+    return *this;
+  }
+  template <typename VT, bool VF>
+  ModalVector& operator+=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ += rhs;
+    return *this;
+  }
+  ModalVector& operator+=(const double& rhs) noexcept {
+    data_ += rhs;
+    return *this;
+  }
+
+  ModalVector& operator-=(const ModalVector& rhs) noexcept {
+    data_ -= rhs.data_;
+    return *this;
+  }
+  template <typename VT, bool VF>
+  ModalVector& operator-=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ -= rhs;
+    return *this;
+  }
+  ModalVector& operator-=(const double& rhs) noexcept {
+    data_ -= rhs;
+    return *this;
+  }
+
+  template <typename VT, bool VF>
+  ModalVector& operator*=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ *= rhs;
+    return *this;
+  }
+  ModalVector& operator*=(const double& rhs) noexcept {
+    data_ *= rhs;
+    return *this;
+  }
+  ModalVector& operator*=(const ModalVector& rhs) noexcept {
+    data_ *= rhs.data_;
+    return *this;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV object in this function
+  ModalVector& operator*=(const DataVector& rhs) noexcept {
+    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    data_ *= _rhs_local;
+    return *this;
+  }
+
+  template <typename VT, bool VF>
+  ModalVector& operator/=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ /= rhs;
+    return *this;
+  }
+  ModalVector& operator/=(const double& rhs) noexcept {
+    data_ /= rhs;
+    return *this;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV object in this function
+  ModalVector& operator/=(const DataVector& rhs) noexcept {
+    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    data_ /= _rhs_local;
+    return *this;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const ModalVector& lhs, const ModalVector& rhs) noexcept {
+    return lhs.data_ + rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
+    return ~lhs + rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ + ~rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ + rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const double& lhs, const ModalVector& rhs) noexcept {
+    return lhs + rhs.data_;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& lhs, const ModalVector& rhs) noexcept {
+    return lhs.data_ - rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
+    return ~lhs - rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ - ~rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ - rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const double& lhs, const ModalVector& rhs) noexcept {
+    return lhs - rhs.data_;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& rhs) noexcept {
+    return -rhs.data_;
+  };
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ * rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const double& lhs, const ModalVector& rhs) noexcept {
+    return lhs * rhs.data_;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const ModalVector& rhs) noexcept {
+    return lhs.data_ * rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
+    return ~lhs * rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ * ~rhs;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV objects in the
+  // next 2 functions?
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const DataVector& rhs) noexcept {
+    blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    return lhs * _rhs_local;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const DataVector& lhs, const ModalVector& rhs) noexcept {
+    blaze::DynamicVector<double> _lhs_local(lhs.size(), lhs.data());
+    return rhs * _lhs_local;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ / rhs;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ / ~rhs;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV object in this function
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
+      const ModalVector& lhs, const DataVector& rhs) noexcept {
+    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    return lhs / _rhs_local;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) min(
+      const ModalVector& t) noexcept {
+    return min(t.data_);
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) max(
+      const ModalVector& t) noexcept {
+    return max(t.data_);
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) abs(
+      const ModalVector& t) noexcept {
+    return abs(t.data_);
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) fabs(
+      const ModalVector& t) noexcept {
+    return abs(t.data_);
+  }
+  // @}
+
+
+ private:
+  /// \cond HIDDEN_SYMBOLS
+  size_t size_ = 0;
+  InternalStorage_t owned_data_;
+  InternalModalVector_t data_;
+  bool owning_{true};
+  /// \endcond
+};
+
+/// Output operator for ModalVector
+std::ostream& operator<<(std::ostream& os, const ModalVector& d);
+
+/// Equivalence operator for ModalVector
+bool operator==(const ModalVector& lhs, const ModalVector& rhs);
+
+/// Inequivalence operator for ModalVector
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs);
+
+template <typename T, size_t Dim>
+std::array<ModalVector, Dim> operator+(
+    const std::array<T, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) + gsl::at(rhs, i);
+  }
+  return result;
+}
+template <typename U, size_t Dim>
+std::array<ModalVector, Dim> operator+(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<U, Dim>& rhs) noexcept {
+  return rhs + lhs;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim> operator+(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) + gsl::at(rhs, i);
+  }
+  return result;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim>& operator+=(
+    std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(lhs, i) += gsl::at(rhs, i);
+  }
+  return lhs;
+}
+
+template <typename T, size_t Dim>
+std::array<ModalVector, Dim> operator-(
+    const std::array<T, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) - gsl::at(rhs, i);
+  }
+  return result;
+}
+template <typename U, size_t Dim>
+std::array<ModalVector, Dim> operator-(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<U, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) - gsl::at(rhs, i);
+  }
+  return result;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim> operator-(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) - gsl::at(rhs, i);
+  }
+  return result;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim>& operator-=(
+    std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(lhs, i) -= gsl::at(rhs, i);
+  }
+  return lhs;
+}
+
+/// \cond HIDDEN_SYMBOLS
+template <typename VT, bool VF>
+ModalVector::ModalVector(const blaze::Vector<VT, VF>& expression)
+    : size_((~expression).size()),
+      owned_data_((~expression).size()),
+      data_(owned_data_.data(), (~expression).size()) {
+  data_ = expression;
+}
+
+template <typename VT, bool VF>
+ModalVector& ModalVector::operator=(
+    const blaze::Vector<VT, VF>& expression) {
+  if (owning_ and (~expression).size() != size()) {
+    size_ = (~expression).size();
+    owned_data_ = InternalStorage_t(size_);
+    data_ = decltype(data_){owned_data_.data(), size_};
+  } else if (not owning_) {
+    ASSERT((~expression).size() == size(), "Must copy into same size");
+  }
+  data_ = expression;
+  return *this;
+}
+/// \endcond
diff --git a/tests/Unit/DataStructures/CMakeLists.txt b/tests/Unit/DataStructures/CMakeLists.txt
index 6c3dfd813afa..c90f8b126982 100644
--- a/tests/Unit/DataStructures/CMakeLists.txt
+++ b/tests/Unit/DataStructures/CMakeLists.txt
@@ -17,6 +17,7 @@ set(LIBRARY_SOURCES
   Test_Index.cpp
   Test_IndexIterator.cpp
   Test_Mesh.cpp
+  Test_ModalVector.cpp
   Test_OrientVariablesOnSlice.cpp
   Test_SliceIterator.cpp
   Test_StripeIterator.cpp
diff --git a/tests/Unit/DataStructures/Test_ModalVector.cpp b/tests/Unit/DataStructures/Test_ModalVector.cpp
new file mode 100644
index 000000000000..47a2e3334a4b
--- /dev/null
+++ b/tests/Unit/DataStructures/Test_ModalVector.cpp
@@ -0,0 +1,455 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "tests/Unit/TestingFramework.hpp"
+
+#include <iostream>
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <numeric>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/ModalVector.hpp"
+#include "ErrorHandling/Error.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/StdArrayHelpers.hpp"
+#include "Utilities/StdHelpers.hpp"  // IWYU pragma: keep
+#include "tests/Unit/TestHelpers.hpp"
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector",
+                  "[DataStructures][Unit]") {
+  ModalVector a{2};
+  CHECK(a.size() == 2);
+  ModalVector b{2, 10.0};
+  CHECK(b.size() == 2);
+  for (size_t i = 0; i < b.size(); ++i) {
+    INFO(i);
+    CHECK(b[i] == 10.0);
+  }
+
+  ModalVector t(5, 10.0);
+  CHECK(t.size() == 5);
+  for (size_t i = 0; i < t.size(); ++i) {
+    INFO(i);
+    CHECK(t[i] == 10.0);
+  }
+  for (const auto& p : t) {
+    CHECK(p == 10.0);
+  }
+  for (auto& p : t) {
+    CHECK(p == 10.0);
+  }
+  ModalVector t2{1.43, 2.83, 3.94, 7.85};
+  CHECK(t2.size() == 4);
+  CHECK(t2.is_owning());
+  CHECK(t2[0] == 1.43);
+  CHECK(t2[1] == 2.83);
+  CHECK(t2[2] == 3.94);
+  CHECK(t2[3] == 7.85);
+  test_copy_semantics(t);
+  auto t_copy = t;
+  CHECK(t_copy.is_owning());
+  test_move_semantics(std::move(t), t_copy);
+  ModalVector t_move_assignment = std::move(t_copy);
+  CHECK(t_move_assignment.is_owning());
+  ModalVector t_move_constructor = std::move(t_move_assignment);
+  CHECK(t_move_constructor.is_owning());
+}
+
+SPECTRE_TEST_CASE("Unit.Serialization.ModalVector",
+                  "[DataStructures][Unit][Serialization]") {
+  const size_t npts = 10;
+  ModalVector t(npts), tgood(npts);
+  std::iota(t.begin(), t.end(), 1.2);
+  std::iota(tgood.begin(), tgood.end(), 1.2);
+  CHECK(tgood == t);
+  CHECK(t.is_owning());
+  CHECK(tgood.is_owning());
+  const ModalVector serialized = serialize_and_deserialize(t);
+  CHECK(tgood == t);
+  CHECK(serialized == tgood);
+  CHECK(serialized.is_owning());
+  CHECK(serialized.data() != t.data());
+  CHECK(t.is_owning());
+}
+
+SPECTRE_TEST_CASE("Unit.Serialization.ModalVector_Ref",
+                  "[DataStructures][Unit][Serialization]") {
+  const size_t npts = 10;
+  ModalVector t(npts);
+  std::iota(t.begin(), t.end(), 4.3);
+  ModalVector t2;
+  t2.set_data_ref(&t);
+  CHECK(t.is_owning());
+  CHECK_FALSE(t2.is_owning());
+  CHECK(t2 == t);
+  const ModalVector serialized = serialize_and_deserialize(t);
+  CHECK(t2 == t);
+  CHECK(serialized == t2);
+  CHECK(serialized.is_owning());
+  CHECK(serialized.data() != t.data());
+  CHECK(t.is_owning());
+  const ModalVector serialized2 = serialize_and_deserialize(t2);
+  CHECK(t2 == t);
+  CHECK(serialized2 == t2);
+  CHECK(serialized2.is_owning());
+  CHECK(serialized2.data() != t2.data());
+  CHECK_FALSE(t2.is_owning());
+}
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
+                  "[DataStructures][Unit]") {
+  ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector t;
+  t.set_data_ref(&data);
+  CHECK(not t.is_owning());
+  CHECK(data.is_owning());
+  CHECK(t.data() == data.data());
+  CHECK(t.size() == 4);
+  CHECK(t[0] == 1.43);
+  CHECK(t[1] == 2.83);
+  CHECK(t[2] == 3.94);
+  CHECK(t[3] == 7.85);
+  test_copy_semantics(t);
+  ModalVector t_copy;
+  t_copy.set_data_ref(&t);
+  test_move_semantics(std::move(t), t_copy);
+  ModalVector t_move_assignment = std::move(t_copy);
+  CHECK(not t_move_assignment.is_owning());
+  ModalVector t_move_constructor = std::move(t_move_assignment);
+  CHECK(not t_move_constructor.is_owning());
+  {
+    ModalVector data_2{1.43, 2.83, 3.94, 7.85};
+    ModalVector data_2_ref;
+    data_2_ref.set_data_ref(&data_2);
+    ModalVector data_3{2.43, 3.83, 4.94, 8.85};
+    data_2_ref = std::move(data_3);
+    CHECK(data_2[0] == 2.43);
+    CHECK(data_2[1] == 3.83);
+    CHECK(data_2[2] == 4.94);
+    CHECK(data_2[3] == 8.85);
+// Intentionally testing self-move
+#ifdef __clang__
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wself-move"
+#endif  // defined(__clang__)
+    data_2_ref = std::move(data_2_ref);
+#ifdef __clang__
+#pragma GCC diagnostic pop
+#endif  // defined(__clang__)
+    CHECK(data_2[0] == 2.43);
+    CHECK(data_2[1] == 3.83);
+    CHECK(data_2[2] == 4.94);
+    CHECK(data_2[3] == 8.85);
+    ModalVector owned_data;
+    // clang-tidy: false positive, used after it was moved
+    owned_data = data_2_ref;  // NOLINT
+    CHECK(owned_data[0] == 2.43);
+    CHECK(owned_data[1] == 3.83);
+    CHECK(owned_data[2] == 4.94);
+    CHECK(owned_data[3] == 8.85);
+    // Test operator!=
+    CHECK_FALSE(owned_data != data_2_ref);
+  }
+}
+
+// [[OutputRegex, Must copy into same size]]
+[[noreturn]] SPECTRE_TEST_CASE(
+                "Unit.DataStructures.ModalVector.ref_diff_size",
+                "[DataStructures][Unit]") {
+  ASSERTION_TEST();
+#ifdef SPECTRE_DEBUG
+  ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector data_ref;
+  data_ref.set_data_ref(&data);
+  ModalVector data2{1.43, 2.83, 3.94};
+  data_ref = data2;
+  ERROR("Failed to trigger ASSERT in an assertion test");
+#endif
+}
+
+// [[OutputRegex, Must copy into same size]]
+[[noreturn]] SPECTRE_TEST_CASE(
+    "Unit.DataStructures.ModalVector.move_ref_diff_size",
+    "[DataStructures][Unit]") {
+  ASSERTION_TEST();
+#ifdef SPECTRE_DEBUG
+  ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector data_ref;
+  data_ref.set_data_ref(&data);
+  ModalVector data2{1.43, 2.83, 3.94};
+  data_ref = std::move(data2);
+  ERROR("Failed to trigger ASSERT in an assertion test");
+#endif
+}
+
+namespace {
+template <typename T1, typename T2>
+void check_vectors(const T1& t1, const T2& t2) {
+  CHECK_ITERABLE_APPROX(t1, t2);
+}
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathAfterMove",
+                  "[Unit][DataStructures]") {
+  ModalVector m0(10, 3.0), m1(10, 9.0);
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{};
+    b = std::move(a);
+    b = m0 + m1;
+    check_vectors(b, ModalVector(10, 12.0));
+    // clang-tidy: use after move (intentional here)
+    CHECK(a.size() == 0);  // NOLINT
+    CHECK(a.is_owning());
+    a = m0 * m1;
+    check_vectors(a, ModalVector(10, 27.0));
+    check_vectors(b, ModalVector(10, 12.0));
+  }
+
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{};
+    b = std::move(a);
+    a = m0 + m1;
+    check_vectors(b, ModalVector(10, 2.0));
+    check_vectors(a, ModalVector(10, 12.0));
+  }
+
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{std::move(a)};
+    b = m0 + m1;
+    CHECK(b.size() == 10);
+    check_vectors(b, ModalVector(10, 12.0));
+    // clang-tidy: use after move (intentional here)
+    CHECK(a.size() == 0);  // NOLINT
+    CHECK(a.is_owning());
+    a = m0 * m1;
+    check_vectors(a, ModalVector(10, 27.0));
+    check_vectors(b, ModalVector(10, 12.0));
+  }
+
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{std::move(a)};
+    a = m0 + m1;
+    check_vectors(b, ModalVector(10, 2.0));
+    check_vectors(a, ModalVector(10, 12.0));
+  }
+}
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
+                  "[Unit][DataStructures]") {
+  const size_t num_pts = 19;
+  ModalVector val{1., -2., 3., -4., 8., 12., -14.};
+  ModalVector nine(num_pts, 9.0);
+  ModalVector one(num_pts, 1.0);
+
+  // Test unary minus
+  check_vectors(-nine, ModalVector(num_pts, -9.0));
+
+  check_vectors(nine + 2.0, ModalVector(num_pts, 11.0));
+  check_vectors(2.0 + nine, ModalVector(num_pts, 11.0));
+  check_vectors(nine - 2.0, ModalVector(num_pts, 7.0));
+  check_vectors(2.0 - nine, ModalVector(num_pts, -7.0));
+  check_vectors(nine + nine, ModalVector(num_pts, 18.0));
+  check_vectors(nine + (one * nine), ModalVector(num_pts, 18.0));
+  check_vectors((one * nine) + nine, ModalVector(num_pts, 18.0));
+  check_vectors(nine - ModalVector(num_pts, 8.0), one);
+  check_vectors(nine - (one * nine), ModalVector(num_pts, 0.0));
+  check_vectors((one * nine) - nine, ModalVector(num_pts, 0.0));
+
+  check_vectors(ModalVector(num_pts, 18.0), (one / 0.5) * nine);
+
+  CHECK(-14 == min(val));
+  CHECK(12 == max(val));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, abs(val));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, fabs(val));
+
+  check_vectors(ModalVector(num_pts, 81.0), nine * nine);
+  check_vectors(ModalVector(num_pts, 81.0), nine * (nine * one));
+  check_vectors(ModalVector(num_pts, 81.0), (nine * nine) * one);
+  check_vectors(ModalVector(num_pts, 81.0), 9.0 * nine);
+  check_vectors(ModalVector(num_pts, 81.0), nine * 9.0);
+  check_vectors(ModalVector(num_pts, 1.0), nine / 9.0);
+
+  ModalVector dummy(nine * nine * 1.0);
+  check_vectors(ModalVector(num_pts, 81.0), dummy);
+
+  // Test assignment
+  ModalVector test_81(num_pts, -1.0);
+  test_81 = nine * nine;
+  check_vectors(ModalVector(num_pts, 81.0), test_81);
+  CHECK(test_81.is_owning());
+  ModalVector test_81_ref(test_81.data(), test_81.size());
+  test_81_ref = 0.0;
+  test_81 = 0.0;
+  test_81_ref = nine * nine;
+  check_vectors(ModalVector(num_pts, 81.0), test_81);
+  CHECK(test_81.is_owning());
+  check_vectors(ModalVector(num_pts, 81.0), test_81_ref);
+  CHECK_FALSE(test_81_ref.is_owning());
+  ModalVector second_81(num_pts);
+  second_81 = test_81;
+  check_vectors(ModalVector(num_pts, 81.0), second_81);
+  CHECK(second_81.is_owning());
+  test_81_ref = 0.0;
+  check_vectors(ModalVector(num_pts, 0.0), test_81_ref);
+  test_81_ref = second_81;
+  check_vectors(ModalVector(num_pts, 81.0), test_81_ref);
+  second_81 = 0.0;
+  test_81_ref.set_data_ref(&test_81);
+  second_81 = std::move(test_81_ref);
+  check_vectors(ModalVector(num_pts, 81.0), second_81);
+  CHECK_FALSE(second_81.is_owning());
+  second_81 = 0.0;
+  check_vectors(ModalVector(num_pts, 0.0), second_81);
+  test_81 = 81.0;
+  check_vectors(ModalVector(num_pts, 81.0), second_81);
+  CHECK_FALSE(second_81.is_owning());
+
+  test_81 = 81.0;
+  ModalVector test_081;
+  test_081.set_data_ref(&test_81);
+  check_vectors(ModalVector(num_pts, 81.0), test_081);
+  CHECK_FALSE(test_081.is_owning());
+
+  ModalVector test_assignment(num_pts, 7.0);
+  test_assignment += nine;
+  check_vectors(ModalVector(num_pts, 16.0), test_assignment);
+  test_assignment += 3.0;
+  check_vectors(ModalVector(num_pts, 19.0), test_assignment);
+  test_assignment += (nine * nine);
+  check_vectors(ModalVector(num_pts, 100.0), test_assignment);
+
+  test_assignment = 7.0;
+  test_assignment -= nine;
+  check_vectors(ModalVector(num_pts, -2.0), test_assignment);
+  test_assignment -= 3.0;
+  check_vectors(ModalVector(num_pts, -5.0), test_assignment);
+  test_assignment -= nine * nine;
+  check_vectors(ModalVector(num_pts, -86.0), test_assignment);
+
+  test_assignment = 2.0;
+  test_assignment *= 3.0;
+  check_vectors(ModalVector(num_pts, 6.0), test_assignment);
+  test_assignment *= nine;
+  check_vectors(ModalVector(num_pts, 54.0), test_assignment);
+  test_assignment = 1.0;
+  test_assignment *= nine * nine;
+  check_vectors(ModalVector(num_pts, 81.0), test_assignment);
+
+  test_assignment = 2.0;
+  test_assignment /= 2.0;
+  check_vectors(ModalVector(num_pts, 1.0), test_assignment);
+
+  // Test assignment where the RHS is an expression that contains the LHS
+  ModalVector x(num_pts, 4.);
+  x *= x;
+  check_vectors(ModalVector(num_pts, 16.0), x);
+
+  // Test addition of arrays of ModalVectors to arrays of doubles.
+  const ModalVector t1{0.0, 1.0, 2.0, 3.0};
+  const ModalVector t2{-0.1, -0.2, -0.3, -0.4};
+  const ModalVector t3{5.0, 4.0, 3.0, 2.0};
+  const ModalVector e1{10.0, 11.0, 12.0, 13.0};
+  const ModalVector e2{19.9, 19.8, 19.7, 19.6};
+  const ModalVector e3{35.0, 34.0, 33.0, 32.0};
+  const ModalVector e4{-10.0, -9.0, -8.0, -7.0};
+  const ModalVector e5{-20.1, -20.2, -20.3, -20.4};
+  const ModalVector e6{-25.0, -26.0, -27.0, -28.0};
+  const ModalVector e7{10.0, 12.0, 14.0, 16.0};
+  const ModalVector e8{19.8, 19.6, 19.4, 19.2};
+  const ModalVector e9{40.0, 38.0, 36.0, 34.0};
+  const ModalVector e10{-10.0, -10.0, -10.0, -10.0};
+  const ModalVector e11{-20.0, -20.0, -20.0, -20.0};
+  const ModalVector e12{-30.0, -30.0, -30.0, -30.0};
+
+  const std::array<double, 3> point{{10.0, 20.0, 30.0}};
+  std::array<ModalVector, 3> points{{t1, t2, t3}};
+  const std::array<ModalVector, 3> expected{{e1, e2, e3}};
+  const std::array<ModalVector, 3> expected2{{e4, e5, e6}};
+  std::array<ModalVector, 3> dvectors1{{t1, t2, t3}};
+  const std::array<ModalVector, 3> dvectors2{{e1, e2, e3}};
+  const std::array<ModalVector, 3> expected3{{e7, e8, e9}};
+  const std::array<ModalVector, 3> expected4{{e10, e11, e12}};
+  CHECK(points + point == expected);
+  CHECK(point + points == expected);
+  CHECK(points - point == expected2);
+  CHECK(point - points == -expected2);
+
+  points += point;
+  CHECK(points == expected);
+  points -= point;
+  points -= point;
+  CHECK(points == expected2);
+
+  CHECK_ITERABLE_APPROX(dvectors1 + dvectors2, expected3);
+  CHECK_ITERABLE_APPROX(dvectors1 - dvectors2, expected4);
+
+  dvectors1 += dvectors2;
+  CHECK_ITERABLE_APPROX(dvectors1, expected3);
+  dvectors1 -= dvectors2;
+  dvectors1 -= dvectors2;
+  CHECK_ITERABLE_APPROX(dvectors1, expected4);
+
+  // Test calculation of magnitude of ModalVector
+  const std::array<ModalVector, 1> d1{{ModalVector{-2.5, 3.4}}};
+  const ModalVector expected_d1{2.5, 3.4};
+  const auto magnitude_d1 = magnitude(d1);
+  CHECK_ITERABLE_APPROX(expected_d1, magnitude_d1);
+  const std::array<ModalVector, 2> d2{{ModalVector(2, 3.),
+                                             ModalVector(2, 4.)}};
+  const ModalVector expected_d2(2, 5.);
+  const auto magnitude_d2 = magnitude(d2);
+  CHECK_ITERABLE_APPROX(expected_d2, magnitude_d2);
+  const std::array<ModalVector, 3> d3{
+      {ModalVector(2, 3.), ModalVector(2, -4.),
+       ModalVector(2, 12.)}};
+  const ModalVector expected_d3(2, 13.);
+  const auto magnitude_d3 = magnitude(d3);
+  CHECK_ITERABLE_APPROX(expected_d3, magnitude_d3);
+}
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathWithDataVector",
+                  "[Unit][DataStructures]") {
+  const size_t num_pts = 23;
+  ModalVector nine(num_pts, 9.0);
+  CHECK(nine.is_owning());
+  DataVector eight(num_pts, 8.0);
+
+  // Test math with DataVector
+  ModalVector test_72(num_pts, -1.0);
+  test_72 = nine * eight;
+  check_vectors(ModalVector(num_pts, 72.0), test_72);
+  CHECK(test_72.is_owning());
+  test_72 = eight * nine;
+  check_vectors(ModalVector(num_pts, 72.0), test_72);
+  CHECK(test_72.is_owning());
+  test_72 = test_72 / eight;
+  check_vectors(ModalVector(num_pts, 9.0), test_72);
+  test_72 = nine;
+  test_72 *= eight;
+  check_vectors(ModalVector(num_pts, 72.0), test_72);
+  CHECK(test_72.is_owning());
+  test_72 /= eight;
+  check_vectors(ModalVector(num_pts, 9.0), test_72);
+  CHECK(test_72.is_owning());
+}
+
+// [[OutputRegex, Must copy into same size]]
+[[noreturn]] SPECTRE_TEST_CASE(
+    "Unit.DataStructures.ModalVector.ExpressionAssignError",
+    "[Unit][DataStructures]") {
+  ASSERTION_TEST();
+#ifdef SPECTRE_DEBUG
+  ModalVector one(10, 1.0);
+  ModalVector one_ref(one.data(), one.size());
+  ModalVector one_b(2, 1.0);
+  one_ref = (one_b * one_b);
+  ERROR("Failed to trigger ASSERT in an assertion test");
+#endif
+}