diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/CMakeLists.txt b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/CMakeLists.txt index 9a3465a765ad..2755f019a9f8 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/CMakeLists.txt @@ -12,6 +12,7 @@ spectre_target_headers( ${LIBRARY} INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src HEADERS + Factory.hpp InitialMagneticField.hpp Poloidal.hpp Toroidal.hpp diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp new file mode 100644 index 000000000000..85c62da5d7be --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp @@ -0,0 +1,8 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp" diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp index b9ffda4eb340..d6262ce68509 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp @@ -6,7 +6,22 @@ #include #include +#include "DataStructures/Tensor/TypeAliases.hpp" #include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/TMPL.hpp" + +/// \cond +class DataVector; +namespace grmhd::AnalyticData::InitialMagneticFields { +class Poloidal; +class Toroidal; +} // namespace grmhd::AnalyticData::InitialMagneticFields + +namespace gsl { +template +class not_null; +} // namespace gsl +/// \endcond namespace grmhd::AnalyticData { @@ -32,21 +47,47 @@ namespace grmhd::AnalyticData { * B^z & = \frac{1}{\sqrt{\gamma}} (\partial_x A_y - \partial_y A_x). * \f} * + * + * \warning The magnetic field classes assume the magnetic field is initialized, + * both in size and value, before being passed into the `variables` function. + * This is so that multiple magnetic fields can be superposed. Each magnetic + * field configuration does a `+=` to make this possible. */ namespace InitialMagneticFields { /*! * \brief The abstract base class for initial magnetic field configurations. + * + * \warning This assumes the magnetic field is initialized, both in size and + * value, before being passed into the `variables` function. This is so that + * multiple magnetic fields can be superposed. Each magnetic field + * configuration does a `+=` to make this possible. */ class InitialMagneticField : public PUP::able { protected: InitialMagneticField() = default; public: + using creatable_classes = tmpl::list; + ~InitialMagneticField() override = default; virtual auto get_clone() const -> std::unique_ptr = 0; + virtual void variables( + gsl::not_null*> result, + const tnsr::I& coords, const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const = 0; + + virtual void variables(gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const = 0; + + virtual bool is_equal(const InitialMagneticField& rhs) const = 0; + /// \cond explicit InitialMagneticField(CkMigrateMessage* msg) : PUP::able(msg) {} WRAPPED_PUPable_abstract(InitialMagneticField); diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.cpp index 0dbaa6f9fe7c..d9b5943b89f7 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.cpp @@ -24,72 +24,103 @@ void Poloidal::pup(PUP::er& p) { p | pressure_exponent_; p | cutoff_pressure_; p | vector_potential_amplitude_; + p | center_; + p | max_distance_from_center_; } // NOLINTNEXTLINE PUP::able::PUP_ID Poloidal::my_PUP_ID = 0; Poloidal::Poloidal(const size_t pressure_exponent, const double cutoff_pressure, - const double vector_potential_amplitude) + const double vector_potential_amplitude, + const std::array center, + const double max_distance_from_center) : pressure_exponent_(pressure_exponent), cutoff_pressure_(cutoff_pressure), - vector_potential_amplitude_(vector_potential_amplitude) {} + vector_potential_amplitude_(vector_potential_amplitude), + center_(center), + max_distance_from_center_(max_distance_from_center) {} + +void Poloidal::variables(const gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const { + ASSERT(result->get(0).size() == get(pressure).size(), + "Result must be of size " << get(pressure).size() << " but got " + << result->get(0).size()); + variables_impl(result, coords, pressure, sqrt_det_spatial_metric, + deriv_pressure); +} + +void Poloidal::variables(const gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const { + variables_impl(result, coords, pressure, sqrt_det_spatial_metric, + deriv_pressure); +} + +bool Poloidal::is_equal(const InitialMagneticField& rhs) const { + const auto& derived_ptr = dynamic_cast(&rhs); + return derived_ptr != nullptr and *derived_ptr == *this; +} template -tuples::TaggedTuple> -Poloidal::variables(const tnsr::I& coords, - const Scalar& pressure, - const Scalar& sqrt_det_spatial_metric, - const tnsr::i& dcoords_pressure) const { - auto magnetic_field = make_with_value>(coords, 0.0); +void Poloidal::variables_impl( + const gsl::not_null*> magnetic_field, + const tnsr::I& coords, const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const { const size_t num_pts = get_size(get(pressure)); for (size_t i = 0; i < num_pts; ++i) { const double pressure_i = get_element(get(pressure), i); - if (pressure_i < cutoff_pressure_) { - get_element(magnetic_field.get(0), i) = 0.0; - get_element(magnetic_field.get(1), i) = 0.0; - get_element(magnetic_field.get(2), i) = 0.0; + const double x = get_element(coords.get(0), i) - center_[0]; + const double y = get_element(coords.get(1), i) - center_[1]; + const double z = get_element(coords.get(2), i) - center_[2]; + const double radius = sqrt(x * x + y * y + z * z); + if (pressure_i < cutoff_pressure_ or radius > max_distance_from_center_) { continue; } // (p - p_c)^{n_s} - const double pressure_term = - pow(pressure_i - cutoff_pressure_, pressure_exponent_); + const double pressure_term = pow(pressure_i - cutoff_pressure_, + static_cast(pressure_exponent_)); // n_s * (p - p_c)^{n_s-1} const double n_times_pressure_to_n_minus_1 = - pressure_exponent_ * pow(pressure_i - cutoff_pressure_, - static_cast(pressure_exponent_) - 1); + static_cast(pressure_exponent_) * + pow(pressure_i - cutoff_pressure_, + static_cast(pressure_exponent_) - 1); - const double x = get_element(coords.get(0), i); - const double y = get_element(coords.get(1), i); - - const auto& dp_dx = get_element(dcoords_pressure.get(0), i); - const auto& dp_dy = get_element(dcoords_pressure.get(1), i); - const auto& dp_dz = get_element(dcoords_pressure.get(2), i); + const auto& dp_dx = get_element(deriv_pressure.get(0), i); + const auto& dp_dy = get_element(deriv_pressure.get(1), i); + const auto& dp_dz = get_element(deriv_pressure.get(2), i); // Assign Bx, By, Bz - get_element(magnetic_field.get(0), i) = - -n_times_pressure_to_n_minus_1 * x * dp_dz; - get_element(magnetic_field.get(1), i) = - -n_times_pressure_to_n_minus_1 * y * dp_dz; - get_element(magnetic_field.get(2), i) = - 2.0 * pressure_term + - n_times_pressure_to_n_minus_1 * (x * dp_dx + y * dp_dy); - } - - for (size_t d = 0; d < 3; ++d) { - magnetic_field.get(d) *= - vector_potential_amplitude_ / get(sqrt_det_spatial_metric); + get_element(magnetic_field->get(0), i) += + vector_potential_amplitude_ * + (-n_times_pressure_to_n_minus_1 * x * dp_dz) / + get_element(get(sqrt_det_spatial_metric), i); + get_element(magnetic_field->get(1), i) += + vector_potential_amplitude_ * + (-n_times_pressure_to_n_minus_1 * y * dp_dz) / + get_element(get(sqrt_det_spatial_metric), i); + get_element(magnetic_field->get(2), i) += + vector_potential_amplitude_ * + (2.0 * pressure_term + + n_times_pressure_to_n_minus_1 * (x * dp_dx + y * dp_dy)) / + get_element(get(sqrt_det_spatial_metric), i); } - - return {std::move(magnetic_field)}; } bool operator==(const Poloidal& lhs, const Poloidal& rhs) { return lhs.pressure_exponent_ == rhs.pressure_exponent_ and lhs.cutoff_pressure_ == rhs.cutoff_pressure_ and - lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_; + lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_ and + lhs.center_ == rhs.center_ and + lhs.max_distance_from_center_ == rhs.max_distance_from_center_; } bool operator!=(const Poloidal& lhs, const Poloidal& rhs) { @@ -98,13 +129,13 @@ bool operator!=(const Poloidal& lhs, const Poloidal& rhs) { #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template tuples::TaggedTuple> \ - Poloidal::variables( \ - const tnsr::I& coords, \ - const Scalar& pressure, \ - const Scalar& sqrt_det_spatial_metric, \ - const tnsr::i& dcoords_pressure) const; +#define INSTANTIATE(_, data) \ + template void Poloidal::variables_impl( \ + gsl::not_null*> magnetic_field, \ + const tnsr::I& coords, \ + const Scalar& pressure, \ + const Scalar& sqrt_det_spatial_metric, \ + const tnsr::i& deriv_pressure) const; GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector)) diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp index dec1dfd99036..67a056b8668b 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp @@ -14,6 +14,13 @@ #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" +/// \cond +namespace gsl { +template +class not_null; +} // namespace gsl +/// \endcond + namespace grmhd::AnalyticData::InitialMagneticFields { /*! @@ -63,6 +70,15 @@ namespace grmhd::AnalyticData::InitialMagneticFields { * 2(p-p_{\mathrm{cut}})^{n_s}. * \f} * + * Note that the coordinates are relative to the `Center` passed in, so the + * field can be centered about any arbitrary point. The field is also zero + * outside of `MaxDistanceFromCenter`, so that compact support can be imposed if + * necessary. + * + * \warning This assumes the magnetic field is initialized, both in size and + * value, before being passed into the `variables` function. This is so that + * multiple magnetic fields can be superposed. Each magnetic field + * configuration does a `+=` to make this possible. */ class Poloidal : public InitialMagneticField { public: @@ -87,8 +103,23 @@ class Poloidal : public InitialMagneticField { static type lower_bound() { return 0.0; } }; + struct Center { + using type = std::array; + static constexpr Options::String help = { + "The center of the magnetic field."}; + }; + + struct MaxDistanceFromCenter { + using type = double; + static constexpr Options::String help = { + "The maximum distance from the center to compute the magnetic field. " + "Everywhere outside the field is set to zero."}; + static type lower_bound() { return 0.0; } + }; + using options = - tmpl::list; + tmpl::list; static constexpr Options::String help = {"Poloidal initial magnetic field"}; @@ -100,7 +131,8 @@ class Poloidal : public InitialMagneticField { ~Poloidal() override = default; Poloidal(size_t pressure_exponent, double cutoff_pressure, - double vector_potential_amplitude); + double vector_potential_amplitude, std::array center, + double max_distance_from_center); auto get_clone() const -> std::unique_ptr override; @@ -114,18 +146,38 @@ class Poloidal : public InitialMagneticField { void pup(PUP::er& p) override; /// Retrieve magnetic fields at `(x)` - template - auto variables(const tnsr::I& coords, - const Scalar& pressure, - const Scalar& sqrt_det_spatial_metric, - const tnsr::i& dcoords_pressure) const - -> tuples::TaggedTuple>; + void variables(gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const override; + + /// Retrieve magnetic fields at `(x)` + void variables(gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const override; + + bool is_equal(const InitialMagneticField& rhs) const override; private: + template + void variables_impl(gsl::not_null*> magnetic_field, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const; + size_t pressure_exponent_ = std::numeric_limits::max(); double cutoff_pressure_ = std::numeric_limits::signaling_NaN(); double vector_potential_amplitude_ = std::numeric_limits::signaling_NaN(); + std::array center_{{std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN()}}; + double max_distance_from_center_ = + std::numeric_limits::signaling_NaN(); friend bool operator==(const Poloidal& lhs, const Poloidal& rhs); friend bool operator!=(const Poloidal& lhs, const Poloidal& rhs); diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.cpp index f050a3e08172..920230781413 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.cpp @@ -25,71 +25,102 @@ void Toroidal::pup(PUP::er& p) { p | pressure_exponent_; p | cutoff_pressure_; p | vector_potential_amplitude_; + p | center_; + p | max_distance_from_center_; } // NOLINTNEXTLINE PUP::able::PUP_ID Toroidal::my_PUP_ID = 0; Toroidal::Toroidal(const size_t pressure_exponent, const double cutoff_pressure, - const double vector_potential_amplitude) + const double vector_potential_amplitude, + const std::array center, + const double max_distance_from_center) : pressure_exponent_(pressure_exponent), cutoff_pressure_(cutoff_pressure), - vector_potential_amplitude_(vector_potential_amplitude) {} + vector_potential_amplitude_(vector_potential_amplitude), + center_(center), + max_distance_from_center_(max_distance_from_center) {} + +void Toroidal::variables(const gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const { + ASSERT(result->get(0).size() == get(pressure).size(), + "Result must be of size " << get(pressure).size() << " but got " + << result->get(0).size()); + variables_impl(result, coords, pressure, sqrt_det_spatial_metric, + deriv_pressure); +} + +void Toroidal::variables(const gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const { + variables_impl(result, coords, pressure, sqrt_det_spatial_metric, + deriv_pressure); +} template -tuples::TaggedTuple> -Toroidal::variables(const tnsr::I& coords, - const Scalar& pressure, - const Scalar& sqrt_det_spatial_metric, - const tnsr::i& deriv_pressure) const { - auto magnetic_field = make_with_value>(coords, 0.0); +void Toroidal::variables_impl( + const gsl::not_null*> magnetic_field, + const tnsr::I& coords, const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const { const size_t num_pts = get_size(get(pressure)); for (size_t i = 0; i < num_pts; ++i) { const double pressure_i = get_element(get(pressure), i); - if (pressure_i < cutoff_pressure_) { - get_element(magnetic_field.get(0), i) = 0.0; - get_element(magnetic_field.get(1), i) = 0.0; - get_element(magnetic_field.get(2), i) = 0.0; + const double x = get_element(coords.get(0), i) - center_[0]; + const double y = get_element(coords.get(1), i) - center_[1]; + const double z = get_element(coords.get(2), i) - center_[2]; + const double radius = sqrt(x * x + y * y + z * z); + if (pressure_i < cutoff_pressure_ or radius > max_distance_from_center_) { continue; } // (p - p_c)^{n_s} const double pressure_term = - 2.0 * pow(pressure_i - cutoff_pressure_, pressure_exponent_); + 2.0 * pow(pressure_i - cutoff_pressure_, + static_cast(pressure_exponent_)); // n_s * (p - p_c)^{n_s-1} const double n_times_pressure_to_n_minus_1 = - pressure_exponent_ * pow(pressure_i - cutoff_pressure_, - static_cast(pressure_exponent_) - 1); + static_cast(pressure_exponent_) * + pow(pressure_i - cutoff_pressure_, + static_cast(pressure_exponent_) - 1); - const double x = get_element(coords.get(0), i); - const double y = get_element(coords.get(1), i); const double varpi_squared = square(x) + square(y); const auto& dp_dx = get_element(deriv_pressure.get(0), i); const auto& dp_dy = get_element(deriv_pressure.get(1), i); // Assign Bx, By - get_element(magnetic_field.get(0), i) = - y * pressure_term + - varpi_squared * n_times_pressure_to_n_minus_1 * dp_dy; - get_element(magnetic_field.get(1), i) = - -(x * pressure_term + - varpi_squared * n_times_pressure_to_n_minus_1 * dp_dx); + get_element(magnetic_field->get(0), i) += + vector_potential_amplitude_ * + (y * pressure_term + + varpi_squared * n_times_pressure_to_n_minus_1 * dp_dy) / + get_element(get(sqrt_det_spatial_metric), i); + get_element(magnetic_field->get(1), i) += + vector_potential_amplitude_ * + (-(x * pressure_term + + varpi_squared * n_times_pressure_to_n_minus_1 * dp_dx)) / + get_element(get(sqrt_det_spatial_metric), i); } +} - magnetic_field.get(0) *= - vector_potential_amplitude_ / get(sqrt_det_spatial_metric); - magnetic_field.get(1) *= - vector_potential_amplitude_ / get(sqrt_det_spatial_metric); - - return {std::move(magnetic_field)}; +bool Toroidal::is_equal(const InitialMagneticField& rhs) const { + const auto& derived_ptr = dynamic_cast(&rhs); + return derived_ptr != nullptr and *derived_ptr == *this; } bool operator==(const Toroidal& lhs, const Toroidal& rhs) { return lhs.pressure_exponent_ == rhs.pressure_exponent_ and lhs.cutoff_pressure_ == rhs.cutoff_pressure_ and - lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_; + lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_ and + lhs.center_ == rhs.center_ and + lhs.max_distance_from_center_ == rhs.max_distance_from_center_; } bool operator!=(const Toroidal& lhs, const Toroidal& rhs) { @@ -98,12 +129,12 @@ bool operator!=(const Toroidal& lhs, const Toroidal& rhs) { #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template tuples::TaggedTuple> \ - Toroidal::variables( \ - const tnsr::I& coords, \ - const Scalar& pressure, \ - const Scalar& sqrt_det_spatial_metric, \ +#define INSTANTIATE(_, data) \ + template void Toroidal::variables_impl( \ + gsl::not_null*> magnetic_field, \ + const tnsr::I& coords, \ + const Scalar& pressure, \ + const Scalar& sqrt_det_spatial_metric, \ const tnsr::i& deriv_pressure) const; GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector)) diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp index 0446a5014e52..412935404524 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp @@ -44,6 +44,15 @@ namespace grmhd::AnalyticData::InitialMagneticFields { * B^z & = 0 . * \f} * + * Note that the coordinates are relative to the `Center` passed in, so the + * field can be centered about any arbitrary point. The field is also zero + * outside of `MaxDistanceFromCenter`, so that compact support can be imposed if + * necessary. + * + * \warning This assumes the magnetic field is initialized, both in size and + * value, before being passed into the `variables` function. This is so that + * multiple magnetic fields can be superposed. Each magnetic field + * configuration does a `+=` to make this possible. */ class Toroidal : public InitialMagneticField { public: @@ -68,8 +77,23 @@ class Toroidal : public InitialMagneticField { static type lower_bound() { return 0.0; } }; + struct Center { + using type = std::array; + static constexpr Options::String help = { + "The center of the magnetic field."}; + }; + + struct MaxDistanceFromCenter { + using type = double; + static constexpr Options::String help = { + "The maximum distance from the center to compute the magnetic field. " + "Everywhere outside the field is set to zero."}; + static type lower_bound() { return 0.0; } + }; + using options = - tmpl::list; + tmpl::list; static constexpr Options::String help = {"Toroidal initial magnetic field"}; @@ -81,7 +105,8 @@ class Toroidal : public InitialMagneticField { ~Toroidal() override = default; Toroidal(size_t pressure_exponent, double cutoff_pressure, - double vector_potential_amplitude); + double vector_potential_amplitude, std::array center, + double max_distance_from_center); auto get_clone() const -> std::unique_ptr override; @@ -95,18 +120,38 @@ class Toroidal : public InitialMagneticField { void pup(PUP::er& p) override; /// Retrieve magnetic fields at `(x)` - template - auto variables(const tnsr::I& coords, - const Scalar& pressure, - const Scalar& sqrt_det_spatial_metric, - const tnsr::i& deriv_pressure) const - -> tuples::TaggedTuple>; + void variables(gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const override; + + /// Retrieve magnetic fields at `(x)` + void variables(gsl::not_null*> result, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const override; + + bool is_equal(const InitialMagneticField& rhs) const override; private: + template + void variables_impl(gsl::not_null*> magnetic_field, + const tnsr::I& coords, + const Scalar& pressure, + const Scalar& sqrt_det_spatial_metric, + const tnsr::i& deriv_pressure) const; + size_t pressure_exponent_ = std::numeric_limits::max(); double cutoff_pressure_ = std::numeric_limits::signaling_NaN(); double vector_potential_amplitude_ = std::numeric_limits::signaling_NaN(); + std::array center_{{std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN()}}; + double max_distance_from_center_ = + std::numeric_limits::signaling_NaN(); friend bool operator==(const Toroidal& lhs, const Toroidal& rhs); }; diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp index 74bd7571213d..91df10162ea8 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp @@ -16,6 +16,8 @@ #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/Algorithm.hpp" +#include "Utilities/CloneUniquePtrs.hpp" #include "Utilities/ConstantExpressions.hpp" #include "Utilities/ContainerHelpers.hpp" #include "Utilities/ErrorHandling/Assert.hpp" @@ -23,20 +25,39 @@ #include "Utilities/MakeWithValue.hpp" namespace grmhd::AnalyticData { +MagnetizedTovStar::MagnetizedTovStar() = default; +MagnetizedTovStar::MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) = default; +MagnetizedTovStar& MagnetizedTovStar::operator=(MagnetizedTovStar&& /*rhs*/) = + default; +MagnetizedTovStar::~MagnetizedTovStar() = default; + MagnetizedTovStar::MagnetizedTovStar( const double central_rest_mass_density, std::unique_ptr equation_of_state, const RelativisticEuler::Solutions::TovCoordinates coordinate_system, - const size_t pressure_exponent, const double cutoff_pressure_fraction, - const double vector_potential_amplitude) + std::vector> + magnetic_fields) : tov_star(central_rest_mass_density, std::move(equation_of_state), coordinate_system), - pressure_exponent_(pressure_exponent), - cutoff_pressure_(cutoff_pressure_fraction * - get(this->equation_of_state().pressure_from_density( - Scalar{central_rest_mass_density}))), - vector_potential_amplitude_(vector_potential_amplitude) {} + magnetic_fields_(std::move(magnetic_fields)) {} + +MagnetizedTovStar::MagnetizedTovStar(const MagnetizedTovStar& rhs) + : evolution::initial_data::InitialData{rhs}, + RelativisticEuler::Solutions::TovStar( + static_cast(rhs)), + magnetic_fields_(clone_unique_ptrs(rhs.magnetic_fields_)) {} + +MagnetizedTovStar& MagnetizedTovStar::operator=(const MagnetizedTovStar& rhs) { + if (this == &rhs) { + return *this; + } + static_cast(*this) = + static_cast(rhs); + magnetic_fields_ = clone_unique_ptrs(rhs.magnetic_fields_); + return *this; +} std::unique_ptr MagnetizedTovStar::get_clone() const { @@ -47,9 +68,7 @@ MagnetizedTovStar::MagnetizedTovStar(CkMigrateMessage* msg) : tov_star(msg) {} void MagnetizedTovStar::pup(PUP::er& p) { tov_star::pup(p); - p | pressure_exponent_; - p | cutoff_pressure_; - p | vector_potential_amplitude_; + p | magnetic_fields_; } namespace magnetized_tov_detail { @@ -58,56 +77,21 @@ void MagnetizedTovVariables::operator()( const gsl::not_null*> magnetic_field, const gsl::not_null cache, hydro::Tags::MagneticField /*meta*/) const { - const size_t num_pts = get_size(get<0>(coords)); - const auto& pressure_profile = - get(cache->get_var(*this, hydro::Tags::Pressure{})); - const auto& dr_pressure_profile = get(cache->get_var( - *this, - RelativisticEuler::Solutions::tov_detail::Tags::DrPressure{})); - const auto& sqrt_det_spatial_metric = - get(cache->get_var(*this, gr::Tags::SqrtDetSpatialMetric{})); - for (size_t i = 0; i < num_pts; ++i) { - const double pressure = get_element(pressure_profile, i); - if (LIKELY(get_element(radius, i) > 1.0e-16)) { - if (pressure < cutoff_pressure) { - get_element(get<0>(*magnetic_field), i) = 0.0; - get_element(get<1>(*magnetic_field), i) = 0.0; - get_element(get<2>(*magnetic_field), i) = 0.0; - continue; - } - - const double x = get_element(get<0>(coords), i); - const double y = get_element(get<1>(coords), i); - const double z = get_element(get<2>(coords), i); - const double radius_i = get_element(radius, i); - const double dr_pressure = get_element(dr_pressure_profile, i); - const double pressure_term = - pow(pressure - cutoff_pressure, pressure_exponent); - const double deriv_pressure_term = - pressure_exponent * - pow(pressure - cutoff_pressure, - static_cast(pressure_exponent) - 1) * - dr_pressure; - - get_element(get<0>(*magnetic_field), i) = - -x * z / radius_i * deriv_pressure_term; - - get_element(get<1>(*magnetic_field), i) = - -y * z / radius_i * deriv_pressure_term; - - get_element(get<2>(*magnetic_field), i) = - 2.0 * pressure_term + - (square(x) + square(y)) / radius_i * deriv_pressure_term; - } else { - get_element(get<0>(*magnetic_field), i) = 0.0; - get_element(get<1>(*magnetic_field), i) = 0.0; - get_element(get<2>(*magnetic_field), i) = - 2.0 * pow(pressure - cutoff_pressure, pressure_exponent); - } - } + const Scalar& sqrt_det_spatial_metric = + cache->get_var(*this, gr::Tags::SqrtDetSpatialMetric{}); + const Scalar& pressure = + cache->get_var(*this, hydro::Tags::Pressure{}); + const auto& deriv_pressure = + cache->get_var(*this, ::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>{}); + set_number_of_grid_points(magnetic_field, get_size(get<0>(coords))); for (size_t i = 0; i < 3; ++i) { - magnetic_field->get(i) *= - vector_potential_amplitude / sqrt_det_spatial_metric; + magnetic_field->get(i) = 0.0; + } + + for (const auto& magnetic_field_compute : magnetic_fields) { + magnetic_field_compute->variables(magnetic_field, coords, pressure, + sqrt_det_spatial_metric, deriv_pressure); } } } // namespace magnetized_tov_detail @@ -115,11 +99,19 @@ void MagnetizedTovVariables::operator()( PUP::able::PUP_ID MagnetizedTovStar::my_PUP_ID = 0; bool operator==(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs) { - return static_cast(lhs) == - static_cast(rhs) and - lhs.pressure_exponent_ == rhs.pressure_exponent_ and - lhs.cutoff_pressure_ == rhs.cutoff_pressure_ and - lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_; + bool equal = + static_cast(lhs) == + static_cast(rhs) and + lhs.magnetic_fields_.size() == rhs.magnetic_fields_.size(); + if (not equal) { + return false; + } + for (size_t i = 0; i < lhs.magnetic_fields_.size(); ++i) { + if (not lhs.magnetic_fields_[i]->is_equal(*rhs.magnetic_fields_[i])) { + return false; + } + } + return true; } bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs) { diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp index ed35ba27019b..5a887b40ac0c 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp @@ -10,6 +10,8 @@ #include "Options/String.hpp" #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp" #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp" #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp" #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp" @@ -42,20 +44,19 @@ struct MagnetizedTovVariables using Base::radial_solution; using Base::radius; - size_t pressure_exponent; - double cutoff_pressure; - double vector_potential_amplitude; + const std::vector>& + magnetic_fields; MagnetizedTovVariables( const tnsr::I& local_x, const DataType& local_radius, const RelativisticEuler::Solutions::TovSolution& local_radial_solution, const EquationsOfState::EquationOfState& local_eos, - size_t local_pressure_exponent, double local_cutoff_pressure, - double local_vector_potential_amplitude) + const std::vector>& + mag_fields) : Base(local_x, local_radius, local_radial_solution, local_eos), - pressure_exponent(local_pressure_exponent), - cutoff_pressure(local_cutoff_pressure), - vector_potential_amplitude(local_vector_potential_amplitude) {} + magnetic_fields(mag_fields) {} void operator()( gsl::not_null*> magnetic_field, @@ -69,86 +70,11 @@ struct MagnetizedTovVariables * \brief Magnetized TOV star initial data, where metric terms only account for * the hydrodynamics not the magnetic fields. * - * Superposes a poloidal magnetic field on top of a TOV solution where the - * vector potential has the form + * Superposes magnetic fields on top of a TOV solution. These can be any of + * the classes derived from + * grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField * - * \f{align*}{ - * A_{\phi} = A_b \varpi^2 \max(p-p_{\mathrm{cut}}, 0)^{n_s} - * \f} - * - * where \f$A_b\f$ controls the amplitude of the magnetic field, - * \f$\varpi^2=x^2+y^2=r^2-z^2\f$ is the cylindrical radius, - * \f$n_s\f$ controls the degree of differentiability, and - * \f$p_{\mathrm{cut}}\f$ controls the pressure cutoff below which the magnetic - * field is zero. - * - * In Cartesian coordinates the vector potential is: - * - * \f{align*}{ - * A_x&=-\frac{y}{\varpi^2}A_\phi, \\ - * A_y&=\frac{x}{\varpi^2}A_\phi, \\ - * A_z&=0, - * \f} - * - * For the poloidal field this means - * - * \f{align*}{ - * A_x &= -y A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\ - * A_y &= x A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}. - * \f} - * - * The magnetic field is computed from the vector potential as - * - * \f{align*}{ - * B^i&=n_a\epsilon^{aijk}\partial_jA_k \\ - * &=\frac{1}{\sqrt{\gamma}}[ijk]\partial_j A_k, - * \f} - * - * where \f$[ijk]\f$ is the total anti-symmetric symbol. This means that - * - * \f{align*}{ - * B^x&=\frac{1}{\sqrt{\gamma}} (\partial_y A_z - \partial_z A_y), \\ - * B^y&=\frac{1}{\sqrt{\gamma}} (\partial_z A_x - \partial_x A_z), \\ - * B^z&=\frac{1}{\sqrt{\gamma}} (\partial_x A_y - \partial_y A_x). - * \f} - * - * Focusing on the region where the field is non-zero we have: - * - * \f{align*}{ - * \partial_x A_y - * &= A_b(p-p_{\mathrm{cut}})^{n_s}+ - * \frac{x^2}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p \\ - * \partial_y A_x - * &= -A_b(p-p_{\mathrm{cut}})^{n_s}- - * \frac{y^2}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p \\ - * \partial_z A_x - * &= -\frac{yz}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\ - * \partial_z A_y - * &= \frac{xz}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp - * \f} - * - * The magnetic field is given by: - * - * \f{align*}{ - * B^x&=-\frac{1}{\sqrt{\gamma}}\frac{xz}{r} - * A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\ - * B^y&=-\frac{1}{\sqrt{\gamma}}\frac{yz}{r} - * A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\ - * B^z&=\frac{A_b}{\sqrt{\gamma}}\left[ - * 2(p-p_{\mathrm{cut}})^{n_s} \phantom{\frac{a}{b}}\right. \\ - * &\left.+\frac{x^2+y^2}{r} - * n_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p - * \right] - * \f} - * - * Taking the small-\f$r\f$ limit gives the \f$r=0\f$ magnetic field: - * - * \f{align*}{ - * B^x&=0, \\ - * B^y&=0, \\ - * B^z&=\frac{A_b}{\sqrt{\gamma}} - * 2(p-p_{\mathrm{cut}})^{n_s}. - * \f} + * ### Conversion to CGS units and values for poloidal magnetic field * * While the amplitude \f$A_b\f$ is specified in the code, it is more natural * to work with the magnetic field strength, which is given by \f$\sqrt{b^2}\f$ @@ -161,7 +87,8 @@ struct MagnetizedTovVariables * &= \sqrt{b^2} \times 8.352\times10^{19}\mathrm{G} \,. * \f} * - * We now give values used for standard tests of magnetized stars. + * We now give values used for standard tests of magnetized stars with a + * poloidal magnetic field. * - \f$\rho_c(0)=1.28\times10^{-3}\f$ * - \f$K=100\f$ * - \f$\Gamma=2\f$ @@ -194,32 +121,14 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData, using tov_star = RelativisticEuler::Solutions::TovStar; public: - struct PressureExponent { - using type = size_t; - static constexpr Options::String help = { - "The exponent n_s controlling the smoothness of the field"}; - }; - - struct VectorPotentialAmplitude { - using type = double; - static constexpr Options::String help = { - "The amplitude A_b of the phi-component of the vector potential. This " - "controls the magnetic field strength."}; - static type lower_bound() { return 0.0; } - }; - - struct CutoffPressureFraction { - using type = double; + struct MagneticFields { + using type = std::vector>; static constexpr Options::String help = { - "The fraction of the central pressure below which there is no magnetic " - "field. p_cut = Fraction * p_max."}; - static type lower_bound() { return 0.0; } - static type upper_bound() { return 1.0; } + "Magnetic fields to superpose on the TOV solution."}; }; - using options = - tmpl::push_back; + using options = tmpl::push_back; static constexpr Options::String help = {"Magnetized TOV star."}; @@ -228,20 +137,21 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData, template using tags = typename tov_star::template tags; - MagnetizedTovStar() = default; - MagnetizedTovStar(const MagnetizedTovStar& /*rhs*/) = default; - MagnetizedTovStar& operator=(const MagnetizedTovStar& /*rhs*/) = default; - MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) = default; - MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/) = default; - ~MagnetizedTovStar() override = default; + MagnetizedTovStar(); + MagnetizedTovStar(const MagnetizedTovStar& rhs); + MagnetizedTovStar& operator=(const MagnetizedTovStar& rhs); + MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/); + MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/); + ~MagnetizedTovStar() override; MagnetizedTovStar( double central_rest_mass_density, std::unique_ptr> equation_of_state, RelativisticEuler::Solutions::TovCoordinates coordinate_system, - size_t pressure_exponent, double cutoff_pressure_fraction, - double vector_potential_amplitude); + std::vector> + magnetic_fields); auto get_clone() const -> std::unique_ptr override; @@ -260,8 +170,7 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData, tuples::TaggedTuple variables(const tnsr::I& x, tmpl::list /*meta*/) const { return variables_impl( - x, tmpl::list{}, pressure_exponent_, cutoff_pressure_, - vector_potential_amplitude_); + x, tmpl::list{}, magnetic_fields_); } // NOLINTNEXTLINE(google-runtime-references) @@ -272,10 +181,9 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData, const MagnetizedTovStar& rhs); protected: - size_t pressure_exponent_ = std::numeric_limits::max(); - double cutoff_pressure_ = std::numeric_limits::signaling_NaN(); - double vector_potential_amplitude_ = - std::numeric_limits::signaling_NaN(); + std::vector> + magnetic_fields_{}; }; bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs); diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp index 639a4069a689..e3a3bd0d3132 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp @@ -36,6 +36,7 @@ #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" #include "Utilities/Gsl.hpp" +#include "Utilities/MakeVector.hpp" #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" @@ -336,9 +337,12 @@ SPECTRE_TEST_CASE( std::make_unique>(100.0, 2.0), RelativisticEuler::Solutions::TovCoordinates::Schwarzschild, - 2, - 0.04, - 2500.0}; + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))}; const auto serialized_and_deserialized_condition = serialize_and_deserialize( *dynamic_cast>(100.0, 2.0), RelativisticEuler::Solutions::TovCoordinates::Schwarzschild, - 2, - 0.04, - 2500.0}; + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))}; const auto serialized_and_deserialized_condition = serialize_and_deserialize( *dynamic_cast& x, const Scalar& pressure, const Scalar& sqrt_det_spatial_metric, const tnsr::i& deriv_pressure) const { - return this->variables(x, pressure, sqrt_det_spatial_metric, - deriv_pressure); + auto mag_field = make_with_value>(pressure, 0.0); + this->variables(make_not_null(&mag_field), x, pressure, + sqrt_det_spatial_metric, deriv_pressure); + return {mag_field}; } }; SPECTRE_TEST_CASE( "Unit.PointwiseFunctions.AnalyticData.GrMhd.InitialMagneticFields.Poloidal", "[Unit][PointwiseFunctions]") { + register_derived_classes_with_charm(); // test creation - const auto solution = TestHelpers::test_creation( - " PressureExponent: 2\n" - " CutoffPressure: 1.0e-5\n" - " VectorPotentialAmplitude: 2500.0\n"); - CHECK(solution == Poloidal(2, 1.0e-5, 2500.0)); + const auto factory_solution = serialize_and_deserialize( + TestHelpers::test_factory_creation( + "Poloidal:\n" + " PressureExponent: 2\n" + " CutoffPressure: 1.0e-5\n" + " VectorPotentialAmplitude: 2500.0\n" + " Center: [0.0, 0.0, 0.0]\n" + " MaxDistanceFromCenter: 100.0\n")); + CHECK(factory_solution->is_equal( + Poloidal(2, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0))); + + const auto& solution = dynamic_cast(*factory_solution); // test serialize test_serialization(solution); // test move { - Poloidal poloidal_field{2, 1.0e-5, 2500.0}; - Poloidal poloidal_field_copy{2, 1.0e-5, 2500.0}; + Poloidal poloidal_field{2, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0}; + Poloidal poloidal_field_copy{2, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0}; test_move_semantics(std::move(poloidal_field), poloidal_field_copy); } @@ -67,12 +78,20 @@ SPECTRE_TEST_CASE( CHECK(dynamic_cast(deserialized_base_ptr.get()) != nullptr); // test equality - const Poloidal field_original{2, 1.0e-5, 2500.0}; + const Poloidal field_original{2, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0}; const auto field = serialize_and_deserialize(field_original); - CHECK(field == Poloidal(2, 1.0e-5, 2500.0)); - CHECK(field != Poloidal(3, 1.0e-5, 2500.0)); - CHECK(field != Poloidal(2, 2.0e-5, 2500.0)); - CHECK(field != Poloidal(2, 1.0e-5, 3500.0)); + CHECK(static_cast(field).is_equal( + Poloidal(2, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Poloidal(3, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Poloidal(2, 2.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Poloidal(2, 1.0e-5, 3500.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Poloidal(2, 1.0e-5, 2500.0, {{1.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Poloidal(2, 1.0e-5, 2500.0, {{0.0, 0.0, 0.0}}, 10.0))); // test solution implementation pypp::SetupLocalPythonEnvironment local_python_env{ @@ -86,7 +105,7 @@ SPECTRE_TEST_CASE( pypp::check_with_random_values<1>( &PoloidalProxy::return_variables, PoloidalProxy(pressure_exponent, cutoff_pressure, - vector_potential_amplitude), + vector_potential_amplitude, {{0.0, 0.0, 0.0}}, 100.0), "Poloidal", {"magnetic_field"}, {{{-10.0, 10.0}}}, std::make_tuple(pressure_exponent, cutoff_pressure, vector_potential_amplitude), @@ -123,9 +142,9 @@ SPECTRE_TEST_CASE( d_pressure.get(1) = 2.0 * y; d_pressure.get(2) = 1.0; - const auto b_field = - get>(solution.variables( - in_coords, pressure, sqrt_det_spatial_metric, d_pressure)); + tnsr::I b_field{num_grid_pts, 0.0}; + solution.variables(&b_field, in_coords, pressure, + sqrt_det_spatial_metric, d_pressure); Scalar mag_b_field{num_grid_pts, 0.0}; for (size_t i = 0; i < 3; ++i) { diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Test_Toroidal.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Test_Toroidal.cpp index b3b85869ad3d..5a71e0a63488 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Test_Toroidal.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Test_Toroidal.cpp @@ -16,6 +16,7 @@ #include "NumericalAlgorithms/LinearOperators/Divergence.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp" #include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp" #include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -33,28 +34,38 @@ struct ToroidalProxy : Toroidal { const tnsr::I& x, const Scalar& pressure, const Scalar& sqrt_det_spatial_metric, const tnsr::i& deriv_pressure) const { - return this->variables(x, pressure, sqrt_det_spatial_metric, - deriv_pressure); + auto mag_field = make_with_value>(pressure, 0.0); + this->variables(make_not_null(&mag_field), x, pressure, + sqrt_det_spatial_metric, deriv_pressure); + return {mag_field}; } }; SPECTRE_TEST_CASE( "Unit.PointwiseFunctions.AnalyticData.GrMhd.InitialMagneticFields.Toroidal", "[Unit][PointwiseFunctions]") { + register_derived_classes_with_charm(); // test creation - const auto solution = TestHelpers::test_creation( - " PressureExponent: 2\n" - " CutoffPressure: 1.0e-5\n" - " VectorPotentialAmplitude: 1000.0\n"); - CHECK(solution == Toroidal(2, 1.0e-5, 1000.0)); + const auto factory_solution = serialize_and_deserialize( + TestHelpers::test_factory_creation( + "Toroidal:\n" + " PressureExponent: 2\n" + " CutoffPressure: 1.0e-5\n" + " VectorPotentialAmplitude: 1000.0\n" + " Center: [0.0, 0.0, 0.0]\n" + " MaxDistanceFromCenter: 100.0\n")); + CHECK(factory_solution->is_equal( + Toroidal(2, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0))); + + const auto& solution = dynamic_cast(*factory_solution); // test serialize test_serialization(solution); // test move { - Toroidal toroidal_field{2, 1.0e-5, 1000.0}; - Toroidal toroidal_field_copy{2, 1.0e-5, 1000.0}; + Toroidal toroidal_field{2, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0}; + Toroidal toroidal_field_copy{2, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0}; test_move_semantics(std::move(toroidal_field), toroidal_field_copy); } @@ -67,12 +78,20 @@ SPECTRE_TEST_CASE( CHECK(dynamic_cast(deserialized_base_ptr.get()) != nullptr); // test equality - const Toroidal field_original{2, 1.0e-5, 1000.0}; + const Toroidal field_original{2, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0}; const auto field = serialize_and_deserialize(field_original); - CHECK(field == Toroidal(2, 1.0e-5, 1000.0)); - CHECK(field != Toroidal(3, 1.0e-5, 1000.0)); - CHECK(field != Toroidal(2, 2.0e-5, 1000.0)); - CHECK(field != Toroidal(2, 1.0e-5, 2000.0)); + CHECK(static_cast(field).is_equal( + Toroidal(2, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Toroidal(3, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Toroidal(2, 2.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Toroidal(2, 1.0e-5, 2000.0, {{0.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Toroidal(2, 1.0e-5, 1000.0, {{1.0, 0.0, 0.0}}, 100.0))); + CHECK_FALSE(static_cast(field).is_equal( + Toroidal(2, 1.0e-5, 1000.0, {{0.0, 0.0, 0.0}}, 10.0))); // test solution implementation pypp::SetupLocalPythonEnvironment local_python_env{ @@ -86,7 +105,7 @@ SPECTRE_TEST_CASE( pypp::check_with_random_values<1>( &ToroidalProxy::return_variables, ToroidalProxy(pressure_exponent, cutoff_pressure, - vector_potential_amplitude), + vector_potential_amplitude, {{0.0, 0.0, 0.0}}, 10000.0), "Toroidal", {"magnetic_field"}, {{{-10.0, 10.0}}}, std::make_tuple(pressure_exponent, cutoff_pressure, vector_potential_amplitude), @@ -96,7 +115,8 @@ SPECTRE_TEST_CASE( // branch for p < p_cutoff pypp::check_with_random_values<1>( &ToroidalProxy::return_variables, - ToroidalProxy(pressure_exponent, 1e5, vector_potential_amplitude), + ToroidalProxy(pressure_exponent, 1e5, vector_potential_amplitude, + {{0.0, 0.0, 0.0}}, 10000.0), "Toroidal", {"magnetic_field"}, {{{-1.0, 1.0}}}, std::make_tuple(pressure_exponent, 1e5, vector_potential_amplitude), used_for_size); @@ -134,9 +154,9 @@ SPECTRE_TEST_CASE( d_pressure.get(1) = 2.0 * y; d_pressure.get(2) = 1.0; - const auto b_field = - get>(solution.variables( - in_coords, pressure, sqrt_det_spatial_metric, d_pressure)); + tnsr::I b_field{num_grid_pts, 0.0}; + solution.variables(&b_field, in_coords, pressure, + sqrt_det_spatial_metric, d_pressure); Scalar mag_b_field{num_grid_pts, 0.0}; for (size_t i = 0; i < 3; ++i) { diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedTovStar.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedTovStar.cpp index 25be91c9e7b9..f55bed905682 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedTovStar.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedTovStar.cpp @@ -5,12 +5,15 @@ #include +#include "DataStructures/DataVector.hpp" #include "Framework/TestCreation.hpp" #include "Framework/TestHelpers.hpp" #include "NumericalAlgorithms/LinearOperators/Divergence.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp" #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp" #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp" @@ -19,6 +22,7 @@ #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp" #include "Utilities/GetOutput.hpp" +#include "Utilities/MakeVector.hpp" #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" namespace grmhd::AnalyticData { @@ -40,49 +44,145 @@ void test_equality() { 1.28e-3, std::make_unique>(100.0, 2.0), TovCoordinates::Schwarzschild, - 2, - 0.04, - 2500.0}; + {}}; const auto mag_tov = serialize_and_deserialize(mag_tov_original); CHECK( mag_tov == MagnetizedTovStar( 1.28e-3, std::make_unique>(100.0, 2.0), - TovCoordinates::Schwarzschild, 2, 0.04, 2500.0)); + TovCoordinates::Schwarzschild, {})); CHECK( mag_tov != MagnetizedTovStar( 2.28e-3, std::make_unique>(100.0, 2.0), - TovCoordinates::Schwarzschild, 2, 0.04, 2500.0)); + TovCoordinates::Schwarzschild, {})); CHECK( mag_tov != MagnetizedTovStar( 1.28e-3, std::make_unique>(100.0, 2.0), - TovCoordinates::Isotropic, 2, 0.04, 2500.0)); + TovCoordinates::Isotropic, {})); CHECK( mag_tov != MagnetizedTovStar( 1.28e-3, std::make_unique>(100.0, 2.0), - TovCoordinates::Schwarzschild, 3, 0.04, 2500.0)); + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)))); CHECK( - mag_tov != MagnetizedTovStar( 1.28e-3, std::make_unique>(100.0, 2.0), - TovCoordinates::Schwarzschild, 2, 0.05, 2500.0)); + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))) == + MagnetizedTovStar( + 1.28e-3, + std::make_unique>(100.0, 2.0), + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)))); CHECK( - mag_tov != MagnetizedTovStar( 1.28e-3, std::make_unique>(100.0, 2.0), - TovCoordinates::Schwarzschild, 2, 0.04, 3500.0)); + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Toroidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))) != + MagnetizedTovStar( + 1.28e-3, + std::make_unique>(100.0, 2.0), + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)))); + CHECK( + MagnetizedTovStar( + 1.28e-3, + std::make_unique>(100.0, 2.0), + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))) != + MagnetizedTovStar( + 1.28e-3, + std::make_unique>(100.0, 2.0), + TovCoordinates::Schwarzschild, + make_vector< + std::unique_ptr>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0), + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Toroidal>( + 2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)))); + // Check order of magnetic fields doesn't matter. + const MagnetizedTovStar toroidal_plus_poloidal( + 1.28e-3, + std::make_unique>(100.0, 2.0), + TovCoordinates::Schwarzschild, + make_vector>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Toroidal>( + 2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0), + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))); + const MagnetizedTovStar poloidal_plus_toroidal( + 1.28e-3, + std::make_unique>(100.0, 2.0), + TovCoordinates::Schwarzschild, + make_vector>( + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Poloidal>( + 2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0), + std::make_unique< + grmhd::AnalyticData::InitialMagneticFields::Toroidal>( + 2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))); + const DataVector coord{2.0}; + const tnsr::I coords{{coord, coord, coord}}; + CHECK_ITERABLE_APPROX( + (get>( + poloidal_plus_toroidal.variables( + coords, + tmpl::list>{}))), + (get>( + toroidal_plus_poloidal.variables( + coords, + tmpl::list>{})))); } void test_magnetized_tov_star(const TovCoordinates coord_system) { + register_derived_classes_with_charm< + grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>(); register_classes_with_charm(); register_classes_with_charm>(); const std::unique_ptr option_solution = @@ -98,9 +198,13 @@ void test_magnetized_tov_star(const TovCoordinates coord_system) { " Coordinates: " + get_output(coord_system) + "\n" - " PressureExponent: 2\n" - " VectorPotentialAmplitude: 2500\n" - " CutoffPressureFraction: 0.04\n") + " MagneticFields:\n" + " - Poloidal:\n" + " PressureExponent: 2\n" + " VectorPotentialAmplitude: 2500\n" + " CutoffPressure: 6.5536e-06\n" // 0.04 * 100 * (1.28e-3)**2 + " Center: [0.0, 0.0, 0.0]\n" + " MaxDistanceFromCenter: 100.0\n") ->get_clone(); const auto deserialized_option_solution = serialize_and_deserialize(option_solution);