diff --git a/src/Elliptic/Systems/Poisson/Equations.cpp b/src/Elliptic/Systems/Poisson/Equations.cpp index f73480ee87cb..11b5912dcce6 100644 --- a/src/Elliptic/Systems/Poisson/Equations.cpp +++ b/src/Elliptic/Systems/Poisson/Equations.cpp @@ -5,6 +5,7 @@ #include +#include "DataStructures/ComplexDataVector.hpp" #include "DataStructures/DataBox/PrefixHelpers.hpp" #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/EagerMath/DotProduct.hpp" @@ -15,111 +16,114 @@ namespace Poisson { -template +template void flat_cartesian_fluxes( - const gsl::not_null*> flux_for_field, - const tnsr::i& field_gradient) { + const gsl::not_null*> flux_for_field, + const tnsr::i& field_gradient) { for (size_t d = 0; d < Dim; d++) { flux_for_field->get(d) = field_gradient.get(d); } } -template -void curved_fluxes( - const gsl::not_null*> flux_for_field, - const tnsr::II& inv_spatial_metric, - const tnsr::i& field_gradient) { +template +void curved_fluxes(const gsl::not_null*> flux_for_field, + const tnsr::II& inv_spatial_metric, + const tnsr::i& field_gradient) { raise_or_lower_index(flux_for_field, field_gradient, inv_spatial_metric); } -template -void fluxes_on_face(gsl::not_null*> flux_for_field, +template +void fluxes_on_face(gsl::not_null*> flux_for_field, const tnsr::I& face_normal_vector, - const Scalar& field) { - *flux_for_field = face_normal_vector; + const Scalar& field) { + std::copy(face_normal_vector.begin(), face_normal_vector.end(), + flux_for_field->begin()); for (size_t d = 0; d < Dim; d++) { flux_for_field->get(d) *= get(field); } } -template -void add_curved_sources( - const gsl::not_null*> source_for_field, - const tnsr::i& christoffel_contracted, - const tnsr::I& flux_for_field) { +template +void add_curved_sources(const gsl::not_null*> source_for_field, + const tnsr::i& christoffel_contracted, + const tnsr::I& flux_for_field) { get(*source_for_field) -= get(dot_product(christoffel_contracted, flux_for_field)); } -template -void Fluxes::apply( - const gsl::not_null*> flux_for_field, - const Scalar& /*field*/, - const tnsr::i& field_gradient) { +template +void Fluxes::apply( + const gsl::not_null*> flux_for_field, + const Scalar& /*field*/, + const tnsr::i& field_gradient) { flat_cartesian_fluxes(flux_for_field, field_gradient); } -template -void Fluxes::apply( - const gsl::not_null*> flux_for_field, +template +void Fluxes::apply( + const gsl::not_null*> flux_for_field, const tnsr::i& /*face_normal*/, const tnsr::I& face_normal_vector, - const Scalar& field) { + const Scalar& field) { fluxes_on_face(flux_for_field, face_normal_vector, field); } -template -void Fluxes::apply( - const gsl::not_null*> flux_for_field, +template +void Fluxes::apply( + const gsl::not_null*> flux_for_field, const tnsr::II& inv_spatial_metric, - const Scalar& /*field*/, - const tnsr::i& field_gradient) { + const Scalar& /*field*/, + const tnsr::i& field_gradient) { curved_fluxes(flux_for_field, inv_spatial_metric, field_gradient); } -template -void Fluxes::apply( - const gsl::not_null*> flux_for_field, +template +void Fluxes::apply( + const gsl::not_null*> flux_for_field, const tnsr::II& /*inv_spatial_metric*/, const tnsr::i& /*face_normal*/, const tnsr::I& face_normal_vector, - const Scalar& field) { + const Scalar& field) { fluxes_on_face(flux_for_field, face_normal_vector, field); } -template -void Sources::apply( - const gsl::not_null*> equation_for_field, +template +void Sources::apply( + const gsl::not_null*> equation_for_field, const tnsr::i& christoffel_contracted, - const Scalar& /*field*/, - const tnsr::I& field_flux) { + const Scalar& /*field*/, + const tnsr::I& field_flux) { add_curved_sources(equation_for_field, christoffel_contracted, field_flux); } } // namespace Poisson -#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) - -#define INSTANTIATE(_, data) \ - template void Poisson::flat_cartesian_fluxes( \ - const gsl::not_null*>, \ - const tnsr::i&); \ - template void Poisson::curved_fluxes( \ - const gsl::not_null*>, \ - const tnsr::II&, \ - const tnsr::i&); \ - template void Poisson::fluxes_on_face( \ - const gsl::not_null*>, \ - const tnsr::I&, const Scalar&); \ - template void Poisson::add_curved_sources( \ - const gsl::not_null*>, \ - const tnsr::i&, \ - const tnsr::I&); \ - template class Poisson::Fluxes; \ - template class Poisson::Fluxes; \ - template class Poisson::Sources; - -GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) +#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) +#define DIM(data) BOOST_PP_TUPLE_ELEM(1, data) + +#define INSTANTIATE(_, data) \ + template void Poisson::flat_cartesian_fluxes( \ + const gsl::not_null*>, \ + const tnsr::i&); \ + template void Poisson::curved_fluxes( \ + const gsl::not_null*>, \ + const tnsr::II&, \ + const tnsr::i&); \ + template void Poisson::fluxes_on_face( \ + const gsl::not_null*>, \ + const tnsr::I&, const Scalar&); \ + template void Poisson::add_curved_sources( \ + const gsl::not_null*>, \ + const tnsr::i&, \ + const tnsr::I&); \ + template class Poisson::Fluxes; \ + template class Poisson::Fluxes; \ + template class Poisson::Sources; + +GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector, ComplexDataVector), (1, 2, 3)) #undef INSTANTIATE #undef DIM diff --git a/src/Elliptic/Systems/Poisson/Equations.hpp b/src/Elliptic/Systems/Poisson/Equations.hpp index 41b04245ae11..c8a1e6c61505 100644 --- a/src/Elliptic/Systems/Poisson/Equations.hpp +++ b/src/Elliptic/Systems/Poisson/Equations.hpp @@ -18,9 +18,11 @@ namespace PUP { class er; } // namespace PUP namespace Poisson { -template +template struct Fluxes; -template +template struct Sources; } // namespace Poisson /// \endcond @@ -31,19 +33,19 @@ namespace Poisson { * \brief Compute the fluxes \f$F^i=\partial_i u(x)\f$ for the Poisson * equation on a flat spatial metric in Cartesian coordinates. */ -template +template void flat_cartesian_fluxes( - gsl::not_null*> flux_for_field, - const tnsr::i& field_gradient); + gsl::not_null*> flux_for_field, + const tnsr::i& field_gradient); /*! * \brief Compute the fluxes \f$F^i=\gamma^{ij}\partial_j u(x)\f$ * for the curved-space Poisson equation on a spatial metric \f$\gamma_{ij}\f$. */ -template -void curved_fluxes(gsl::not_null*> flux_for_field, +template +void curved_fluxes(gsl::not_null*> flux_for_field, const tnsr::II& inv_spatial_metric, - const tnsr::i& field_gradient); + const tnsr::i& field_gradient); /*! * \brief Compute the fluxes $F^i=\gamma^{ij} n_j u$ where $n_j$ is the @@ -51,10 +53,10 @@ void curved_fluxes(gsl::not_null*> flux_for_field, * * The `face_normal_vector` is $\gamma^{ij} n_j$. */ -template -void fluxes_on_face(gsl::not_null*> flux_for_field, +template +void fluxes_on_face(gsl::not_null*> flux_for_field, const tnsr::I& face_normal_vector, - const Scalar& field); + const Scalar& field); /*! * \brief Add the sources \f$S=-\Gamma^i_{ij}v^j\f$ @@ -63,10 +65,10 @@ void fluxes_on_face(gsl::not_null*> flux_for_field, * These sources arise from the non-principal part of the Laplacian on a * non-Euclidean background. */ -template -void add_curved_sources(gsl::not_null*> source_for_field, +template +void add_curved_sources(gsl::not_null*> source_for_field, const tnsr::i& christoffel_contracted, - const tnsr::I& flux_for_field); + const tnsr::I& flux_for_field); /*! * \brief Compute the fluxes \f$F^i\f$ for the Poisson equation on a flat @@ -74,20 +76,20 @@ void add_curved_sources(gsl::not_null*> source_for_field, * * \see Poisson::FirstOrderSystem */ -template -struct Fluxes { +template +struct Fluxes { using argument_tags = tmpl::list<>; using volume_tags = tmpl::list<>; using const_global_cache_tags = tmpl::list<>; static constexpr bool is_trivial = true; static constexpr bool is_discontinuous = false; - static void apply(gsl::not_null*> flux_for_field, - const Scalar& field, - const tnsr::i& field_gradient); - static void apply(gsl::not_null*> flux_for_field, + static void apply(gsl::not_null*> flux_for_field, + const Scalar& field, + const tnsr::i& field_gradient); + static void apply(gsl::not_null*> flux_for_field, const tnsr::i& face_normal, const tnsr::I& face_normal_vector, - const Scalar& field); + const Scalar& field); }; /*! @@ -96,23 +98,23 @@ struct Fluxes { * * \see Poisson::FirstOrderSystem */ -template -struct Fluxes { +template +struct Fluxes { using argument_tags = tmpl::list>; using volume_tags = tmpl::list<>; using const_global_cache_tags = tmpl::list<>; static constexpr bool is_trivial = true; static constexpr bool is_discontinuous = false; - static void apply(gsl::not_null*> flux_for_field, + static void apply(gsl::not_null*> flux_for_field, const tnsr::II& inv_spatial_metric, - const Scalar& field, - const tnsr::i& field_gradient); - static void apply(gsl::not_null*> flux_for_field, + const Scalar& field, + const tnsr::i& field_gradient); + static void apply(gsl::not_null*> flux_for_field, const tnsr::II& inv_spatial_metric, const tnsr::i& face_normal, const tnsr::I& face_normal_vector, - const Scalar& field); + const Scalar& field); }; /*! @@ -121,14 +123,14 @@ struct Fluxes { * * \see Poisson::FirstOrderSystem */ -template -struct Sources { +template +struct Sources { using argument_tags = tmpl::list< gr::Tags::SpatialChristoffelSecondKindContracted>; - static void apply(gsl::not_null*> equation_for_field, + static void apply(gsl::not_null*> equation_for_field, const tnsr::i& christoffel_contracted, - const Scalar& field, - const tnsr::I& field_flux); + const Scalar& field, + const tnsr::I& field_flux); }; } // namespace Poisson diff --git a/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp b/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp index 04e005ce8be9..43282d66da88 100644 --- a/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp @@ -53,20 +53,30 @@ namespace Poisson { * The fluxes and sources simplify significantly when the background metric is * flat and we employ Cartesian coordinates so \f$\gamma_{ij} = \delta_{ij}\f$ * and \f$\Gamma^i_{jk} = 0\f$. Set the template parameter `BackgroundGeometry` - * to `Poisson::Geometry::FlatCartesian` to specialise the system for this case. + * to `Poisson::Geometry::FlatCartesian` to specialize the system for this case. * Set it to `Poisson::Geometry::Curved` for the general case. + * + * ## Complex Poisson equation + * + * This system can also be used to solve the complex Poisson equation where + * $u(x)$ and $f(x)$ are complex-valued, by setting the `DataType` template + * parameter to `ComplexDataVector`. Note that the real and imaginary sectors of + * the equations decouple, so they are essentially two independent Poisson + * equations. This is useful for testing the elliptic solver with complex-valued + * fields, and also as building blocks for other Poisson-like systems of + * equations that have additional complex-valued source terms. */ -template +template struct FirstOrderSystem : tt::ConformsTo { static constexpr size_t volume_dim = Dim; - using primal_fields = tmpl::list>; + using primal_fields = tmpl::list>; // We just use the standard `Flux` prefix because the fluxes don't have // symmetries and we don't need to give them a particular meaning. - using primal_fluxes = - tmpl::list<::Tags::Flux, tmpl::size_t, - Frame::Inertial>>; + using primal_fluxes = tmpl::list< + ::Tags::Flux, tmpl::size_t, Frame::Inertial>>; using background_fields = tmpl::conditional_t< BackgroundGeometry == Geometry::FlatCartesian, tmpl::list<>, @@ -77,10 +87,10 @@ struct FirstOrderSystem tmpl::conditional_t>; - using fluxes_computer = Fluxes; + using fluxes_computer = Fluxes; using sources_computer = tmpl::conditional_t>; + Sources>; using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition; diff --git a/tests/Unit/Elliptic/Systems/Poisson/Test_Equations.cpp b/tests/Unit/Elliptic/Systems/Poisson/Test_Equations.cpp index 794dae6c4444..6807d229b47c 100644 --- a/tests/Unit/Elliptic/Systems/Poisson/Test_Equations.cpp +++ b/tests/Unit/Elliptic/Systems/Poisson/Test_Equations.cpp @@ -22,23 +22,23 @@ namespace helpers = TestHelpers::elliptic; namespace { -template +template void test_equations(const DataVector& used_for_size) { - pypp::check_with_random_values<1>(&Poisson::flat_cartesian_fluxes, - "Equations", {"flat_cartesian_fluxes"}, + pypp::check_with_random_values<1>( + &Poisson::flat_cartesian_fluxes, "Equations", + {"flat_cartesian_fluxes"}, {{{0., 1.}}}, used_for_size); + pypp::check_with_random_values<1>(&Poisson::curved_fluxes, + "Equations", {"curved_fluxes"}, {{{0., 1.}}}, used_for_size); - pypp::check_with_random_values<1>(&Poisson::curved_fluxes, "Equations", - {"curved_fluxes"}, {{{0., 1.}}}, - used_for_size); pypp::check_with_random_values<1>( - &Poisson::add_curved_sources, "Equations", {"add_curved_sources"}, - {{{0., 1.}}}, used_for_size, 1.e-12, {}, 0.); + &Poisson::add_curved_sources, "Equations", + {"add_curved_sources"}, {{{0., 1.}}}, used_for_size, 1.e-12, {}, 0.); } -template +template void test_computers(const DataVector& used_for_size) { CAPTURE(Dim); - using system = Poisson::FirstOrderSystem; + using system = Poisson::FirstOrderSystem; static_assert( tt::assert_conforms_to_v); helpers::test_first_order_fluxes_computer(used_for_size); @@ -52,8 +52,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.Systems.Poisson", "[Unit][Elliptic]") { "Elliptic/Systems/Poisson"}; GENERATE_UNINITIALIZED_DATAVECTOR; - CHECK_FOR_DATAVECTORS(test_equations, (1, 2, 3)); + CHECK_FOR_DATAVECTORS(test_equations, (DataVector, ComplexDataVector), + (1, 2, 3)); CHECK_FOR_DATAVECTORS( - test_computers, (1, 2, 3), + test_computers, (DataVector, ComplexDataVector), (1, 2, 3), (Poisson::Geometry::FlatCartesian, Poisson::Geometry::Curved)); }