4138821718

Build Type

push

github

Committed by GitHub

Commit Message

Merge pull request #4714 from nilsdeppe/fix_jac_diag_missing_include

Run Details

63323 of 66009 relevant lines covered (95.93%)

403761.71 hits per line

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

67.12

/src/PointwiseFunctions/Hydro/EquationsOfState/Tabulated3d.cpp

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

#include "PointwiseFunctions/Hydro/EquationsOfState/Tabulated3d.hpp"

#include <limits>

#include "DataStructures/DataVector.hpp"  // IWYU pragma: keep
#include "DataStructures/Tensor/Tensor.hpp"
#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
#include "Utilities/ConstantExpressions.hpp"

// IWYU pragma: no_forward_declare Tensor

namespace EquationsOfState {

EQUATION_OF_STATE_MEMBER_DEFINITIONS(template <bool IsRelativistic>,
                                     Tabulated3D<IsRelativistic>, double, 3)
EQUATION_OF_STATE_MEMBER_DEFINITIONS(template <bool IsRelativistic>,
                                     Tabulated3D<IsRelativistic>, DataVector, 3)

template <bool IsRelativistic>
template <class DataType>
Scalar<DataType> Tabulated3D<IsRelativistic>::
    equilibrium_electron_fraction_from_density_temperature_impl(
        const Scalar<DataType>& rest_mass_density,
        const Scalar<DataType>& temperature) const {
  Scalar<DataType> electron_fraction = make_with_value<Scalar<DataType>>(
      rest_mass_density, electron_fraction_lower_bound());

  Scalar<DataType> converted_electron_fraction;
  Scalar<DataType> log_rest_mass_density;
  Scalar<DataType> log_temperature;

  convert_to_table_quantities(
      make_not_null(&converted_electron_fraction),
      make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
      electron_fraction, rest_mass_density, temperature);

  // Compute free-streaming beta-eq. electron fraction (from DeltaMu==0)

  if constexpr (std::is_same_v<DataType, double>) {
    const auto& log_rho = get(log_rest_mass_density);
    const auto& log_T = get(log_temperature);

    const auto f = [this, log_rho, log_T](const double ye) {

      const auto weights = interpolator_.get_weights(log_T, log_rho, ye);
      const auto interpolated_values =
          interpolator_.template interpolate<DeltaMu>(weights);

      return interpolated_values[0];
    };

    const auto root_from_lambda =
        RootFinder::toms748(f, electron_fraction_lower_bound(),
                            electron_fraction_upper_bound(), 1.0e-14, 1.0e-15);

    get(converted_electron_fraction) = root_from_lambda;

  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      const auto& log_rho = get(log_rest_mass_density)[s];
      const auto& log_T = get(log_temperature)[s];

      const auto f = [this, log_rho, log_T](const double ye) {

        const auto weights = interpolator_.get_weights(log_T, log_rho, ye);
        const auto interpolated_values =
            interpolator_.template interpolate<DeltaMu>(weights);

        return interpolated_values[0];
      };

      const auto root_from_lambda = RootFinder::toms748(
          f, electron_fraction_lower_bound(), electron_fraction_upper_bound(),
          1.0e-14, 1.0e-15);

      get(converted_electron_fraction)[s] = root_from_lambda;
    }
  }
  return converted_electron_fraction;
}

template <bool IsRelativistic>
std::unique_ptr<EquationOfState<IsRelativistic, 3>>
Tabulated3D<IsRelativistic>::get_clone() const {
  auto clone = std::make_unique<Tabulated3D<IsRelativistic>>(*this);
  return std::unique_ptr<EquationOfState<IsRelativistic, 3>>(std::move(clone));
}

template <bool IsRelativistic>
void Tabulated3D<IsRelativistic>::initialize(
    std::vector<double> electron_fraction, std::vector<double> log_density,
    std::vector<double> log_temperature, std::vector<double> table_data,
    double energy_shift, double enthalpy_minimum) {
  energy_shift_ = energy_shift;
  enthalpy_minimum_ = enthalpy_minimum;
  table_electron_fraction_ = std::move(electron_fraction);
  table_log_density_ = std::move(log_density);
  table_log_temperature_ = std::move(log_temperature);
  table_data_ = std::move(table_data);
  // Need to table

  Index<3> num_x_points;

  // The order is T, rho, Ye
  num_x_points[0] = table_log_temperature_.size();
  num_x_points[1] = table_log_density_.size();
  num_x_points[2] = table_electron_fraction_.size();

  std::array<gsl::span<double const>, 3> independent_data_view;

  independent_data_view[0] =
      gsl::span<double const>{table_log_temperature_.data(), num_x_points[0]};

  independent_data_view[1] =
      gsl::span<double const>{table_log_density_.data(), num_x_points[1]};

  independent_data_view[2] =
      gsl::span<double const>{table_electron_fraction_.data(), num_x_points[2]};

  interpolator_ = intrp::UniformMultiLinearSpanInterpolation<3, NumberOfVars>(
      independent_data_view, {table_data_.data(), table_data_.size()},
      num_x_points);
}

template <bool IsRelativistic>
bool Tabulated3D<IsRelativistic>::is_equal(
    const EquationOfState<IsRelativistic, 3>& rhs) const {
  const auto& derived_ptr =
      dynamic_cast<const Tabulated3D<IsRelativistic>* const>(&rhs);
  return derived_ptr != nullptr and *derived_ptr == *this;
}

