diff --git a/src/DataStructures/CMakeLists.txt b/src/DataStructures/CMakeLists.txt
index 744a928e2e2e..9874fbff7e4a 100644
--- a/src/DataStructures/CMakeLists.txt
+++ b/src/DataStructures/CMakeLists.txt
@@ -8,6 +8,7 @@ set(LIBRARY_SOURCES
     Index.cpp
     IndexIterator.cpp
     LeviCivitaIterator.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..cd99f18f76ff
--- /dev/null
+++ b/src/DataStructures/ModalVector.cpp
@@ -0,0 +1,115 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "DataStructures/ModalVector.hpp"
+
+#include <algorithm>
+#include <pup.h>
+
+#include "Utilities/StdHelpers.hpp"
+
+ModalVector::ModalVector(const size_t size, const double value)
+    : owned_data_(size, value) {
+  reset_pointer_vector();
+}
+
+ModalVector::ModalVector(double* start, size_t size) noexcept
+    : BaseType(start, size), owned_data_(0), owning_(false) {}
+
+template <class T, Requires<cpp17::is_same_v<T, double>>>
+ModalVector::ModalVector(std::initializer_list<T> list)
+    : owned_data_(std::move(list)) {
+  reset_pointer_vector();
+}
+
+/// \cond HIDDEN_SYMBOLS
+// clang-tidy: calling a base constructor other than the copy constructor.
+//             We reset the base class in reset_pointer_vector after calling its
+//             default constructor
+ModalVector::ModalVector(const ModalVector& rhs) : BaseType{} {  // NOLINT
+  if (rhs.is_owning()) {
+    owned_data_ = rhs.owned_data_;
+  } else {
+    owned_data_.assign(rhs.begin(), rhs.end());
+  }
+  reset_pointer_vector();
+}
+
+ModalVector& ModalVector::operator=(const ModalVector& rhs) {
+  if (this != &rhs) {
+    if (owning_) {
+      if (rhs.is_owning()) {
+        owned_data_ = rhs.owned_data_;
+      } else {
+        owned_data_.assign(rhs.begin(), rhs.end());
+      }
+      reset_pointer_vector();
+    } 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 {
+  owned_data_ = std::move(rhs.owned_data_);
+  ~*this = ~rhs;  // PointerVector is trivially copyable
+  owning_ = rhs.owning_;
+
+  rhs.owning_ = true;
+  rhs.reset();
+}
+
+ModalVector& ModalVector::operator=(ModalVector&& rhs) noexcept {
+  if (this != &rhs) {
+    if (owning_) {
+      owned_data_ = std::move(rhs.owned_data_);
+      ~*this = ~rhs;  // PointerVector is trivially copyable
+      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.reset();
+  }
+  return *this;
+}
+/// \endcond
+
+void ModalVector::pup(PUP::er& p) noexcept {  // NOLINT
+  auto my_size = size();
+  p | my_size;
+  if (my_size > 0) {
+    if (p.isUnpacking()) {
+      owning_ = true;
+      owned_data_.resize(my_size);
+      reset_pointer_vector();
+    }
+    PUParray(p, data(), size());
+  }
+}
+
+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) noexcept {
+  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) noexcept {
+  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..5e3d80fc77e9
--- /dev/null
+++ b/src/DataStructures/ModalVector.hpp
@@ -0,0 +1,469 @@
+// 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 "ErrorHandling/Assert.hpp"
+#include "Utilities/ForceInline.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/MakeWithValue.hpp"
+#include "Utilities/PointerVector.hpp" // IWYU pragma: keep
+#include "Utilities/Requires.hpp"
+#include "Utilities/TMPL.hpp"          // for list
+
+/// \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
+
+// IWYU doesn't like that we want PointerVector.hpp to expose Blaze and also
+// have ModalVector.hpp to expose PointerVector.hpp without including Blaze
+// directly in ModalVector.hpp
+//
+// IWYU pragma: no_include <blaze/math/dense/DenseVector.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecAddExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecDivExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecMultExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecSubExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecMapExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecScalarDivExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecScalarMultExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DenseVector.h>
+// IWYU pragma: no_include <blaze/math/expressions/Vector.h>
+// IWYU pragma: no_include <blaze/math/typetraits/IsVector.h>
+// IWYU pragma: no_include <blaze/math/expressions/Forward.h>
+// IWYU pragma: no_include <blaze/math/AlignmentFlag.h>
+// IWYU pragma: no_include <blaze/math/PaddingFlag.h>
+// IWYU pragma: no_include <blaze/math/traits/AddTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/DivTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/MultTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/SubTrait.h>
+// IWYU pragma: no_include <blaze/system/TransposeFlag.h>
+// IWYU pragma: no_include <blaze/math/traits/BinaryMapTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/UnaryMapTrait.h>
+// IWYU pragma: no_include <blaze/math/typetraits/TransposeFlag.h>
+// IWYU pragma: no_include "DataStructures/DataVector.hpp"
+
+// IWYU pragma: no_forward_declare blaze::DenseVector
+// IWYU pragma: no_forward_declare blaze::UnaryMapTrait
+// IWYU pragma: no_forward_declare blaze::BinaryMapTrait
+// IWYU pragma: no_forward_declare blaze::IsVector
+// IWYU pragma: no_forward_declare blaze::TransposeFlag
+
+/*!
+ * \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/fabs/magnitude
+ * - max
+ * - min
+ *
+ * In order to allow filtering, multiplication (*, *=) and division (/, /=)
+ * operations with a DenseVectors (holding filters) is supported.
+ *
+ */
+class ModalVector
+    : public PointerVector<double, blaze::unaligned, blaze::unpadded,
+                           blaze::defaultTransposeFlag, 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;
+  using BaseType = PointerVector<value_type, blaze::unaligned, blaze::unpadded,
+                                 blaze::defaultTransposeFlag, ModalVector>;
+  static constexpr bool transpose_flag = blaze::defaultTransposeFlag;
+
+  using BaseType::ElementType;
+  using TransposeType = ModalVector;
+  using CompositeType = const ModalVector&;
+
+  /// Iterators etc obtained from base class PointerVector
+  using BaseType::operator[];
+  using BaseType::begin;
+  using BaseType::cbegin;
+  using BaseType::cend;
+  using BaseType::data;
+  using BaseType::end;
+  using BaseType::size;
+
+  /// Do we need an upcast from ModalVector to PointerVector?
+  // @{
+  // Upcast to `BaseType`
+  const BaseType& operator~() const noexcept {
+    return static_cast<const BaseType&>(*this);
+  }
+  BaseType& operator~() noexcept { return static_cast<BaseType&>(*this); }
+  // @}
+
+  /// 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) noexcept;
+
+  /// Create from an initializer list of doubles (must have decimal points)
+  template <class T, Requires<cpp17::is_same_v<T, double>> = nullptr>
+  ModalVector(std::initializer_list<T> list);
+
+  /// Empty ModalVector
+  ModalVector() noexcept = default;
+  /// \cond HIDDEN_SYMBOLS
+  ~ModalVector() = default;
+
+  // Initialize ModalVector with expressions involving itself
+  ModalVector(const ModalVector& rhs);
+  ModalVector(ModalVector&& rhs) noexcept;
+  ModalVector& operator=(const ModalVector& rhs);
+  ModalVector& operator=(ModalVector&& rhs) noexcept;
+
+  /// Catch DenseVector expressions into ModalVectors, only if their return
+  /// types are defined as ModalVectors
+  //
+  // 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::DenseVector<VT, VF>& expression) noexcept;  // NOLINT
+
+  template <typename VT, bool VF>
+  ModalVector& operator=(const blaze::DenseVector<VT, VF>& expression) noexcept;
+  /// \endcond
+
+  ModalVector& operator=(const double& rhs) noexcept {
+    ~*this = rhs;
+    return *this;
+  }
+
+  // Using expression macro, defined in PointerVector.hpp, these lines overload
+  // assignment operators { +=, -=, *=, /= } for LHS = ModalVector and RHS =
+  // (i) ModalVector, (ii) blaze::DenseVector (to hold spectral filters), and
+  // (iii) double (i.e. the ElementType)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(+=, ModalVector)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(-=, ModalVector)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(*=, ModalVector)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(/=, ModalVector)
+
+  // @{
+  /// Set 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 {
+    owned_data_ = decltype(owned_data_){};
+    (~*this).reset(start, size);
+    owning_ = false;
+  }
+  // @}
+
+  /// Returns true if the class owns the data
+  bool is_owning() const noexcept { return owning_; }
+
+  /// Serialization for Charm++
+  // clang-tidy: google-runtime-references
+  void pup(PUP::er& p) noexcept;  // NOLINT
+
+ private:
+  SPECTRE_ALWAYS_INLINE void reset_pointer_vector() noexcept {
+    reset(owned_data_.data(), owned_data_.size());
+  }
+
+  /// \cond HIDDEN_SYMBOLS
+  std::vector<value_type, allocator_type> owned_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) noexcept;
+
+/// Inequivalence operator for ModalVector
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs) noexcept;
+
+/// \cond
+// Used for comparing ModalVector to an expression
+template <typename VT, bool VF>
+bool operator==(const ModalVector& lhs,
+                const blaze::DenseVector<VT, VF>& rhs) noexcept {
+  return lhs == ModalVector(rhs);
+}
+
+template <typename VT, bool VF>
+bool operator!=(const ModalVector& lhs,
+                const blaze::DenseVector<VT, VF>& rhs) noexcept {
+  return not(lhs == rhs);
+}
+
+template <typename VT, bool VF>
+bool operator==(const blaze::DenseVector<VT, VF>& lhs,
+                const ModalVector& rhs) noexcept {
+  return ModalVector(lhs) == rhs;
+}
+
+template <typename VT, bool VF>
+bool operator!=(const blaze::DenseVector<VT, VF>& lhs,
+                const ModalVector& rhs) noexcept {
+  return not(lhs == rhs);
+}
+/// \endcond
+
+// Specialize the Blaze type traits to correctly handle ModalVector
+namespace blaze {
+template <>
+struct IsVector<ModalVector> : std::true_type {};
+
+template <>
+struct TransposeFlag<ModalVector> : BoolConstant<
+                ModalVector::transpose_flag> {};
+
+template <>
+struct AddTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct AddTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
+
+template <>
+struct AddTrait<double, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct SubTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct SubTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
+
+template <>
+struct SubTrait<double, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct MultTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct MultTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
+
+template <>
+struct MultTrait<double, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct DivTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct DivTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
+
+// Forbid math operations in this specialization of UnaryMap traits for
+// ModalVector that are unlikely to be used on spectral coefficients
+template <typename Operator>
+struct UnaryMapTrait<ModalVector, Operator> {
+  static_assert(not tmpl::list_contains_v<tmpl::list<//blaze::Max, blaze::Min,
+                //blaze::Sqrt,
+                blaze::Cbrt,
+                blaze::InvSqrt, blaze::InvCbrt,
+                blaze::Acos, blaze::Acosh, blaze::Cos, blaze::Cosh,
+                blaze::Asin, blaze::Asinh, blaze::Sin, blaze::Sinh,
+                blaze::Atan, blaze::Atan2, blaze::Atanh,
+                blaze::Tan, blaze::Tanh, blaze::Hypot,
+                blaze::Exp, blaze::Exp2, blaze::Exp10,
+                blaze::Log, blaze::Log2, blaze::Log10,
+                blaze::Erf, blaze::Erfc, blaze::StepFunction
+                >, Operator>,
+                "This operation is not permitted on a ModalVector."
+                "Only unary operation permitted are: abs, fabs, sqrt.");
+  using Type = ModalVector;
+};
+
+// Forbid math operations in this specialization of BinaryMap traits for
+// ModalVector that are unlikely to be used on spectral coefficients
+template <typename Operator>
+struct BinaryMapTrait<ModalVector, ModalVector, Operator> {
+  static_assert(not tmpl::list_contains_v<tmpl::list<blaze::Max, blaze::Min
+                >, Operator>,
+                "This operation is not permitted on a ModalVector."
+                "Only unary operation are permitted: abs, fabs, sqrt.");
+  using Type = ModalVector;
+};
+
+}  // namespace blaze
+
+SPECTRE_ALWAYS_INLINE decltype(auto) fabs(const ModalVector& t) noexcept {
+  return abs(~t);
+}
+
+SPECTRE_ALWAYS_INLINE decltype(auto) abs(const ModalVector& t) noexcept {
+  return abs(~t);
+}
+
+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::DenseVector<VT, VF>& expression) noexcept
+    : owned_data_((~expression).size()) {
+  static_assert(cpp17::is_same_v<typename VT::ResultType, ModalVector>,
+                "You are attempting to assign the result of an expression that "
+                "is not a ModalVector to a ModalVector.");
+  reset_pointer_vector();
+  ~*this = expression;
+}
+
+template <typename VT, bool VF>
+ModalVector& ModalVector::operator=(
+    const blaze::DenseVector<VT, VF>& expression) noexcept {
+  static_assert(cpp17::is_same_v<typename VT::ResultType, ModalVector>,
+                "You are attempting to assign the result of an expression that "
+                "is not a ModalVector to a ModalVector.");
+  if (owning_ and (~expression).size() != size()) {
+    owned_data_.resize((~expression).size());
+    reset_pointer_vector();
+  } else if (not owning_) {
+    ASSERT((~expression).size() == size(), "Must copy into same size");
+  }
+  ~*this = expression;
+  return *this;
+}
+/// \endcond
+
+namespace MakeWithValueImpls {
+/// \brief Returns a ModalVector the same size as `input`, with each element
+/// equal to `value`.
+template <>
+SPECTRE_ALWAYS_INLINE ModalVector
+MakeWithValueImpl<ModalVector, ModalVector>::apply(const ModalVector& input,
+                                                 const double value) {
+  return ModalVector(input.size(), value);
+}
+}  // namespace MakeWithValueImpls
diff --git a/tests/Unit/DataStructures/CMakeLists.txt b/tests/Unit/DataStructures/CMakeLists.txt
index 85c05e6cbf84..020f1bbbd6f3 100644
--- a/tests/Unit/DataStructures/CMakeLists.txt
+++ b/tests/Unit/DataStructures/CMakeLists.txt
@@ -12,6 +12,7 @@ set(LIBRARY_SOURCES
   Test_Index.cpp
   Test_IndexIterator.cpp
   Test_LeviCivitaIterator.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..288e43b0a8be
--- /dev/null
+++ b/tests/Unit/DataStructures/Test_ModalVector.cpp
@@ -0,0 +1,533 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "tests/Unit/TestingFramework.hpp"
+
+#include <algorithm>                         // for move
+#include <array>                             // for array, operator==
+#include <cmath>                             // for abs
+#include <cstddef>                           // for size_t
+#include <functional>                        // for std::reference_wrapper
+#include <numeric>                           // for iota
+
+#include "DataStructures/ModalVector.hpp"
+#include "ErrorHandling/Error.hpp"
+#include "Utilities/DereferenceWrapper.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/Requires.hpp"
+#include "Utilities/StdArrayHelpers.hpp"
+#include "Utilities/StdHelpers.hpp"  // IWYU pragma: keep
+#include "tests/Unit/TestHelpers.hpp"
+
+// IWYU wants to include DataVector.hpp to access the function blaze::smpAssign
+//
+// IWYU pragma: no_include "DataStructures/DataVector.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 = 100;
+  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 = 11;
+  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};
+  CHECK(data.is_owning());
+  ModalVector t;
+  t.set_data_ref(&data);
+  CHECK(not t.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 that data_2 now contain data_3's values
+    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.is_owning());
+    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 data2{1.43, 2.83, 5.94};
+  ModalVector data_ref;
+  data_ref.set_data_ref(&data);
+  CHECK(data_ref[2] == 3.94);
+  data_ref = data2;
+  CHECK(data_ref[2] == 5.94);
+  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 data2{1.43, 2.83, 3.94};
+  ModalVector data_ref;
+  data_ref.set_data_ref(&data);
+  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
+
+/// Tests of ModalVector math
+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));
+  }
+}
+
+namespace {
+enum class UseRefWrap { None, Cref, Ref };
+
+template <UseRefWrap Wrap, class T,
+          Requires<Wrap == UseRefWrap::Cref> = nullptr>
+decltype(auto) wrap(const T& t) noexcept {
+  return std::cref(t);
+}
+template <UseRefWrap Wrap, class T,
+          Requires<Wrap == UseRefWrap::Ref> = nullptr>
+decltype(auto) wrap(T& t) noexcept {
+  return std::ref(t);
+}
+template <UseRefWrap Wrap, class T,
+          Requires<Wrap == UseRefWrap::None> = nullptr>
+decltype(auto) wrap(const T& t) noexcept {
+  return t;
+}
+
+// Wrap is used to wrap values in a std::reference_wrapper using std::cref and
+// std::ref, or to not wrap at all. This is done to verify that all math
+// operations work transparently with a `std::reference_wrapper` too.
+template <UseRefWrap WrapLeftOp, UseRefWrap WrapRightOp>
+void test_ModalVector_math() noexcept {
+  constexpr size_t num_pts = 19;
+  ModalVector val{1., 2., 3., -4., 8., 12., -14.};
+  ModalVector one(num_pts, 1.0);
+  ModalVector eight(num_pts, 8.0);
+  ModalVector nine(num_pts, 9.0);
+
+  // Test unary minus
+  check_vectors(-wrap<WrapLeftOp>(nine), ModalVector(num_pts, -9.0));
+
+  check_vectors(wrap<WrapLeftOp>(nine) + 2.0, ModalVector(num_pts, 11.0));
+  check_vectors(2.0 + wrap<WrapLeftOp>(nine), ModalVector(num_pts, 11.0));
+  check_vectors(wrap<WrapLeftOp>(nine) - 2.0, ModalVector(num_pts, 7.0));
+  check_vectors(2.0 - wrap<WrapLeftOp>(nine), ModalVector(num_pts, -7.0));
+  check_vectors(wrap<WrapLeftOp>(nine) + wrap<WrapRightOp>(nine),
+                ModalVector(num_pts, 18.0));
+  check_vectors(wrap<WrapLeftOp>(nine) +
+                    (wrap<WrapLeftOp>(one) * wrap<WrapRightOp>(nine)),
+                ModalVector(num_pts, 18.0));
+  check_vectors((wrap<WrapLeftOp>(one) * wrap<WrapLeftOp>(nine)) +
+                    wrap<WrapRightOp>(nine),
+                ModalVector(num_pts, 18.0));
+  check_vectors(wrap<WrapLeftOp>(nine) - ModalVector(num_pts, 8.0),
+                ModalVector(num_pts, 1.0));
+  check_vectors(wrap<WrapLeftOp>(nine) -
+                    (wrap<WrapRightOp>(one) * wrap<WrapLeftOp>(nine)),
+                ModalVector(num_pts, 0.0));
+  check_vectors((wrap<WrapLeftOp>(one) * wrap<WrapLeftOp>(nine)) -
+                    wrap<WrapRightOp>(nine),
+                ModalVector(num_pts, 0.0));
+
+  check_vectors(ModalVector(num_pts, -1.0 / 9.0),
+                -one / wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, -8.0 / 9.0),
+                -(wrap<WrapLeftOp>(nine) - wrap<WrapRightOp>(one)) /
+                    wrap<WrapLeftOp>(nine));
+  check_vectors(ModalVector(num_pts, 18.0),
+                (wrap<WrapLeftOp>(one) / 0.5) * wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, 1.0), 9.0 / wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, 1.0),
+                (wrap<WrapLeftOp>(one) * 9.0) / wrap<WrapRightOp>(nine));
+
+  check_vectors(ModalVector(num_pts, 81.0),
+                wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, 81.0),
+                wrap<WrapLeftOp>(nine) *
+                    (wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(one)));
+  check_vectors(ModalVector(num_pts, 81.0),
+                (wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(nine)) *
+                    wrap<WrapLeftOp>(one));
+  check_vectors(ModalVector(num_pts, 81.0), 9.0 * wrap<WrapLeftOp>(nine));
+  check_vectors(ModalVector(num_pts, 81.0), wrap<WrapLeftOp>(nine) * 9.0);
+  check_vectors(ModalVector(num_pts, 1.0), wrap<WrapLeftOp>(nine) / 9.0);
+  check_vectors(ModalVector(num_pts, 1.0),
+                wrap<WrapLeftOp>(nine) / (wrap<WrapRightOp>(one) * 9.0));
+  check_vectors(ModalVector(num_pts, 9.0),
+                (wrap<WrapLeftOp>(one) * wrap<WrapLeftOp>(nine)) /
+                    wrap<WrapRightOp>(one));
+
+  CHECK(-14 == min(wrap<WrapLeftOp>(val)));
+  CHECK(12 == max(wrap<WrapLeftOp>(val)));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.},
+                abs(wrap<WrapLeftOp>(val)));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.},
+                fabs(wrap<WrapLeftOp>(val)));
+
+  check_vectors(sqrt(wrap<WrapLeftOp>(nine)), ModalVector(num_pts, 3.0));
+
+  ModalVector dummy(wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(nine) * 1.0);
+  check_vectors(ModalVector(num_pts, 81.0), dummy);
+
+  // Test assignment
+  ModalVector test_81(num_pts, -1.0);
+  // cannot assign to a reference_wrapper using an expr :(
+  test_81 = wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(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 = wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(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 = wrap<WrapLeftOp>(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 = wrap<WrapLeftOp>(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 += wrap<WrapLeftOp>(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 += (wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(nine));
+  check_vectors(ModalVector(num_pts, 100.0), test_assignment);
+
+  test_assignment = 7.0;
+  test_assignment -= wrap<WrapLeftOp>(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 -= wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(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 *= wrap<WrapLeftOp>(nine);
+  check_vectors(ModalVector(num_pts, 54.0), test_assignment);
+  test_assignment = 1.0;
+  test_assignment *= wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(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 /= ModalVector(num_pts, 0.5);
+  check_vectors(ModalVector(num_pts, 2.0), test_assignment);
+  test_assignment /= (ModalVector(num_pts, 2.0) * ModalVector(num_pts, 3.0));
+  check_vectors(ModalVector(num_pts, 1.0 / 3.0), test_assignment);
+
+  // Test assignment where the RHS is an expression that contains the LHS
+  ModalVector x(num_pts, 4.);
+  x += sqrt(wrap<WrapLeftOp>(x));
+  check_vectors(ModalVector(num_pts, 6.0), x);
+  x -= sqrt(wrap<WrapLeftOp>(x) - 2.0);
+  check_vectors(ModalVector(num_pts, 4.0), x);
+  x = sqrt(wrap<WrapLeftOp>(x));
+  check_vectors(ModalVector(num_pts, 2.0), x);
+  x *= wrap<WrapLeftOp>(x);
+  check_vectors(ModalVector(num_pts, 4.0), x);
+  x /= wrap<WrapLeftOp>(x);
+  check_vectors(ModalVector(num_pts, 1.0), x);
+
+  // Test composition of constant expressions with ModalVector math member
+  // functions
+  x = ModalVector(num_pts, 2.);
+  check_vectors(ModalVector(num_pts, 1.4142135623730951),
+                sqrt(abs(wrap<WrapLeftOp>(x))));
+  check_vectors(ModalVector(num_pts, 2.),
+                sqrt(wrap<WrapLeftOp>(x) * wrap<WrapRightOp>(x)));
+}
+
+void test_ModalVector_array_math() noexcept {
+  // 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);
+}
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
+                  "[Unit][DataStructures]") {
+  test_ModalVector_math<UseRefWrap::Cref, UseRefWrap::Cref>();
+  test_ModalVector_math<UseRefWrap::Cref, UseRefWrap::Ref>();
+  test_ModalVector_math<UseRefWrap::Cref, UseRefWrap::None>();
+  test_ModalVector_math<UseRefWrap::Ref, UseRefWrap::Cref>();
+  test_ModalVector_math<UseRefWrap::Ref, UseRefWrap::Ref>();
+  test_ModalVector_math<UseRefWrap::Ref, UseRefWrap::None>();
+  test_ModalVector_math<UseRefWrap::None, UseRefWrap::Cref>();
+  test_ModalVector_math<UseRefWrap::None, UseRefWrap::Ref>();
+  test_ModalVector_math<UseRefWrap::None, UseRefWrap::None>();
+
+  test_ModalVector_array_math();
+}
+
+// [[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
+}
+