Test complex elliptic DG operator with solutions

src/PointwiseFunctions/AnalyticSolutions/Poisson/Lorentzian.cpp

@@ -4,8 +4,10 @@
 #include "PointwiseFunctions/AnalyticSolutions/Poisson/Lorentzian.hpp"
 
 #include <array>
+#include <complex>
 #include <cstddef>
 
+#include "DataStructures/ComplexDataVector.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
@@ -23,6 +25,9 @@ void LorentzianVariables<DataType, Dim>::operator()(
     const gsl::not_null<Cache*> /*cache*/,
     Tags::Field<DataType> /*meta*/) const {
   get(*field) = 1. / sqrt(1. + get(dot_product(x, x))) + constant;
+  if constexpr (std::is_same_v<DataType, ComplexDataVector>) {
+    get(*field) *= std::complex<double>{cos(complex_phase), sin(complex_phase)};
+  }
 }
 
 template <typename DataType, size_t Dim>
@@ -31,7 +36,10 @@ void LorentzianVariables<DataType, Dim>::operator()(
     const gsl::not_null<Cache*> /*cache*/,
     ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>,
                   Frame::Inertial> /*meta*/) const {
-  const DataType prefactor = -1. / cube(sqrt(1. + get(dot_product(x, x))));
+  DataType prefactor = -1. / cube(sqrt(1. + get(dot_product(x, x))));
+  if constexpr (std::is_same_v<DataType, ComplexDataVector>) {
+    prefactor *= std::complex<double>{cos(complex_phase), sin(complex_phase)};
+  }
   get<0>(*field_gradient) = prefactor * get<0>(x);
   get<1>(*field_gradient) = prefactor * get<1>(x);
   get<2>(*field_gradient) = prefactor * get<2>(x);
@@ -57,6 +65,10 @@ void LorentzianVariables<DataType, Dim>::operator()(
     const gsl::not_null<Cache*> /*cache*/,
     ::Tags::FixedSource<Tags::Field<DataType>> /*meta*/) const {
   get(*fixed_source_for_field) = 3. / pow<5>(sqrt(1. + get(dot_product(x, x))));
+  if constexpr (std::is_same_v<DataType, ComplexDataVector>) {
+    get(*fixed_source_for_field) *=
+        std::complex<double>{cos(complex_phase), sin(complex_phase)};
+  }
 }
 
 #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
@@ -65,7 +77,7 @@ void LorentzianVariables<DataType, Dim>::operator()(
 #define INSTANTIATE(_, data) \
   template class LorentzianVariables<DTYPE(data), DIM(data)>;
 
-GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector), (3))
+GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector, ComplexDataVector), (3))
 
 #undef DTYPE
 #undef DIM

src/PointwiseFunctions/AnalyticSolutions/Poisson/Lorentzian.hpp

@@ -7,7 +7,9 @@
 #include <pup.h>
 
 #include "DataStructures/CachedTempBuffer.hpp"
+#include "DataStructures/ComplexDataVector.hpp"
 #include "DataStructures/DataBox/Prefixes.hpp"
