Merge pull request #1074 from carmaza/bondi-hoyle-initial-data

Add Bondi-Hoyle Accretion initial data

src/PointwiseFunctions/AnalyticData/CMakeLists.txt

@@ -0,0 +1,4 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+add_subdirectory(GrMhd)

src/PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp

@@ -0,0 +1,12 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+namespace grmhd {
+/*!
+ * \ingroup AnalyticDataGroup
+ * \brief Holds classes implementing analytic data for the GrMhd system.
+ */
+namespace AnalyticData {}
+}  // namespace grmhd

src/PointwiseFunctions/AnalyticData/GrMhd/BondiHoyleAccretion.cpp

@@ -0,0 +1,255 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "PointwiseFunctions/AnalyticData/GrMhd/BondiHoyleAccretion.hpp"
+
+#include <cmath>
+#include <pup.h>
+
+#include "DataStructures/DataBox/DataBoxTag.hpp"
+#include "DataStructures/DataVector.hpp"  // IWYU pragma: keep
+#include "DataStructures/Tensor/Tensor.hpp"  // IWYU pragma: keep
+#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"  // IWYU pragma: keep
+#include "Utilities/ConstantExpressions.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/MakeWithValue.hpp"
+
+// IWYU pragma:  no_include "DataStructures/Tensor/TypeAliases.hpp"
+
+/// \cond
+namespace grmhd {
+namespace AnalyticData {
+
+BondiHoyleAccretion::BondiHoyleAccretion(
+    const double bh_mass, const double bh_dimless_spin,
+    const double rest_mass_density, const double flow_speed,
+    const double magnetic_field_strength, const double polytropic_constant,
+    const double polytropic_exponent) noexcept
+    : bh_mass_(bh_mass),
+      bh_spin_a_(bh_mass * bh_dimless_spin),
+      rest_mass_density_(rest_mass_density),
+      flow_speed_(flow_speed),
+      magnetic_field_strength_(magnetic_field_strength),
+      polytropic_constant_(polytropic_constant),
+      polytropic_exponent_(polytropic_exponent),
+      equation_of_state_{polytropic_constant_, polytropic_exponent_},
+      background_spacetime_{
+          bh_mass, {{0.0, 0.0, bh_dimless_spin}}, {{0.0, 0.0, 0.0}}},
+      kerr_schild_coords_{bh_mass, bh_dimless_spin} {}
+
+void BondiHoyleAccretion::pup(PUP::er& p) noexcept {
+  p | bh_mass_;
+  p | bh_spin_a_;
+  p | rest_mass_density_;
+  p | flow_speed_;
+  p | magnetic_field_strength_;
+  p | polytropic_constant_;
+  p | polytropic_exponent_;
+  p | equation_of_state_;
+  p | background_spacetime_;
+  p | kerr_schild_coords_;
+}
+
+template <typename DataType>
+typename hydro::Tags::SpatialVelocity<DataType, 3, Frame::NoFrame>::type
+BondiHoyleAccretion::spatial_velocity(const DataType& r_squared,
+                                      const DataType& cos_theta,
+                                      const DataType& sin_theta) const
+    noexcept {
+  auto result = make_with_value<
+      typename hydro::Tags::SpatialVelocity<DataType, 3, Frame::NoFrame>::type>(
+      r_squared, 0.0);
+
+  const DataType sigma = r_squared + square(bh_spin_a_) * square(cos_theta);
+  get<0>(result) = flow_speed_ * cos_theta /
+                   sqrt(1.0 + 2.0 * bh_mass_ * sqrt(r_squared) / sigma);
+  get<1>(result) = -flow_speed_ * sin_theta / sqrt(sigma);
+  // get<2>(result) is identically zero
+
+  return result;
+}
+
+template <typename DataType>
+typename hydro::Tags::MagneticField<DataType, 3, Frame::NoFrame>::type
+BondiHoyleAccretion::magnetic_field(const DataType& r_squared,
+                                    const DataType& cos_theta,
+                                    const DataType& sin_theta) const noexcept {
+  const double a_squared = bh_spin_a_ * bh_spin_a_;
+  const DataType cos_theta_squared = square(cos_theta);
+  const DataType two_m_r = 2.0 * bh_mass_ * sqrt(r_squared);
+
+  // The square root of the determinant of the spatial metric is proportional to
+  // sin(theta) squared. That factor cancels out after multiplying by Faraday.
+  DataType sigma = r_squared + a_squared * cos_theta_squared;
+  const DataType prefactor =
+      magnetic_field_strength_ / sqrt(sigma * (sigma + two_m_r));
+
+  auto result = make_with_value<
+      typename hydro::Tags::MagneticField<DataType, 3, Frame::NoFrame>::type>(
+      r_squared, 0.0);
+
+  // Redefine sigma to save an allocation.
+  sigma = 1.0 / square(sigma);
+  get<0>(result) = prefactor *
+                   (r_squared - two_m_r + a_squared +
+                    sigma * two_m_r * (square(r_squared) - square(a_squared))) *
+                   cos_theta;
+
+  get<1>(result) = -prefactor *
+                   (sqrt(r_squared) +
+                    bh_mass_ * a_squared * sigma *
+                        (r_squared - a_squared * cos_theta_squared) *
+                        (1.0 + cos_theta_squared)) *
+                   sin_theta;
+
+  get<2>(result) = bh_spin_a_ * prefactor *
+                   (1.0 + sigma * two_m_r * (r_squared - a_squared)) *
+                   cos_theta;
+
+  return result;
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/) const
+    noexcept {
+  return {
+      make_with_value<db::item_type<hydro::Tags::RestMassDensity<DataType>>>(
+          x, rest_mass_density_)};
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/) const
+    noexcept {
+  const DataType r_squared = get(kerr_schild_coords_.r_coord_squared(x));
+  return kerr_schild_coords_.cartesian_from_spherical_ks(
+      spatial_velocity(r_squared, DataType{get<2>(x) / sqrt(r_squared)},
+                       DataType{sqrt((square(get<0>(x)) + square(get<1>(x))) /
+                                     (r_squared + square(bh_spin_a_)))}),
+      x);
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
+    noexcept {
+  return equation_of_state_.specific_internal_energy_from_density(
+      get<hydro::Tags::RestMassDensity<DataType>>(
+          variables(x, tmpl::list<hydro::Tags::RestMassDensity<DataType>>{})));
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const noexcept {
+  return equation_of_state_.pressure_from_density(
+      get<hydro::Tags::RestMassDensity<DataType>>(
+          variables(x, tmpl::list<hydro::Tags::RestMassDensity<DataType>>{})));
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/) const
+    noexcept {
+  const DataType r_squared = get(kerr_schild_coords_.r_coord_squared(x));
+  return kerr_schild_coords_.cartesian_from_spherical_ks(
+      magnetic_field(r_squared, DataType{get<2>(x) / sqrt(r_squared)},
+                     DataType{sqrt((square(get<0>(x)) + square(get<1>(x))) /
+                                   (r_squared + square(bh_spin_a_)))}),
+      x);
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
+    noexcept {
+  return {make_with_value<
+      db::item_type<hydro::Tags::DivergenceCleaningField<DataType>>>(x, 0.0)};
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/) const noexcept {
+  return {make_with_value<db::item_type<hydro::Tags::LorentzFactor<DataType>>>(
+      x, 1.0 / sqrt(1.0 - square(flow_speed_)))};
+}
+
+template <typename DataType>
+tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>
+BondiHoyleAccretion::variables(
+    const tnsr::I<DataType, 3, Frame::Inertial>& x,
+    tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/) const
+    noexcept {
+  return equation_of_state_.specific_enthalpy_from_density(
+      get<hydro::Tags::RestMassDensity<DataType>>(
+          variables(x, tmpl::list<hydro::Tags::RestMassDensity<DataType>>{})));
+}
+
+bool operator==(const BondiHoyleAccretion& lhs,
+                const BondiHoyleAccretion& rhs) noexcept {
+  // There is no comparison operator for the `equation_of_state` and the
+  // `background_spacetime`, but should be okay as the `bh_mass`s,
+  // `bh_dimless_spin`s, `polytropic_exponent`s and `polytropic_constant`s are
+  // compared.
+  return lhs.bh_mass_ == rhs.bh_mass_ and lhs.bh_spin_a_ == rhs.bh_spin_a_ and
+         lhs.rest_mass_density_ == rhs.rest_mass_density_ and
+         lhs.flow_speed_ == rhs.flow_speed_ and
+         lhs.magnetic_field_strength_ == rhs.magnetic_field_strength_ and
+         lhs.polytropic_constant_ == rhs.polytropic_constant_ and
+         lhs.polytropic_exponent_ == rhs.polytropic_exponent_;
+}
+
+bool operator!=(const BondiHoyleAccretion& lhs,
+                const BondiHoyleAccretion& rhs) noexcept {
+  return not(lhs == rhs);
+}
+
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define TAG(data) BOOST_PP_TUPLE_ELEM(1, data)
+
+#define INSTANTIATE_SCALARS(_, data)                         \
+  template tuples::TaggedTuple<TAG(data) < DTYPE(data)>>     \
+      BondiHoyleAccretion::variables(                        \
+          const tnsr::I<DTYPE(data), 3, Frame::Inertial>& x, \
+          tmpl::list<TAG(data) < DTYPE(data)>> /*meta*/) const noexcept;
+
+GENERATE_INSTANTIATIONS(
+    INSTANTIATE_SCALARS, (double, DataVector),
+    (hydro::Tags::RestMassDensity, hydro::Tags::SpecificInternalEnergy,
+     hydro::Tags::Pressure, hydro::Tags::DivergenceCleaningField,
+     hydro::Tags::LorentzFactor, hydro::Tags::SpecificEnthalpy))
+
+#define INSTANTIATE_VECTORS(_, data)                                         \
+  template tuples::TaggedTuple<TAG(data) < DTYPE(data), 3>>                  \
+      BondiHoyleAccretion::variables(                                        \
+          const tnsr::I<DTYPE(data), 3, Frame::Inertial>& x,                 \
+          tmpl::list<TAG(data) < DTYPE(data), 3, Frame::Inertial>> /*meta*/) \
+          const noexcept;
+
+GENERATE_INSTANTIATIONS(INSTANTIATE_VECTORS, (double, DataVector),
+                        (hydro::Tags::SpatialVelocity,
+                         hydro::Tags::MagneticField))
+
+#undef DTYPE
+#undef TAG
+#undef INSTANTIATE_SCALARS
+#undef INSTANTIATE_VECTORS
+
+}  // namespace AnalyticData
+}  // namespace grmhd
+/// \endcond