4175800929

Build Type

push

github

Committed by GitHub

Commit Message

Merge pull request #4727 from knelli2/fix_shape_naming

Run Details

5 of 5 new or added lines in 1 file covered. (100.0%)

63479 of 66187 relevant lines covered (95.91%)

365273.36 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.36

/src/PointwiseFunctions/AnalyticSolutions/Xcts/ConstantDensityStar.cpp

// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "PointwiseFunctions/AnalyticSolutions/Xcts/ConstantDensityStar.hpp"

#include <array>
#include <cmath>
#include <cstddef>
#include <ostream>
#include <pup.h>
#include <utility>

#include "DataStructures/DataVector.hpp"
#include "DataStructures/Tensor/EagerMath/Magnitude.hpp"  // IWYU pragma: keep
#include "DataStructures/Tensor/Tensor.hpp"               // IWYU pragma: keep
#include "Elliptic/Systems/Xcts/Tags.hpp"
#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
#include "Options/Options.hpp"
#include "Utilities/ConstantExpressions.hpp"
#include "Utilities/ContainerHelpers.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/MakeWithValue.hpp"

// IWYU pragma: no_forward_declare Tensor

namespace {

// Find the alpha parameter that corresponds to the weak-field solution,
// since this is the solution we get when we set \psi = 1 initially
double compute_alpha(const double density, const double radius) {
  const double alpha_source = sqrt(2. * M_PI * density / 3.) * radius;
  return RootFinder::toms748(
      [alpha_source](const double a) {
        const double a_square = pow<2>(a);
        const double pow_2_one_plus_a_square = pow<2>(1. + a_square);
        return alpha_source - a_square * pow<3>(a) /
                                  (pow_2_one_plus_a_square * (1. + a_square));
      },
      sqrt(5.), 1.0 / alpha_source,
      // Choose a precision of 14 base-10 digits for no particular reason
      1.0e-14, 1.0e-15);
}

template <typename DataType>
Scalar<DataType> compute_piecewise(const Scalar<DataType>& r, double radius,
                                   double inner_value, double outer_value);
template <>
Scalar<double> compute_piecewise(const Scalar<double>& r, const double radius,
                                 const double inner_value,
                                 const double outer_value) {
  return Scalar<double>(get(r) < radius ? inner_value : outer_value);
}
template <>
Scalar<DataVector> compute_piecewise(const Scalar<DataVector>& r,
                                     const double radius,
                                     const double inner_value,
                                     const double outer_value) {
  return Scalar<DataVector>(inner_value - (inner_value - outer_value) *
                                              step_function(get(r) - radius));
}

}  // namespace