+#include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Elliptic/Systems/Poisson/Tags.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
@@ -29,8 +31,9 @@ struct LorentzianVariables {
       ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
       ::Tags::FixedSource<Tags::Field<DataType>>>;
 
-  const tnsr::I<DataType, Dim>& x;
+  const tnsr::I<DataVector, Dim>& x;
   const double constant;
+  const double complex_phase;
 
   void operator()(gsl::not_null<Scalar<DataType>*> field,
                   gsl::not_null<Cache*> cache,
@@ -59,9 +62,13 @@ struct LorentzianVariables {
  * \f$r^2=x^2+y^2+z^2\f$. The corresponding source is
  * \f$f(\boldsymbol{x})=3\left(1+r^2\right)^{-\frac{5}{2}}\f$.
  *
+ * If `DataType` is `ComplexDataVector`, the solution is multiplied by
+ * `exp(i * complex_phase)` to rotate it in the complex plane. This allows to
+ * use this solution for the complex Poisson equation.
+ *
  * \note Corresponding 1D and 2D solutions are not implemented yet.
  */
-template <size_t Dim>
+template <size_t Dim, typename DataType = DataVector>
 class Lorentzian : public elliptic::analytic_data::AnalyticSolution {
   static_assert(
       Dim == 3,
@@ -73,7 +80,17 @@ class Lorentzian : public elliptic::analytic_data::AnalyticSolution {
     static constexpr Options::String help{"Constant added to the solution."};
   };
 
-  using options = tmpl::list<PlusConstant>;
+  struct ComplexPhase {
+    using type = double;
+    static constexpr Options::String help{
+        "Phase 'phi' of a complex exponential 'exp(i phi)' that rotates the "
+        "solution in the complex plane."};
+  };
+
+  using options = tmpl::flatten<tmpl::list<
+      PlusConstant,
+      tmpl::conditional_t<std::is_same_v<DataType, ComplexDataVector>,
+                          ComplexPhase, tmpl::list<>>>>;
   static constexpr Options::String help{
       "A Lorentzian solution to the Poisson equation."};
 
@@ -84,9 +101,14 @@ class Lorentzian : public elliptic::analytic_data::AnalyticSolution {
   Lorentzian& operator=(Lorentzian&&) = default;
   ~Lorentzian() override = default;
 
-  explicit Lorentzian(const double constant) : constant_(constant) {}
+  explicit Lorentzian(const double constant, const double complex_phase = 0.)
+      : constant_(constant), complex_phase_(complex_phase) {
+    ASSERT((std::is_same_v<DataType, ComplexDataVector> or complex_phase == 0.),
+           "The complex phase is only supported for ComplexDataVector.");
+  }
 
   double constant() const { return constant_; }
+  double complex_phase() const { return complex_phase_; }
 
   std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
       const override {
@@ -100,38 +122,42 @@ class Lorentzian : public elliptic::analytic_data::AnalyticSolution {
   WRAPPED_PUPable_decl_template(Lorentzian);  // NOLINT
   /// \endcond
 
-  template <typename DataType, typename... RequestedTags>
+  template <typename... RequestedTags>
   tuples::TaggedTuple<RequestedTags...> variables(
-      const tnsr::I<DataType, Dim>& x,
+      const tnsr::I<DataVector, Dim>& x,
       tmpl::list<RequestedTags...> /*meta*/) const {
     using VarsComputer = detail::LorentzianVariables<DataType, Dim>;
     typename VarsComputer::Cache cache{get_size(*x.begin())};
-    const VarsComputer computer{x, constant_};
+    const VarsComputer computer{x, constant_, complex_phase_};
     return {cache.get_var(computer, RequestedTags{})...};
   }
 
   void pup(PUP::er& p) override {
     elliptic::analytic_data::AnalyticSolution::pup(p);
     p | constant_;
+    p | complex_phase_;
   }
 
  private:
   double constant_ = std::numeric_limits<double>::signaling_NaN();
+  double complex_phase_ = std::numeric_limits<double>::signaling_NaN();
 };
 
 /// \cond
-template <size_t Dim>
-PUP::able::PUP_ID Lorentzian<Dim>::my_PUP_ID = 0;  // NOLINT
+template <size_t Dim, typename DataType>
+PUP::able::PUP_ID Lorentzian<Dim, DataType>::my_PUP_ID = 0;  // NOLINT
 /// \endcond
 
-template <size_t Dim>
-bool operator==(const Lorentzian<Dim>& lhs,
-                const Lorentzian<Dim>& rhs) {
-  return lhs.constant() == rhs.constant();
+template <size_t Dim, typename DataType>
+bool operator==(const Lorentzian<Dim, DataType>& lhs,
+                const Lorentzian<Dim, DataType>& rhs) {
+  return lhs.constant() == rhs.constant() and
+         lhs.complex_phase() == rhs.complex_phase();
 }
 
-template <size_t Dim>
-bool operator!=(const Lorentzian<Dim>& lhs, const Lorentzian<Dim>& rhs) {
+template <size_t Dim, typename DataType>
+bool operator!=(const Lorentzian<Dim, DataType>& lhs,
+                const Lorentzian<Dim, DataType>& rhs) {
   return not(lhs == rhs);
 }
 

src/PointwiseFunctions/AnalyticSolutions/Poisson/ProductOfSinusoids.cpp

@@ -4,8 +4,10 @@
 #include "PointwiseFunctions/AnalyticSolutions/Poisson/ProductOfSinusoids.hpp"
 
 #include <cmath>
+#include <complex>
 #include <cstddef>
 
+#include "DataStructures/ComplexDataVector.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
@@ -26,6 +28,9 @@ void ProductOfSinusoidsVariables<DataType, Dim>::operator()(
   for (size_t d = 0; d < Dim; d++) {
     get(*field) *= sin(gsl::at(wave_numbers, d) * x.get(d));
   }
+  if constexpr (std::is_same_v<DataType, ComplexDataVector>) {
+    get(*field) *= std::complex<double>{cos(complex_phase), sin(complex_phase)};
+  }
 }
 
 template <typename DataType, size_t Dim>
@@ -43,6 +48,10 @@ void ProductOfSinusoidsVariables<DataType, Dim>::operator()(
             sin(gsl::at(wave_numbers, other_d) * x.get(other_d));
       }
     }
+    if constexpr (std::is_same_v<DataType, ComplexDataVector>) {
+      field_gradient->get(d) *=
+          std::complex<double>{cos(complex_phase), sin(complex_phase)};
+    }
   }
 }
 
@@ -75,7 +84,7 @@ void ProductOfSinusoidsVariables<DataType, Dim>::operator()(
 #define INSTANTIATE(_, data) \
   template class ProductOfSinusoidsVariables<DTYPE(data), DIM(data)>;
 
-GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector), (1, 2, 3))
+GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector, ComplexDataVector), (1, 2, 3))
 
 #undef DTYPE
 #undef DIM

src/PointwiseFunctions/AnalyticSolutions/Poisson/ProductOfSinusoids.hpp

@@ -9,7 +9,9 @@
 #include <pup.h>
 
 #include "DataStructures/CachedTempBuffer.hpp"
+#include "DataStructures/ComplexDataVector.hpp"
 #include "DataStructures/DataBox/Prefixes.hpp"
+#include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Elliptic/Systems/Poisson/Tags.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
@@ -31,8 +33,9 @@ struct ProductOfSinusoidsVariables {
       ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
       ::Tags::FixedSource<Tags::Field<DataType>>>;
 
-  const tnsr::I<DataType, Dim>& x;
+  const tnsr::I<DataVector, Dim>& x;
   const std::array<double, Dim>& wave_numbers;
+  const double complex_phase;
 
   void operator()(gsl::not_null<Scalar<DataType>*> field,
                   gsl::not_null<Cache*> cache,
@@ -56,16 +59,30 @@ struct ProductOfSinusoidsVariables {
  *
  * \details Solves the Poisson equation \f$-\Delta u(x)=f(x)\f$ for a source
  * \f$f(x)=\boldsymbol{k}^2\prod_i \sin(k_i x_i)\f$.
+ *
+ * If `DataType` is `ComplexDataVector`, the solution is multiplied by
+ * `exp(i * complex_phase)` to rotate it in the complex plane. This allows to
+ * use this solution for the complex Poisson equation.
  */
-template <size_t Dim>
+template <size_t Dim, typename DataType = DataVector>
 class ProductOfSinusoids : public elliptic::analytic_data::AnalyticSolution {
  public:
   struct WaveNumbers {
     using type = std::array<double, Dim>;
     static constexpr Options::String help{"The wave numbers of the sinusoids"};
   };
 
-  using options = tmpl::list<WaveNumbers>;
+  struct ComplexPhase {
+    using type = double;
+    static constexpr Options::String help{
+        "Phase 'phi' of a complex exponential 'exp(i phi)' that rotates the "
+        "solution in the complex plane."};
+  };
+
+  using options = tmpl::flatten<tmpl::list<
+      WaveNumbers,
+      tmpl::conditional_t<std::is_same_v<DataType, ComplexDataVector>,
+                          ComplexPhase, tmpl::list<>>>>;
   static constexpr Options::String help{
       "A product of sinusoids that are taken of a wave number times the "
       "coordinate in each dimension."};
@@ -88,46 +105,54 @@ class ProductOfSinusoids : public elliptic::analytic_data::AnalyticSolution {
   WRAPPED_PUPable_decl_template(ProductOfSinusoids);  // NOLINT
   /// \endcond
 
-  explicit ProductOfSinusoids(const std::array<double, Dim>& wave_numbers)
-      : wave_numbers_(wave_numbers) {}
+  explicit ProductOfSinusoids(const std::array<double, Dim>& wave_numbers,
+                              const double complex_phase = 0.)
+      : wave_numbers_(wave_numbers), complex_phase_(complex_phase) {
+    ASSERT((std::is_same_v<DataType, ComplexDataVector> or complex_phase == 0.),
+           "The complex phase is only supported for ComplexDataVector.");
+  }
 
-  template <typename DataType, typename... RequestedTags>
+  template <typename... RequestedTags>
   tuples::TaggedTuple<RequestedTags...> variables(
-      const tnsr::I<DataType, Dim>& x,
+      const tnsr::I<DataVector, Dim>& x,
       tmpl::list<RequestedTags...> /*meta*/) const {
     using VarsComputer = detail::ProductOfSinusoidsVariables<DataType, Dim>;
     typename VarsComputer::Cache cache{get_size(*x.begin())};
-    const VarsComputer computer{x, wave_numbers_};
+    const VarsComputer computer{x, wave_numbers_, complex_phase_};
     return {cache.get_var(computer, RequestedTags{})...};
   }
 
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p) override {
     elliptic::analytic_data::AnalyticSolution::pup(p);
     p | wave_numbers_;
+    p | complex_phase_;
   }
 
   const std::array<double, Dim>& wave_numbers() const { return wave_numbers_; }
+  double complex_phase() const { return complex_phase_; }
 
  private:
   std::array<double, Dim> wave_numbers_{
       {std::numeric_limits<double>::signaling_NaN()}};
+  double complex_phase_ = std::numeric_limits<double>::signaling_NaN();
 };
 
 /// \cond
-template <size_t Dim>
-PUP::able::PUP_ID ProductOfSinusoids<Dim>::my_PUP_ID = 0;  // NOLINT
+template <size_t Dim, typename DataType>
+PUP::able::PUP_ID ProductOfSinusoids<Dim, DataType>::my_PUP_ID = 0;  // NOLINT
 /// \endcond
 
-template <size_t Dim>
-bool operator==(const ProductOfSinusoids<Dim>& lhs,
-                const ProductOfSinusoids<Dim>& rhs) {
-  return lhs.wave_numbers() == rhs.wave_numbers();
+template <size_t Dim, typename DataType>
+bool operator==(const ProductOfSinusoids<Dim, DataType>& lhs,
+                const ProductOfSinusoids<Dim, DataType>& rhs) {
+  return lhs.wave_numbers() == rhs.wave_numbers() and
+         lhs.complex_phase() == rhs.complex_phase();
 }
 
-template <size_t Dim>
-bool operator!=(const ProductOfSinusoids<Dim>& lhs,
-                const ProductOfSinusoids<Dim>& rhs) {
+template <size_t Dim, typename DataType>
+bool operator!=(const ProductOfSinusoids<Dim, DataType>& lhs,
+                const ProductOfSinusoids<Dim, DataType>& rhs) {
   return not(lhs == rhs);
 }