template <bool IsRelativistic>
bool Tabulated3D<IsRelativistic>::operator==(
    const Tabulated3D<IsRelativistic>& rhs) const {
  bool result = true;
  result &= (rhs.enthalpy_minimum_ == this->enthalpy_minimum_);
  result &= (rhs.energy_shift_ == this->energy_shift_);
  result &= (rhs.table_electron_fraction_ == this->table_electron_fraction_);
  result &= (rhs.table_log_density_ == this->table_log_density_);
  result &= (rhs.table_log_temperature_ == this->table_log_temperature_);
  result &= (rhs.table_data_ == this->table_data_);

  return result;
}

template <bool IsRelativistic>
bool Tabulated3D<IsRelativistic>::operator!=(
    const Tabulated3D<IsRelativistic>& rhs) const {
  return not(*this == rhs);
}

template <bool IsRelativistic>
Tabulated3D<IsRelativistic>::Tabulated3D(CkMigrateMessage* msg)
    : EquationOfState<IsRelativistic, 3>(msg) {}
template <bool IsRelativistic>
template <class DataType>
void Tabulated3D<IsRelativistic>::enforce_physicality(
    Scalar<DataType>& electron_fraction, Scalar<DataType>& rest_mass_density,
    Scalar<DataType>& temperature) const {
  if constexpr (std::is_same_v<DataType, double>) {
    get(rest_mass_density) = std::max(
        std::min(get(rest_mass_density), rest_mass_density_upper_bound()),
        rest_mass_density_lower_bound());

    get(electron_fraction) = std::max(
        std::min(get(electron_fraction), electron_fraction_upper_bound()),
        electron_fraction_lower_bound());

    get(temperature) =
        std::max(std::min(get(temperature), temperature_upper_bound()),
                 temperature_lower_bound());

  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      get(rest_mass_density)[s] = std::max(
          std::min(get(rest_mass_density)[s], rest_mass_density_upper_bound()),
          rest_mass_density_lower_bound());

      get(electron_fraction)[s] = std::max(
          std::min(get(electron_fraction)[s], electron_fraction_upper_bound()),
          electron_fraction_lower_bound());

      get(temperature)[s] =
          std::max(std::min(get(temperature)[s], temperature_upper_bound()),
                   temperature_lower_bound());
    }
  }
}


template <bool IsRelativistic>
template <class DataType>
Scalar<DataType>
Tabulated3D<IsRelativistic>::pressure_from_density_and_energy_impl(
    const Scalar<DataType>& rest_mass_density,
    const Scalar<DataType>& specific_internal_energy,
    const Scalar<DataType>& electron_fraction) const {
  auto temperature = temperature_from_density_and_energy_impl(
      rest_mass_density, specific_internal_energy, electron_fraction);

  return pressure_from_density_and_temperature_impl(
      rest_mass_density, temperature, electron_fraction);
}

template <bool IsRelativistic>
template <class DataType>
Scalar<DataType>
Tabulated3D<IsRelativistic>::pressure_from_density_and_temperature_impl(
    const Scalar<DataType>& rest_mass_density,
    const Scalar<DataType>& temperature,
    const Scalar<DataType>& electron_fraction) const {
  Scalar<DataType> converted_electron_fraction;
  Scalar<DataType> log_rest_mass_density;
  Scalar<DataType> log_temperature;

  convert_to_table_quantities(
      make_not_null(&converted_electron_fraction),
      make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
      electron_fraction, rest_mass_density, temperature);

  Scalar<DataType> pressure =
      make_with_value<Scalar<DataType>>(get(rest_mass_density), 0.0);

  if constexpr (std::is_same_v<DataType, double>) {
    auto weights = interpolator_.get_weights(get(log_temperature),
                                             get(log_rest_mass_density),
                                             get(converted_electron_fraction));
    auto interpolated_state =
        interpolator_.template interpolate<Pressure>(weights);
    get(pressure) = std::exp(interpolated_state[0]);

  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      auto weights = interpolator_.get_weights(
          get(log_temperature)[s], get(log_rest_mass_density)[s],
          get(converted_electron_fraction)[s]);
      auto interpolated_state =
          interpolator_.template interpolate<Pressure>(weights);
      get(pressure)[s] = std::exp(interpolated_state[0]);
    }
  }

  return pressure;
}


template <bool IsRelativistic>
void Tabulated3D<IsRelativistic>::pup(PUP::er& p) {
  EquationOfState<IsRelativistic, 3>::pup(p);
  p | energy_shift_;
  p | enthalpy_minimum_;
  p | table_electron_fraction_;
  p | table_log_density_;
  p | table_log_temperature_;
  p | table_data_;
}