namespace Xcts::Solutions {

ConstantDensityStar::ConstantDensityStar(const double density,
                                         const double radius,
                                         const Options::Context& context)
    : density_(density), radius_(radius) {
  const double critical_density =
      3. * pow<5>(5.) / (2. * pow<6>(6.) * M_PI * pow<2>(radius));
  if (density <= critical_density) {
    alpha_ = compute_alpha(density, radius);
  } else {
    PARSE_ERROR(context,
                "A ConstantDensityStar has no solutions for a density below "
                "the critical density ("
                    << critical_density << ").");
  }
}

void ConstantDensityStar::pup(PUP::er& p) {
  elliptic::analytic_data::AnalyticSolution::pup(p);
  p | density_;
  p | radius_;
  if (p.isUnpacking()) {
    alpha_ = compute_alpha(density_, radius_);
  }
}

template <typename DataType>
tuples::TaggedTuple<Xcts::Tags::ConformalFactor<DataType>>
ConstantDensityStar::variables(
    const tnsr::I<DataType, 3, Frame::Inertial>& x,
    tmpl::list<Xcts::Tags::ConformalFactor<DataType>> /*meta*/) const {
  const DataType r = get(magnitude(x));
  const double inner_prefactor =
      sqrt(alpha_ * radius_) / std::pow(2. * M_PI * density_ / 3., 0.25);
  const double alpha_times_radius_square = square(alpha_ * radius_);
  const double beta = inner_prefactor / sqrt(1. + square(alpha_)) - radius_;
  auto conformal_factor = make_with_value<Scalar<DataType>>(r, 0.);
  for (size_t i = 0; i < get_size(r); i++) {
    if (get_element(r, i) <= radius_) {
      get_element(get(conformal_factor), i) =
          inner_prefactor /
          sqrt(square(get_element(r, i)) + alpha_times_radius_square);
    } else {
      get_element(get(conformal_factor), i) = beta / get_element(r, i) + 1.;
    }
  }
  return {std::move(conformal_factor)};
}

template <typename DataType>
tuples::TaggedTuple<::Tags::Initial<Xcts::Tags::ConformalFactor<DataType>>>
ConstantDensityStar::variables(
    const tnsr::I<DataType, 3, Frame::Inertial>& x,
    tmpl::list<::Tags::Initial<Xcts::Tags::ConformalFactor<DataType>>> /*meta*/)
    const {
  return {make_with_value<Scalar<DataType>>(x, 1.)};
}

template <typename DataType>
tuples::TaggedTuple<::Tags::Initial<::Tags::deriv<
    Xcts::Tags::ConformalFactor<DataType>, tmpl::size_t<3>, Frame::Inertial>>>
ConstantDensityStar::variables(
    const tnsr::I<DataType, 3, Frame::Inertial>& x,
    tmpl::list<::Tags::Initial<
        ::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>, tmpl::size_t<3>,
                      Frame::Inertial>>> /*meta*/) const {
  return {make_with_value<tnsr::i<DataType, 3, Frame::Inertial>>(x, 0.)};
}

template <typename DataType>
tuples::TaggedTuple<::Tags::FixedSource<Xcts::Tags::ConformalFactor<DataType>>>
ConstantDensityStar::variables(
    const tnsr::I<DataType, 3, Frame::Inertial>& x,
    tmpl::list<
        ::Tags::FixedSource<Xcts::Tags::ConformalFactor<DataType>>> /*meta*/)
    const {
  return {make_with_value<Scalar<DataType>>(x, 0.)};
}

template <typename DataType>
tuples::TaggedTuple<gr::Tags::EnergyDensity<DataType>>
ConstantDensityStar::variables(
    const tnsr::I<DataType, 3, Frame::Inertial>& x,
    tmpl::list<gr::Tags::EnergyDensity<DataType>> /*meta*/) const {
  return {compute_piecewise(magnitude(x), radius_, density_, 0.)};
}

PUP::able::PUP_ID ConstantDensityStar::my_PUP_ID = 0;  // NOLINT

bool operator==(const ConstantDensityStar& lhs,
                const ConstantDensityStar& rhs) {
  return lhs.density() == rhs.density() and lhs.radius() == rhs.radius();
}

bool operator!=(const ConstantDensityStar& lhs,
                const ConstantDensityStar& rhs) {
  return not(lhs == rhs);
}

#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)

#define INSTANTIATE(_, data)                                                  \
  template tuples::TaggedTuple<Xcts::Tags::ConformalFactor<DTYPE(data)>>      \
  ConstantDensityStar::variables(                                             \
      const tnsr::I<DTYPE(data), 3, Frame::Inertial>&,                        \
      tmpl::list<Xcts::Tags::ConformalFactor<DTYPE(data)>>) const;            \
  template tuples::TaggedTuple<                                               \
      ::Tags::Initial<Xcts::Tags::ConformalFactor<DTYPE(data)>>>              \
  ConstantDensityStar::variables(                                             \
      const tnsr::I<DTYPE(data), 3, Frame::Inertial>&,                        \
      tmpl::list<::Tags::Initial<Xcts::Tags::ConformalFactor<DTYPE(data)>>>)  \
      const;                                                                  \
  template tuples::TaggedTuple<                                               \
      ::Tags::Initial<::Tags::deriv<Xcts::Tags::ConformalFactor<DTYPE(data)>, \
                                    tmpl::size_t<3>, Frame::Inertial>>>       \
  ConstantDensityStar::variables(                                             \
      const tnsr::I<DTYPE(data), 3, Frame::Inertial>&,                        \
      tmpl::list<::Tags::Initial<                                             \
          ::Tags::deriv<Xcts::Tags::ConformalFactor<DTYPE(data)>,             \
                        tmpl::size_t<3>, Frame::Inertial>>>) const;           \
  template tuples::TaggedTuple<                                               \
      ::Tags::FixedSource<Xcts::Tags::ConformalFactor<DTYPE(data)>>>          \
  ConstantDensityStar::variables(                                             \
      const tnsr::I<DTYPE(data), 3, Frame::Inertial>&,                        \
      tmpl::list<                                                             \
          ::Tags::FixedSource<Xcts::Tags::ConformalFactor<DTYPE(data)>>>)     \
      const;                                                                  \
  template tuples::TaggedTuple<gr::Tags::EnergyDensity<DTYPE(data)>>          \
  ConstantDensityStar::variables(                                             \
      const tnsr::I<DTYPE(data), 3, Frame::Inertial>&,                        \
      tmpl::list<gr::Tags::EnergyDensity<DTYPE(data)>>) const;

GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))

#undef DTYPE
#undef INSTANTIATE

}  // namespace Xcts::Solutions

