Add m1 analytic data & couple GreyM1 to fixed hydro

docs/References.bib

@@ -2278,6 +2278,29 @@ @article{Porth2016rfi
   SLACcitation   = "%%CITATION = ARXIV:1611.09720;%%"
 }
 
+@ARTICLE{radice2022,
+       author = {{Radice}, David and {Bernuzzi}, Sebastiano and {Perego},
+       Albino and {Haas}, Roland},
+        title = "{A new moment-based general-relativistic neutrino-radiation
+        transport code: Methods and first applications to neutron star
+        mergers}",
+      journal = {\mnras},
+     keywords = {neutrinos, methods: numerical, neutron star mergers,
+     Astrophysics - High Energy Astrophysical Phenomena, General Relativity
+     and Quantum Cosmology},
+         year = 2022,
+        month = may,
+       volume = {512},
+       number = {1},
+        pages = {1499-1521},
+          doi = {10.1093/mnras/stac589},
+archivePrefix = {arXiv},
+       eprint = {2111.14858},
+ primaryClass = {astro-ph.HE},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2022MNRAS.512.1499R},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
 @article{Renkhoff2023,
   author        = {Renkhoff, Sarah and Cors, Daniela and Hilditch, David and
                    Br\"ugmann, Bernd},

src/PointwiseFunctions/AnalyticData/RadiationTransport/M1Grey/HomogeneousSphere.cpp

@@ -0,0 +1,226 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "PointwiseFunctions/AnalyticData/RadiationTransport/M1Grey/HomogeneousSphere.hpp"
+
+#include <boost/preprocessor/punctuation/comma_if.hpp>
+#include <boost/preprocessor/repetition/repeat.hpp>
+#include <cmath>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Evolution/Systems/RadiationTransport/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/MakeWithValue.hpp"
+#include "Utilities/Math.hpp"
+
+namespace RadiationTransport::M1Grey::AnalyticData {
+HomogeneousSphere::HomogeneousSphere(const double radius,
+                                     const double emissivity_and_opacity,
+                                     const double outer_radius,
+                                     const double outer_opacity)
+    : radius_(radius),
+      emissivity_and_opacity_(emissivity_and_opacity),
+      outer_radius_(outer_radius),
+      outer_opacity_(outer_opacity) {
+  if (UNLIKELY(radius_ > outer_radius_)) {
+    ERROR("Radius " << radius_
+                    << " is greater than the  "
+                       "outer radius: "
+                    << outer_radius_);
+  }
+}
+
+namespace {
+// This function is used to round the edges of the homogeneous sphere, as
+// opposed to a pure, rectangular step function.
+Scalar<DataVector> rounded_step_function(const DataVector& x,
+                                         const double inner_value,
+                                         const double outer_value,
+                                         const double sphere_radius) {
+  // the closer to 0 this becomes, the sharper the discontinuity
+  const double sharpness = -0.03;
+  return Scalar<DataVector>{
+      (inner_value - outer_value) / M_PI *
+          atan((sqrt(square(x)) - sphere_radius) / sharpness) +
+      0.5 * (inner_value + outer_value)};
+}
+}  // namespace
+
+template <typename NeutrinoSpecies>
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<RadiationTransport::M1Grey::Tags::TildeE<
+        Frame::Inertial, NeutrinoSpecies>> /*meta*/) const
+    -> tuples::TaggedTuple<RadiationTransport::M1Grey::Tags::TildeE<
+        Frame::Inertial, NeutrinoSpecies>> {
+  const DataVector r = sqrt(radius_squared(x));
+
+  const double inner_energy = 1.0;
+  const double outer_energy = 1.0e-12;
+  const double sphere_radius = radius_;
+
+  return rounded_step_function(r, inner_energy, outer_energy, sphere_radius);
+}
+
+template <typename NeutrinoSpecies>
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<RadiationTransport::M1Grey::Tags::TildeS<
+        Frame::Inertial, NeutrinoSpecies>> /*meta*/) const
+    -> tuples::TaggedTuple<RadiationTransport::M1Grey::Tags::TildeS<
+        Frame::Inertial, NeutrinoSpecies>> {
+  return {make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(x, 0.0)};
+}
+
+template <typename NeutrinoSpecies>
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<RadiationTransport::M1Grey::Tags::GreyEmissivity<
+        NeutrinoSpecies>> /*meta*/) const
+    -> tuples::TaggedTuple<
+        RadiationTransport::M1Grey::Tags::GreyEmissivity<NeutrinoSpecies>> {
+  const DataVector r = sqrt(radius_squared(x));
+
+  const double inner_emissivity = emissivity_and_opacity_;
+  const double outer_emissivity = 0.0;
+  const double sphere_radius = radius_;
+
+  return rounded_step_function(r, inner_emissivity, outer_emissivity,
+                               sphere_radius);
+}
+
+template <typename NeutrinoSpecies>
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<RadiationTransport::M1Grey::Tags::GreyAbsorptionOpacity<
+        NeutrinoSpecies>> /*meta*/) const
+    -> tuples::TaggedTuple<RadiationTransport::M1Grey::Tags::
+                               GreyAbsorptionOpacity<NeutrinoSpecies>> {
+  const DataVector r = sqrt(radius_squared(x));
+
+  const double sphere_radius = radius_;
+
+  return rounded_step_function(r, emissivity_and_opacity_, outer_opacity_,
+                               sphere_radius);
+}
+
+template <typename NeutrinoSpecies>
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<RadiationTransport::M1Grey::Tags::GreyScatteringOpacity<
+        NeutrinoSpecies>> /*meta*/) const
+    -> tuples::TaggedTuple<RadiationTransport::M1Grey::Tags::
+                               GreyScatteringOpacity<NeutrinoSpecies>> {
+  return {make_with_value<Scalar<DataVector>>(x, 0.0)};
+}
+
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<hydro::Tags::LorentzFactor<DataVector>> /*meta*/)
+    -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataVector>> {
+  return {make_with_value<Scalar<DataVector>>(x, 1.0)};
+}
+
+auto HomogeneousSphere::variables(
+    const tnsr::I<DataVector, 3>& x,
+    tmpl::list<hydro::Tags::SpatialVelocity<DataVector, 3>> /*meta*/)
+    -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataVector, 3>> {
+  return {make_with_value<tnsr::I<DataVector, 3, Frame::Inertial>>(x, 0.0)};
+}
+
+std::unique_ptr<evolution::initial_data::InitialData>
+HomogeneousSphere::get_clone() const {
+  return std::make_unique<HomogeneousSphere>(*this);
+}
+
+void HomogeneousSphere::pup(PUP::er& p) {
+  evolution::initial_data::InitialData::pup(p);
+  p | radius_;
+  p | emissivity_and_opacity_;
+  p | outer_radius_;
+  p | outer_opacity_;
+}
+// NOLINTNEXTLINE
+PUP::able::PUP_ID HomogeneousSphere::my_PUP_ID = 0;
+
+bool operator!=(const HomogeneousSphere& lhs, const HomogeneousSphere& rhs) {
+  return not(lhs == rhs);
+}
+
+DataVector HomogeneousSphere::radius_squared(const tnsr::I<DataVector, 3>& x) {
+  return square(get<0>(x)) + square(get<1>(x)) + square(get<2>(x));
+}
+
+#define DERIVED_CLASSES (HomogeneousSphere)
+
+#define DERIVED(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define TAG(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define NTYPE(data) BOOST_PP_TUPLE_ELEM(2, data)
+#define EBIN(data) BOOST_PP_TUPLE_ELEM(3, data)
+#define GENERATE_LIST(z, n, _) BOOST_PP_COMMA_IF(n) n
+
+#define INSTANTIATE_M1_FUNCTION_WITH_FRAME(_, data)                            \
+  template tuples::TaggedTuple<TAG(data) < Frame::Inertial,                    \
+                               NTYPE(data) < EBIN(data)> >>                    \
+      DERIVED(data)::variables(                                                \
+          const tnsr::I<DataVector, 3>& x,                                     \
+          tmpl::list<TAG(data) < Frame::Inertial, NTYPE(data) < EBIN(data)> >> \
+          /*meta*/) const;
+
+#define temp_list \
+  (BOOST_PP_REPEAT(MAX_NUMBER_OF_NEUTRINO_ENERGY_BINS, GENERATE_LIST, _))
+
+GENERATE_INSTANTIATIONS(INSTANTIATE_M1_FUNCTION_WITH_FRAME, DERIVED_CLASSES,
+                        (RadiationTransport::M1Grey::Tags::TildeE,
+                         RadiationTransport::M1Grey::Tags::TildeS),
+                        (neutrinos::ElectronNeutrinos,
+                         neutrinos::ElectronAntiNeutrinos,
+                         neutrinos::HeavyLeptonNeutrinos),
+                        temp_list)
+
+#undef temp_list
+#undef INSTANTIATE_M1_FUNCTION_WITH_FRAME
+#undef DERIVED
+#undef TAG
+#undef NTYPE
+#undef EBIN
+#undef GENERATE_LIST
+
+#define DERIVED(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define TAG(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define NTYPE(data) BOOST_PP_TUPLE_ELEM(2, data)
+#define EBIN(data) BOOST_PP_TUPLE_ELEM(3, data)
+#define GENERATE_LIST(z, n, _) BOOST_PP_COMMA_IF(n) n
+
+#define INSTANTIATE_M1_FUNCTION(_, data)                                \
+  template tuples::TaggedTuple<TAG(data) < NTYPE(data) < EBIN(data)> >> \
+      DERIVED(data)::variables(                                         \
+          const tnsr::I<DataVector, 3>& x,                              \
+          tmpl::list<TAG(data) < NTYPE(data) < EBIN(data)> >>           \
+          /*meta*/) const;
+
+#define temp_list \
+  (BOOST_PP_REPEAT(MAX_NUMBER_OF_NEUTRINO_ENERGY_BINS, GENERATE_LIST, _))
+
+GENERATE_INSTANTIATIONS(
+    INSTANTIATE_M1_FUNCTION, DERIVED_CLASSES,
+    (RadiationTransport::M1Grey::Tags::GreyEmissivity,
+     RadiationTransport::M1Grey::Tags::GreyAbsorptionOpacity,
+     RadiationTransport::M1Grey::Tags::GreyScatteringOpacity),
+    (neutrinos::ElectronNeutrinos, neutrinos::ElectronAntiNeutrinos,
+     neutrinos::HeavyLeptonNeutrinos),
+    temp_list)
+
+#undef INSTANTIATE_M1_FUNCTION
+#undef temp_list
+#undef DERIVED
+#undef TAG
+#undef NTYPE
+#undef EBIN
+#undef GENERATE_LIST
+
+// template class HomogeneousSphere;
+
+}  // namespace RadiationTransport::M1Grey::AnalyticData