template <bool IsRelativistic>
template <class DataType>
Scalar<DataType>
Tabulated3D<IsRelativistic>::temperature_from_density_and_energy_impl(
    const Scalar<DataType>& rest_mass_density,
    const Scalar<DataType>& specific_internal_energy,
    const Scalar<DataType>& electron_fraction) const {
  Scalar<DataType> converted_electron_fraction;
  Scalar<DataType> log_rest_mass_density;

  Scalar<DataType> log_temperature;
  Scalar<DataType> temperature;

  temperature = make_with_value<Scalar<DataType>>(rest_mass_density,
                                                  temperature_lower_bound());

  convert_to_table_quantities(
      make_not_null(&converted_electron_fraction),
      make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
      electron_fraction, rest_mass_density, temperature);

  // Check bounds on eps, note that eps may be negative
  Scalar<DataType> log_specific_internal_energy = specific_internal_energy;

  if constexpr (std::is_same_v<DataType, double>) {
    get(log_specific_internal_energy) =
        std::max(std::min(get(specific_internal_energy),
                          specific_internal_energy_upper_bound(
                              get(rest_mass_density), get(electron_fraction))),
                 specific_internal_energy_lower_bound(get(rest_mass_density),
                                                      get(electron_fraction)));
  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      get(log_specific_internal_energy)[s] = std::max(
          std::min(get(specific_internal_energy)[s],
                   specific_internal_energy_upper_bound(
                       get(rest_mass_density)[s], get(electron_fraction)[s])),
          specific_internal_energy_lower_bound(get(rest_mass_density)[s],
                                               get(electron_fraction)[s]));
    }
  }

  // Correct for negative eps
  get(log_specific_internal_energy) -= energy_shift_;
  get(log_specific_internal_energy) = log(get(log_specific_internal_energy));

  if constexpr (std::is_same_v<DataType, double>) {
    const auto& log_eps = get(log_specific_internal_energy);
    const auto& log_rho = get(log_rest_mass_density);
    const auto& ye = get(converted_electron_fraction);

    // Root-finding appropriate between reference density and maximum density
    // We can use x=0 and x=x_max as bounds
    const auto f = [this, log_eps, log_rho, ye](const double log_T) {

      const auto weights = interpolator_.get_weights(log_T, log_rho, ye);
      const auto interpolated_values =
          interpolator_.template interpolate<Epsilon>(weights);

      return log_eps - interpolated_values[0];
    };

    const auto root_from_lambda = RootFinder::toms748(
        f, table_log_temperature_.front(),
        upper_bound_tolerance_ * table_log_temperature_.back(), 1.0e-14,
        1.0e-15);

    get(temperature) = exp(root_from_lambda);

  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      const auto& log_eps = get(log_specific_internal_energy)[s];
      const auto& log_rho = get(log_rest_mass_density)[s];
      const auto& ye = get(converted_electron_fraction)[s];

      // Root-finding appropriate between reference density and maximum density
      // We can use x=0 and x=x_max as bounds
      const auto f = [this, log_eps, log_rho, ye](const double log_T) {

        const auto weights = interpolator_.get_weights(log_T, log_rho, ye);
        const auto interpolated_values =
            interpolator_.template interpolate<Epsilon>(weights);

        return log_eps - interpolated_values[0];
      };
      const auto root_from_lambda = RootFinder::toms748(
          f, table_log_temperature_.front(),
          upper_bound_tolerance_ * table_log_temperature_.back(), 1.0e-14,
          1.0e-15);

      get(temperature)[s] = exp(root_from_lambda);
    }
  }
  return temperature;
}

template <bool IsRelativistic>
template <class DataType>
Scalar<DataType> Tabulated3D<IsRelativistic>::
    specific_internal_energy_from_density_and_temperature_impl(
        const Scalar<DataType>& rest_mass_density,
        const Scalar<DataType>& temperature,
        const Scalar<DataType>& electron_fraction) const {
  Scalar<DataType> converted_electron_fraction;
  Scalar<DataType> log_rest_mass_density;
  Scalar<DataType> log_temperature;

  convert_to_table_quantities(
      make_not_null(&converted_electron_fraction),
      make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
      electron_fraction, rest_mass_density, temperature);

  Scalar<DataType> specific_internal_energy =
      make_with_value<Scalar<DataType>>(get(rest_mass_density), 0.0);

  if constexpr (std::is_same_v<DataType, double>) {
    auto weights = interpolator_.get_weights(get(log_temperature),
                                             get(log_rest_mass_density),
                                             get(converted_electron_fraction));
    auto interpolated_state =
        interpolator_.template interpolate<Epsilon>(weights);
    get(specific_internal_energy) =
        std::exp(interpolated_state[0]) + energy_shift_;
  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      auto weights = interpolator_.get_weights(
          get(log_temperature)[s], get(log_rest_mass_density)[s],
          get(converted_electron_fraction)[s]);
      auto interpolated_state =
          interpolator_.template interpolate<Epsilon>(weights);
      get(specific_internal_energy)[s] =
          std::exp(interpolated_state[0]) + energy_shift_;
    }
  }

  return specific_internal_energy;
}

template <bool IsRelativistic>
template <class DataType>
Scalar<DataType> Tabulated3D<IsRelativistic>::
    sound_speed_squared_from_density_and_temperature_impl(
        const Scalar<DataType>& rest_mass_density,
        const Scalar<DataType>& temperature,
        const Scalar<DataType>& electron_fraction) const {
  Scalar<DataType> converted_electron_fraction;
  Scalar<DataType> log_rest_mass_density;
  Scalar<DataType> log_temperature;

  convert_to_table_quantities(
      make_not_null(&converted_electron_fraction),
      make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
      electron_fraction, rest_mass_density, temperature);

  Scalar<DataType> cs2 =
      make_with_value<Scalar<DataType>>(get(rest_mass_density), 0.0);

  if constexpr (std::is_same_v<DataType, double>) {
    auto weights = interpolator_.get_weights(get(log_temperature),
                                             get(log_rest_mass_density),
                                             get(converted_electron_fraction));
    auto interpolated_state =
        interpolator_.template interpolate<CsSquared>(weights);
    get(cs2) = interpolated_state[0];

  } else if constexpr (std::is_same_v<DataType, DataVector>) {
    for (size_t s = 0; s < electron_fraction.size(); ++s) {
      auto weights = interpolator_.get_weights(
          get(log_temperature)[s], get(log_rest_mass_density)[s],
          get(converted_electron_fraction)[s]);
      auto interpolated_state =
          interpolator_.template interpolate<CsSquared>(weights);
      get(cs2)[s] = interpolated_state[0];
    }
  }

  return cs2;
}

template <bool IsRelativistic>
double Tabulated3D<IsRelativistic>::specific_internal_energy_lower_bound(
    const double rest_mass_density, const double electron_fraction) const {
  double converted_electron_fraction =
      std::min(std::max(electron_fraction_lower_bound(), electron_fraction),
               electron_fraction_upper_bound());

  double log_rest_mass_density =
      std::min(std::max(rest_mass_density_lower_bound(), rest_mass_density),
               rest_mass_density_upper_bound());

  log_rest_mass_density = log(log_rest_mass_density);

  auto weights = interpolator_.get_weights(log(temperature_lower_bound()),
                                           log_rest_mass_density,
                                           converted_electron_fraction);
  auto interpolated_state =
      interpolator_.template interpolate<Epsilon>(weights);

  return exp(interpolated_state[0]) + energy_shift_;
}

template <bool IsRelativistic>
double Tabulated3D<IsRelativistic>::specific_internal_energy_upper_bound(
    const double rest_mass_density, const double electron_fraction) const {
  double converted_electron_fraction =
      std::min(std::max(electron_fraction_lower_bound(), electron_fraction),
               electron_fraction_upper_bound());

  double log_rest_mass_density =
      std::min(std::max(rest_mass_density_lower_bound(), rest_mass_density),
               rest_mass_density_upper_bound());

  log_rest_mass_density = log(log_rest_mass_density);

  auto weights = interpolator_.get_weights(
      log(upper_bound_tolerance_ * temperature_upper_bound()),
      log_rest_mass_density, converted_electron_fraction);
  auto interpolated_state =
      interpolator_.template interpolate<Epsilon>(weights);

  return exp(interpolated_state[0]) + energy_shift_;
}

template <bool IsRelativistic>
Tabulated3D<IsRelativistic>::
    Tabulated3D(  // NOLINTNEXTLINE(performance-unnecessary-value-param)
        std::vector<double> electron_fraction,
        // NOLINTNEXTLINE(performance-unnecessary-value-param)
        std::vector<double> log_density,
        // NOLINTNEXTLINE(performance-unnecessary-value-param)
        std::vector<double> log_temperature,
        // NOLINTNEXTLINE(performance-unnecessary-value-param)
        std::vector<double> table_data, double energy_shift,
        double enthalpy_minimum) {
  initialize(std::move(electron_fraction), std::move(log_density),
             std::move(log_temperature), std::move(table_data), energy_shift,
             enthalpy_minimum);
}
}  // namespace EquationsOfState

template class EquationsOfState::Tabulated3D<true>;
template class EquationsOfState::Tabulated3D<false>;


1	// Distributed under the MIT License.
2	// See LICENSE.txt for details.
3
4	#include "PointwiseFunctions/Hydro/EquationsOfState/Tabulated3d.hpp"
5
6	#include <limits>
7
8	#include "DataStructures/DataVector.hpp" // IWYU pragma: keep
9	#include "DataStructures/Tensor/Tensor.hpp"
10	#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
11	#include "Utilities/ConstantExpressions.hpp"
12
13	// IWYU pragma: no_forward_declare Tensor
14
15	namespace EquationsOfState {
16
17	EQUATION_OF_STATE_MEMBER_DEFINITIONS(template <bool IsRelativistic>,	3✔
18	Tabulated3D<IsRelativistic>, double, 3)
19	EQUATION_OF_STATE_MEMBER_DEFINITIONS(template <bool IsRelativistic>,	7✔
20	Tabulated3D<IsRelativistic>, DataVector, 3)
21
22	template <bool IsRelativistic>
23	template <class DataType>
24	Scalar<DataType> Tabulated3D<IsRelativistic>::	×
25	equilibrium_electron_fraction_from_density_temperature_impl(
26	const Scalar<DataType>& rest_mass_density,
27	const Scalar<DataType>& temperature) const {
28	Scalar<DataType> electron_fraction = make_with_value<Scalar<DataType>>(	×
29	rest_mass_density, electron_fraction_lower_bound());
30
31	Scalar<DataType> converted_electron_fraction;	×
32	Scalar<DataType> log_rest_mass_density;	×
33	Scalar<DataType> log_temperature;	×
34
35	convert_to_table_quantities(	×
36	make_not_null(&converted_electron_fraction),
37	make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
38	electron_fraction, rest_mass_density, temperature);
39
40	// Compute free-streaming beta-eq. electron fraction (from DeltaMu==0)
41
42	if constexpr (std::is_same_v<DataType, double>) {
43	const auto& log_rho = get(log_rest_mass_density);	×
44	const auto& log_T = get(log_temperature);	×
45
46	const auto f = [this, log_rho, log_T](const double ye) {	×
47
48	const auto weights = interpolator_.get_weights(log_T, log_rho, ye);	×
49	const auto interpolated_values =	×
50	interpolator_.template interpolate<DeltaMu>(weights);
51
52	return interpolated_values[0];	×
53	};
54
55	const auto root_from_lambda =
56	RootFinder::toms748(f, electron_fraction_lower_bound(),	×
57	electron_fraction_upper_bound(), 1.0e-14, 1.0e-15);
58
59	get(converted_electron_fraction) = root_from_lambda;	×
60
61	} else if constexpr (std::is_same_v<DataType, DataVector>) {
62	for (size_t s = 0; s < electron_fraction.size(); ++s) {	×
63	const auto& log_rho = get(log_rest_mass_density)[s];	×
64	const auto& log_T = get(log_temperature)[s];	×
65
66	const auto f = [this, log_rho, log_T](const double ye) {	×
67
68	const auto weights = interpolator_.get_weights(log_T, log_rho, ye);	×
69	const auto interpolated_values =	×
70	interpolator_.template interpolate<DeltaMu>(weights);
71
72	return interpolated_values[0];	×
73	};
74
75	const auto root_from_lambda = RootFinder::toms748(	×
76	f, electron_fraction_lower_bound(), electron_fraction_upper_bound(),
77	1.0e-14, 1.0e-15);
78
79	get(converted_electron_fraction)[s] = root_from_lambda;	×
80	}
81	}
82	return converted_electron_fraction;	×
83	}
84
85	template <bool IsRelativistic>
86	std::unique_ptr<EquationOfState<IsRelativistic, 3>>
87	Tabulated3D<IsRelativistic>::get_clone() const {	×
88	auto clone = std::make_unique<Tabulated3D<IsRelativistic>>(*this);	×
89	return std::unique_ptr<EquationOfState<IsRelativistic, 3>>(std::move(clone));	×
90	}
91
92	template <bool IsRelativistic>
93	void Tabulated3D<IsRelativistic>::initialize(	1✔
94	std::vector<double> electron_fraction, std::vector<double> log_density,
95	std::vector<double> log_temperature, std::vector<double> table_data,
96	double energy_shift, double enthalpy_minimum) {
97	energy_shift_ = energy_shift;	1✔
98	enthalpy_minimum_ = enthalpy_minimum;	1✔
99	table_electron_fraction_ = std::move(electron_fraction);	1✔
100	table_log_density_ = std::move(log_density);	1✔
101	table_log_temperature_ = std::move(log_temperature);	1✔
102	table_data_ = std::move(table_data);	1✔
103	// Need to table
104
105	Index<3> num_x_points;	1✔
106
107	// The order is T, rho, Ye
108	num_x_points[0] = table_log_temperature_.size();	1✔
109	num_x_points[1] = table_log_density_.size();	1✔
110	num_x_points[2] = table_electron_fraction_.size();	1✔
111
112	std::array<gsl::span<double const>, 3> independent_data_view;	2✔
113
114	independent_data_view[0] =	1✔
115	gsl::span<double const>{table_log_temperature_.data(), num_x_points[0]};	1✔
116
117	independent_data_view[1] =	1✔
118	gsl::span<double const>{table_log_density_.data(), num_x_points[1]};	1✔
119
120	independent_data_view[2] =	1✔
121	gsl::span<double const>{table_electron_fraction_.data(), num_x_points[2]};	1✔
122
123	interpolator_ = intrp::UniformMultiLinearSpanInterpolation<3, NumberOfVars>(	1✔
124	independent_data_view, {table_data_.data(), table_data_.size()},	1✔
125	num_x_points);
126	}	1✔
127
128	template <bool IsRelativistic>
129	bool Tabulated3D<IsRelativistic>::is_equal(	×
130	const EquationOfState<IsRelativistic, 3>& rhs) const {
131	const auto& derived_ptr =	×
132	dynamic_cast<const Tabulated3D<IsRelativistic>* const>(&rhs);	×
133	return derived_ptr != nullptr and derived_ptr == this;	×
134	}
135
136	template <bool IsRelativistic>
137	bool Tabulated3D<IsRelativistic>::operator==(	×
138	const Tabulated3D<IsRelativistic>& rhs) const {
139	bool result = true;	×
140	result &= (rhs.enthalpy_minimum_ == this->enthalpy_minimum_);	×
141	result &= (rhs.energy_shift_ == this->energy_shift_);	×
142	result &= (rhs.table_electron_fraction_ == this->table_electron_fraction_);	×
143	result &= (rhs.table_log_density_ == this->table_log_density_);	×
144	result &= (rhs.table_log_temperature_ == this->table_log_temperature_);	×
145	result &= (rhs.table_data_ == this->table_data_);	×
146
147	return result;	×
148	}
149
150	template <bool IsRelativistic>
151	bool Tabulated3D<IsRelativistic>::operator!=(	×
152	const Tabulated3D<IsRelativistic>& rhs) const {
153	return not(*this == rhs);	×
154	}
155
156	template <bool IsRelativistic>
157	Tabulated3D<IsRelativistic>::Tabulated3D(CkMigrateMessage* msg)	×
158	: EquationOfState<IsRelativistic, 3>(msg) {}	×
159	template <bool IsRelativistic>
160	template <class DataType>
161	void Tabulated3D<IsRelativistic>::enforce_physicality(	10✔
162	Scalar<DataType>& electron_fraction, Scalar<DataType>& rest_mass_density,
163	Scalar<DataType>& temperature) const {
164	if constexpr (std::is_same_v<DataType, double>) {
165	get(rest_mass_density) = std::max(	6✔
166	std::min(get(rest_mass_density), rest_mass_density_upper_bound()),	6✔
167	rest_mass_density_lower_bound());
168
169	get(electron_fraction) = std::max(	6✔
170	std::min(get(electron_fraction), electron_fraction_upper_bound()),	6✔
171	electron_fraction_lower_bound());
172
173	get(temperature) =	3✔
174	std::max(std::min(get(temperature), temperature_upper_bound()),	6✔
175	temperature_lower_bound());
176
177	} else if constexpr (std::is_same_v<DataType, DataVector>) {
178	for (size_t s = 0; s < electron_fraction.size(); ++s) {	21✔
179	get(rest_mass_density)[s] = std::max(	21✔
180	std::min(get(rest_mass_density)[s], rest_mass_density_upper_bound()),	21✔
181	rest_mass_density_lower_bound());
182
183	get(electron_fraction)[s] = std::max(	21✔
184	std::min(get(electron_fraction)[s], electron_fraction_upper_bound()),	21✔
185	electron_fraction_lower_bound());
186
187	get(temperature)[s] =	14✔
188	std::max(std::min(get(temperature)[s], temperature_upper_bound()),	28✔
189	temperature_lower_bound());
190	}
191	}
192	}	10✔
193
194
195	template <bool IsRelativistic>
196	template <class DataType>
197	Scalar<DataType>
198	Tabulated3D<IsRelativistic>::pressure_from_density_and_energy_impl(	×
199	const Scalar<DataType>& rest_mass_density,
200	const Scalar<DataType>& specific_internal_energy,
201	const Scalar<DataType>& electron_fraction) const {
202	auto temperature = temperature_from_density_and_energy_impl(	×
203	rest_mass_density, specific_internal_energy, electron_fraction);
204
205	return pressure_from_density_and_temperature_impl(
206	rest_mass_density, temperature, electron_fraction);	×
207	}
208
209	template <bool IsRelativistic>
210	template <class DataType>
211	Scalar<DataType>
212	Tabulated3D<IsRelativistic>::pressure_from_density_and_temperature_impl(	2✔
213	const Scalar<DataType>& rest_mass_density,
214	const Scalar<DataType>& temperature,
215	const Scalar<DataType>& electron_fraction) const {
216	Scalar<DataType> converted_electron_fraction;	2✔
217	Scalar<DataType> log_rest_mass_density;	2✔
218	Scalar<DataType> log_temperature;	2✔
219
220	convert_to_table_quantities(	2✔
221	make_not_null(&converted_electron_fraction),
222	make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
223	electron_fraction, rest_mass_density, temperature);
224
225	Scalar<DataType> pressure =
226	make_with_value<Scalar<DataType>>(get(rest_mass_density), 0.0);	4✔
227
228	if constexpr (std::is_same_v<DataType, double>) {
229	auto weights = interpolator_.get_weights(get(log_temperature),	2✔
230	get(log_rest_mass_density),	1✔
231	get(converted_electron_fraction));	1✔
232	auto interpolated_state =	1✔
233	interpolator_.template interpolate<Pressure>(weights);
234	get(pressure) = std::exp(interpolated_state[0]);	1✔
235
236	} else if constexpr (std::is_same_v<DataType, DataVector>) {
237	for (size_t s = 0; s < electron_fraction.size(); ++s) {	3✔
238	auto weights = interpolator_.get_weights(	1✔
239	get(log_temperature)[s], get(log_rest_mass_density)[s],	6✔
240	get(converted_electron_fraction)[s]);	2✔
241	auto interpolated_state =	1✔
242	interpolator_.template interpolate<Pressure>(weights);
243	get(pressure)[s] = std::exp(interpolated_state[0]);	3✔
244	}
245	}
246
247	return pressure;	3✔
248	}
249
250
251	template <bool IsRelativistic>
252	void Tabulated3D<IsRelativistic>::pup(PUP::er& p) {	×
253	EquationOfState<IsRelativistic, 3>::pup(p);	×
254	p \| energy_shift_;	×
255	p \| enthalpy_minimum_;	×
256	p \| table_electron_fraction_;	×
257	p \| table_log_density_;	×
258	p \| table_log_temperature_;	×
259	p \| table_data_;	×
260	}	×
261
262	template <bool IsRelativistic>
263	template <class DataType>
264	Scalar<DataType>
265	Tabulated3D<IsRelativistic>::temperature_from_density_and_energy_impl(	2✔
266	const Scalar<DataType>& rest_mass_density,
267	const Scalar<DataType>& specific_internal_energy,
268	const Scalar<DataType>& electron_fraction) const {
269	Scalar<DataType> converted_electron_fraction;	4✔
270	Scalar<DataType> log_rest_mass_density;	4✔
271
272	Scalar<DataType> log_temperature;	4✔
273	Scalar<DataType> temperature;	2✔
274
275	temperature = make_with_value<Scalar<DataType>>(rest_mass_density,	4✔
276	temperature_lower_bound());
277
278	convert_to_table_quantities(	2✔
279	make_not_null(&converted_electron_fraction),
280	make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
281	electron_fraction, rest_mass_density, temperature);
282
283	// Check bounds on eps, note that eps may be negative
284	Scalar<DataType> log_specific_internal_energy = specific_internal_energy;	4✔
285
286	if constexpr (std::is_same_v<DataType, double>) {
287	get(log_specific_internal_energy) =	×
288	std::max(std::min(get(specific_internal_energy),	×
289	specific_internal_energy_upper_bound(
290	get(rest_mass_density), get(electron_fraction))),	×
291	specific_internal_energy_lower_bound(get(rest_mass_density),	×
292	get(electron_fraction)));	×
293	} else if constexpr (std::is_same_v<DataType, DataVector>) {
294	for (size_t s = 0; s < electron_fraction.size(); ++s) {	6✔
295	get(log_specific_internal_energy)[s] = std::max(	6✔
296	std::min(get(specific_internal_energy)[s],	8✔
297	specific_internal_energy_upper_bound(
298	get(rest_mass_density)[s], get(electron_fraction)[s])),	8✔
299	specific_internal_energy_lower_bound(get(rest_mass_density)[s],	4✔
300	get(electron_fraction)[s]));	4✔
301	}
302	}
303
304	// Correct for negative eps
305	get(log_specific_internal_energy) -= energy_shift_;	4✔
306	get(log_specific_internal_energy) = log(get(log_specific_internal_energy));	4✔
307
308	if constexpr (std::is_same_v<DataType, double>) {
309	const auto& log_eps = get(log_specific_internal_energy);	×
310	const auto& log_rho = get(log_rest_mass_density);	×
311	const auto& ye = get(converted_electron_fraction);	×
312
313	// Root-finding appropriate between reference density and maximum density
314	// We can use x=0 and x=x_max as bounds
315	const auto f = [this, log_eps, log_rho, ye](const double log_T) {	×
316
317	const auto weights = interpolator_.get_weights(log_T, log_rho, ye);	×
318	const auto interpolated_values =	×
319	interpolator_.template interpolate<Epsilon>(weights);
320
321	return log_eps - interpolated_values[0];	×
322	};
323
324	const auto root_from_lambda = RootFinder::toms748(	×
325	f, table_log_temperature_.front(),	×
326	upper_bound_tolerance_ * table_log_temperature_.back(), 1.0e-14,	×
327	1.0e-15);
328
329	get(temperature) = exp(root_from_lambda);	×
330
331	} else if constexpr (std::is_same_v<DataType, DataVector>) {
332	for (size_t s = 0; s < electron_fraction.size(); ++s) {	6✔
333	const auto& log_eps = get(log_specific_internal_energy)[s];	4✔
334	const auto& log_rho = get(log_rest_mass_density)[s];	4✔
335	const auto& ye = get(converted_electron_fraction)[s];	2✔
336
337	// Root-finding appropriate between reference density and maximum density
338	// We can use x=0 and x=x_max as bounds
339	const auto f = [this, log_eps, log_rho, ye](const double log_T) {	20✔
340
341	const auto weights = interpolator_.get_weights(log_T, log_rho, ye);	6✔
342	const auto interpolated_values =	6✔
343	interpolator_.template interpolate<Epsilon>(weights);
344
345	return log_eps - interpolated_values[0];	6✔
346	};
347	const auto root_from_lambda = RootFinder::toms748(	2✔
348	f, table_log_temperature_.front(),	2✔
349	upper_bound_tolerance_ * table_log_temperature_.back(), 1.0e-14,	2✔
350	1.0e-15);
351
352	get(temperature)[s] = exp(root_from_lambda);	6✔
353	}
354	}
355	return temperature;	4✔
356	}
357
358	template <bool IsRelativistic>
359	template <class DataType>
360	Scalar<DataType> Tabulated3D<IsRelativistic>::	4✔
361	specific_internal_energy_from_density_and_temperature_impl(
362	const Scalar<DataType>& rest_mass_density,
363	const Scalar<DataType>& temperature,
364	const Scalar<DataType>& electron_fraction) const {
365	Scalar<DataType> converted_electron_fraction;	6✔
366	Scalar<DataType> log_rest_mass_density;	6✔
367	Scalar<DataType> log_temperature;	6✔
368
369	convert_to_table_quantities(	4✔
370	make_not_null(&converted_electron_fraction),
371	make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
372	electron_fraction, rest_mass_density, temperature);
373
374	Scalar<DataType> specific_internal_energy =
375	make_with_value<Scalar<DataType>>(get(rest_mass_density), 0.0);	8✔
376
377	if constexpr (std::is_same_v<DataType, double>) {
378	auto weights = interpolator_.get_weights(get(log_temperature),	2✔
379	get(log_rest_mass_density),	1✔
380	get(converted_electron_fraction));	1✔
381	auto interpolated_state =	1✔
382	interpolator_.template interpolate<Epsilon>(weights);
383	get(specific_internal_energy) =	1✔
384	std::exp(interpolated_state[0]) + energy_shift_;	1✔
385	} else if constexpr (std::is_same_v<DataType, DataVector>) {
386	for (size_t s = 0; s < electron_fraction.size(); ++s) {	9✔
387	auto weights = interpolator_.get_weights(	3✔
388	get(log_temperature)[s], get(log_rest_mass_density)[s],	12✔
389	get(converted_electron_fraction)[s]);	6✔
390	auto interpolated_state =	3✔
391	interpolator_.template interpolate<Epsilon>(weights);
392	get(specific_internal_energy)[s] =	6✔
393	std::exp(interpolated_state[0]) + energy_shift_;	3✔
394	}
395	}
396
397	return specific_internal_energy;	7✔
398	}
399
400	template <bool IsRelativistic>
401	template <class DataType>
402	Scalar<DataType> Tabulated3D<IsRelativistic>::	2✔
403	sound_speed_squared_from_density_and_temperature_impl(
404	const Scalar<DataType>& rest_mass_density,
405	const Scalar<DataType>& temperature,
406	const Scalar<DataType>& electron_fraction) const {
407	Scalar<DataType> converted_electron_fraction;	2✔
408	Scalar<DataType> log_rest_mass_density;	2✔
409	Scalar<DataType> log_temperature;	2✔
410
411	convert_to_table_quantities(	2✔
412	make_not_null(&converted_electron_fraction),
413	make_not_null(&log_rest_mass_density), make_not_null(&log_temperature),
414	electron_fraction, rest_mass_density, temperature);
415
416	Scalar<DataType> cs2 =
417	make_with_value<Scalar<DataType>>(get(rest_mass_density), 0.0);	4✔
418
419	if constexpr (std::is_same_v<DataType, double>) {
420	auto weights = interpolator_.get_weights(get(log_temperature),	2✔
421	get(log_rest_mass_density),	1✔
422	get(converted_electron_fraction));	1✔
423	auto interpolated_state =	1✔
424	interpolator_.template interpolate<CsSquared>(weights);
425	get(cs2) = interpolated_state[0];	1✔
426
427	} else if constexpr (std::is_same_v<DataType, DataVector>) {
428	for (size_t s = 0; s < electron_fraction.size(); ++s) {	3✔
429	auto weights = interpolator_.get_weights(	1✔
430	get(log_temperature)[s], get(log_rest_mass_density)[s],	4✔
431	get(converted_electron_fraction)[s]);	2✔
432	auto interpolated_state =	1✔
433	interpolator_.template interpolate<CsSquared>(weights);
434	get(cs2)[s] = interpolated_state[0];	3✔
435	}
436	}
437
438	return cs2;	3✔
439	}
440
441	template <bool IsRelativistic>
442	double Tabulated3D<IsRelativistic>::specific_internal_energy_lower_bound(	2✔
443	const double rest_mass_density, const double electron_fraction) const {
444	double converted_electron_fraction =	2✔
445	std::min(std::max(electron_fraction_lower_bound(), electron_fraction),	2✔
446	electron_fraction_upper_bound());
447
448	double log_rest_mass_density =	2✔
449	std::min(std::max(rest_mass_density_lower_bound(), rest_mass_density),	2✔
450	rest_mass_density_upper_bound());
451
452	log_rest_mass_density = log(log_rest_mass_density);	2✔
453
454	auto weights = interpolator_.get_weights(log(temperature_lower_bound()),	2✔
455	log_rest_mass_density,
456	converted_electron_fraction);
457	auto interpolated_state =	2✔
458	interpolator_.template interpolate<Epsilon>(weights);
459
460	return exp(interpolated_state[0]) + energy_shift_;	2✔
461	}
462
463	template <bool IsRelativistic>
464	double Tabulated3D<IsRelativistic>::specific_internal_energy_upper_bound(	2✔
465	const double rest_mass_density, const double electron_fraction) const {
466	double converted_electron_fraction =	2✔
467	std::min(std::max(electron_fraction_lower_bound(), electron_fraction),	2✔
468	electron_fraction_upper_bound());
469
470	double log_rest_mass_density =	2✔
471	std::min(std::max(rest_mass_density_lower_bound(), rest_mass_density),	2✔
472	rest_mass_density_upper_bound());
473
474	log_rest_mass_density = log(log_rest_mass_density);	2✔
475
476	auto weights = interpolator_.get_weights(	4✔
477	log(upper_bound_tolerance_ * temperature_upper_bound()),	2✔
478	log_rest_mass_density, converted_electron_fraction);
479	auto interpolated_state =	2✔
480	interpolator_.template interpolate<Epsilon>(weights);
481
482	return exp(interpolated_state[0]) + energy_shift_;	2✔
483	}
484
485	template <bool IsRelativistic>
486	Tabulated3D<IsRelativistic>::	1✔
487	Tabulated3D( // NOLINTNEXTLINE(performance-unnecessary-value-param)
488	std::vector<double> electron_fraction,
489	// NOLINTNEXTLINE(performance-unnecessary-value-param)
490	std::vector<double> log_density,
491	// NOLINTNEXTLINE(performance-unnecessary-value-param)
492	std::vector<double> log_temperature,
493	// NOLINTNEXTLINE(performance-unnecessary-value-param)
494	std::vector<double> table_data, double energy_shift,
495	double enthalpy_minimum) {	1✔
496	initialize(std::move(electron_fraction), std::move(log_density),	2✔
497	std::move(log_temperature), std::move(table_data), energy_shift,	2✔
498	enthalpy_minimum);
499	}	1✔
500	} // namespace EquationsOfState
501
502	template class EquationsOfState::Tabulated3D<true>;
503	template class EquationsOfState::Tabulated3D<false>;
504

sxs-collaboration / spectre / 4138821718

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