1	// Distributed under the MIT License.
2	// See LICENSE.txt for details.
3
4	#include "PointwiseFunctions/AnalyticSolutions/Xcts/ConstantDensityStar.hpp"
5
6	#include <array>
7	#include <cmath>
8	#include <cstddef>
9	#include <ostream>
10	#include <pup.h>
11	#include <utility>
12
13	#include "DataStructures/DataVector.hpp"
14	#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" // IWYU pragma: keep
15	#include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep
16	#include "Elliptic/Systems/Xcts/Tags.hpp"
17	#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
18	#include "Options/Options.hpp"
19	#include "Utilities/ConstantExpressions.hpp"
20	#include "Utilities/ContainerHelpers.hpp"
21	#include "Utilities/GenerateInstantiations.hpp"
22	#include "Utilities/MakeWithValue.hpp"
23
24	// IWYU pragma: no_forward_declare Tensor
25
26	namespace {
27
28	// Find the alpha parameter that corresponds to the weak-field solution,
29	// since this is the solution we get when we set \psi = 1 initially
30	double compute_alpha(const double density, const double radius) {	3✔
31	const double alpha_source = sqrt(2. * M_PI * density / 3.) * radius;	3✔
32	return RootFinder::toms748(	6✔
33	[alpha_source](const double a) {	27✔
34	const double a_square = pow<2>(a);	27✔
35	const double pow_2_one_plus_a_square = pow<2>(1. + a_square);	54✔
36	return alpha_source - a_square * pow<3>(a) /	27✔
37	(pow_2_one_plus_a_square * (1. + a_square));	27✔
38	},
39	sqrt(5.), 1.0 / alpha_source,
40	// Choose a precision of 14 base-10 digits for no particular reason
41	1.0e-14, 1.0e-15);	6✔
42	}
43
44	template <typename DataType>
45	Scalar<DataType> compute_piecewise(const Scalar<DataType>& r, double radius,
46	double inner_value, double outer_value);
47	template <>
48	Scalar<double> compute_piecewise(const Scalar<double>& r, const double radius,	×
49	const double inner_value,
50	const double outer_value) {
51	return Scalar<double>(get(r) < radius ? inner_value : outer_value);	×
52	}
53	template <>
54	Scalar<DataVector> compute_piecewise(const Scalar<DataVector>& r,	×
55	const double radius,
56	const double inner_value,
57	const double outer_value) {
58	return Scalar<DataVector>(inner_value - (inner_value - outer_value) *	×
59	step_function(get(r) - radius));	×
60	}
61
62	} // namespace
63
64	namespace Xcts::Solutions {
65
66	ConstantDensityStar::ConstantDensityStar(const double density,	2✔
67	const double radius,
68	const Options::Context& context)	1✔
69	: density_(density), radius_(radius) {	2✔
70	const double critical_density =
71	3. * pow<5>(5.) / (2. * pow<6>(6.) * M_PI * pow<2>(radius));	6✔
72	if (density <= critical_density) {	2✔
73	alpha_ = compute_alpha(density, radius);	2✔
74	} else {
75	PARSE_ERROR(context,	×
76	"A ConstantDensityStar has no solutions for a density below "
77	"the critical density ("
78	<< critical_density << ").");
79	}
80	}	2✔
81
82	void ConstantDensityStar::pup(PUP::er& p) {	3✔
83	elliptic::analytic_data::AnalyticSolution::pup(p);	3✔
84	p \| density_;	3✔
85	p \| radius_;	3✔
86	if (p.isUnpacking()) {	3✔
87	alpha_ = compute_alpha(density_, radius_);	1✔
88	}
89	}	3✔
90
91	template <typename DataType>
92	tuples::TaggedTuple<Xcts::Tags::ConformalFactor<DataType>>
93	ConstantDensityStar::variables(	2✔
94	const tnsr::I<DataType, 3, Frame::Inertial>& x,
95	tmpl::list<Xcts::Tags::ConformalFactor<DataType>> /meta/) const {
96	const DataType r = get(magnitude(x));	6✔
97	const double inner_prefactor =	2✔
98	sqrt(alpha_ * radius_) / std::pow(2. * M_PI * density_ / 3., 0.25);	2✔
99	const double alpha_times_radius_square = square(alpha_ * radius_);	2✔
100	const double beta = inner_prefactor / sqrt(1. + square(alpha_)) - radius_;	2✔
101	auto conformal_factor = make_with_value<Scalar<DataType>>(r, 0.);	4✔
102	for (size_t i = 0; i < get_size(r); i++) {	20✔
103	if (get_element(r, i) <= radius_) {	8✔
104	get_element(get(conformal_factor), i) =	×
105	inner_prefactor /	×
106	sqrt(square(get_element(r, i)) + alpha_times_radius_square);	×
107	} else {
108	get_element(get(conformal_factor), i) = beta / get_element(r, i) + 1.;	24✔
109	}
110	}
111	return {std::move(conformal_factor)};	4✔
112	}
113
114	template <typename DataType>
115	tuples::TaggedTuple<::Tags::Initial<Xcts::Tags::ConformalFactor<DataType>>>
116	ConstantDensityStar::variables(	2✔
117	const tnsr::I<DataType, 3, Frame::Inertial>& x,
118	tmpl::list<::Tags::Initial<Xcts::Tags::ConformalFactor<DataType>>> /meta/)
119	const {
120	return {make_with_value<Scalar<DataType>>(x, 1.)};	4✔
121	}
122
123	template <typename DataType>
124	tuples::TaggedTuple<::Tags::Initial<::Tags::deriv<
125	Xcts::Tags::ConformalFactor<DataType>, tmpl::size_t<3>, Frame::Inertial>>>
126	ConstantDensityStar::variables(	2✔
127	const tnsr::I<DataType, 3, Frame::Inertial>& x,
128	tmpl::list<::Tags::Initial<
129	::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>, tmpl::size_t<3>,
130	Frame::Inertial>>> /meta/) const {
131	return {make_with_value<tnsr::i<DataType, 3, Frame::Inertial>>(x, 0.)};	4✔
132	}
133
134	template <typename DataType>
135	tuples::TaggedTuple<::Tags::FixedSource<Xcts::Tags::ConformalFactor<DataType>>>
136	ConstantDensityStar::variables(	1✔
137	const tnsr::I<DataType, 3, Frame::Inertial>& x,
138	tmpl::list<
139	::Tags::FixedSource<Xcts::Tags::ConformalFactor<DataType>>> /meta/)
140	const {
141	return {make_with_value<Scalar<DataType>>(x, 0.)};	2✔
142	}
143
144	template <typename DataType>
145	tuples::TaggedTuple<gr::Tags::EnergyDensity<DataType>>
146	ConstantDensityStar::variables(	×
147	const tnsr::I<DataType, 3, Frame::Inertial>& x,
148	tmpl::list<gr::Tags::EnergyDensity<DataType>> /meta/) const {
149	return {compute_piecewise(magnitude(x), radius_, density_, 0.)};	×
150	}
151
152	PUP::able::PUP_ID ConstantDensityStar::my_PUP_ID = 0; // NOLINT
153
154	bool operator==(const ConstantDensityStar& lhs,	2✔
155	const ConstantDensityStar& rhs) {
156	return lhs.density() == rhs.density() and lhs.radius() == rhs.radius();	2✔
157	}
158
159	bool operator!=(const ConstantDensityStar& lhs,	×
160	const ConstantDensityStar& rhs) {
161	return not(lhs == rhs);	×
162	}
163
164	#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
165
166	#define INSTANTIATE(_, data) \
167	template tuples::TaggedTuple<Xcts::Tags::ConformalFactor<DTYPE(data)>> \
168	ConstantDensityStar::variables( \
169	const tnsr::I<DTYPE(data), 3, Frame::Inertial>&, \
170	tmpl::list<Xcts::Tags::ConformalFactor<DTYPE(data)>>) const; \
171	template tuples::TaggedTuple< \
172	::Tags::Initial<Xcts::Tags::ConformalFactor<DTYPE(data)>>> \
173	ConstantDensityStar::variables( \
174	const tnsr::I<DTYPE(data), 3, Frame::Inertial>&, \
175	tmpl::list<::Tags::Initial<Xcts::Tags::ConformalFactor<DTYPE(data)>>>) \
176	const; \
177	template tuples::TaggedTuple< \
178	::Tags::Initial<::Tags::deriv<Xcts::Tags::ConformalFactor<DTYPE(data)>, \
179	tmpl::size_t<3>, Frame::Inertial>>> \
180	ConstantDensityStar::variables( \
181	const tnsr::I<DTYPE(data), 3, Frame::Inertial>&, \
182	tmpl::list<::Tags::Initial< \
183	::Tags::deriv<Xcts::Tags::ConformalFactor<DTYPE(data)>, \
184	tmpl::size_t<3>, Frame::Inertial>>>) const; \
185	template tuples::TaggedTuple< \
186	::Tags::FixedSource<Xcts::Tags::ConformalFactor<DTYPE(data)>>> \
187	ConstantDensityStar::variables( \
188	const tnsr::I<DTYPE(data), 3, Frame::Inertial>&, \
189	tmpl::list< \
190	::Tags::FixedSource<Xcts::Tags::ConformalFactor<DTYPE(data)>>>) \
191	const; \
192	template tuples::TaggedTuple<gr::Tags::EnergyDensity<DTYPE(data)>> \
193	ConstantDensityStar::variables( \
194	const tnsr::I<DTYPE(data), 3, Frame::Inertial>&, \
195	tmpl::list<gr::Tags::EnergyDensity<DTYPE(data)>>) const;
196
197	GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))
198
199	#undef DTYPE
200	#undef INSTANTIATE
201
202	} // namespace Xcts::Solutions

sxs-collaboration / spectre / 4175800929

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous