Skip to content

Commit

Permalink
Add tabulated equation of state class.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ermost
ermost committed Jan 30, 2023
1 parent 77a7a66 commit 4dd2ab1
Show file tree
Hide file tree
Showing 7 changed files with 1,032 additions and 18 deletions.
4 changes: 3 additions & 1 deletion src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,8 @@ class MultiLinearSpanInterpolation {
std::make_index_sequence<Dimension>{});
}

MultiLinearSpanInterpolation() = default;

MultiLinearSpanInterpolation(
std::array<gsl::span<const double>, Dimension> x_,
gsl::span<const double> y_, Index<Dimension> number_of_points__);
Expand Down Expand Up @@ -221,7 +223,7 @@ size_t MultiLinearSpanInterpolation<
// Use linear extrapolation based of the lowest
// two points in the table
ASSERT(allow_extrapolation_below_data_[which_dimension] or
UNLIKELY(relative_coordinate > 0.),
UNLIKELY(relative_coordinate >= 0.),
"Interpolation exceeds lower table bounds.");

// We are exceeding the table bounds:
Expand Down
2 changes: 2 additions & 0 deletions src/PointwiseFunctions/Hydro/EquationsOfState/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ spectre_target_sources(
PolytropicFluid.cpp
RegisterDerivedWithCharm.cpp
Spectral.cpp
Tabulated3d.cpp
)

spectre_target_headers(
Expand All @@ -28,6 +29,7 @@ spectre_target_headers(
PolytropicFluid.hpp
RegisterDerivedWithCharm.hpp
Spectral.hpp
Tabulated3d.hpp
)

add_subdirectory(Python)
204 changes: 187 additions & 17 deletions src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ class PiecewisePolytropicFluid;
class Spectral;
template <typename LowDensityEoS>
class Enthalpy;
template <bool IsRelativistic>
class Tabulated3D;
} // namespace EquationsOfState
/// \endcond

Expand Down Expand Up @@ -67,6 +69,16 @@ struct DerivedClasses<false, 2> {
using type = tmpl::list<IdealFluid<false>, HybridEos<PolytropicFluid<false>>>;
};

template <>
struct DerivedClasses<true, 3> {
using type = tmpl::list<Tabulated3D<true>>;
};

template <>
struct DerivedClasses<false, 3> {
using type = tmpl::list<Tabulated3D<false>>;
};

} // namespace detail

/*!
Expand Down Expand Up @@ -392,6 +404,155 @@ class EquationOfState<IsRelativistic, 2> : public PUP::able {
/// The lower bound of the specific enthalpy that is valid for this EOS
virtual double specific_enthalpy_lower_bound() const = 0;
};

/*!
* \ingroup EquationsOfStateGroup
* \brief Base class for equations of state which need three independent
* thermodynamic variables in order to determine the pressure.
*
* The template parameter `IsRelativistic` is `true` for relativistic equations
* of state and `false` for non-relativistic equations of state.
*/
template <bool IsRelativistic>
class EquationOfState<IsRelativistic, 3> : public PUP::able {
public:
static constexpr bool is_relativistic = IsRelativistic;
static constexpr size_t thermodynamic_dim = 3;
using creatable_classes =
typename detail::DerivedClasses<IsRelativistic, 3>::type;

EquationOfState() = default;
EquationOfState(const EquationOfState&) = default;
EquationOfState& operator=(const EquationOfState&) = default;
EquationOfState(EquationOfState&&) = default;
EquationOfState& operator=(EquationOfState&&) = default;
~EquationOfState() override = default;

explicit EquationOfState(CkMigrateMessage* msg) : PUP::able(msg) {}

WRAPPED_PUPable_abstract(EquationOfState); // NOLINT

virtual inline std::unique_ptr<EquationOfState<IsRelativistic, 3>> get_clone()
const = 0;

virtual bool is_equal(
const EquationOfState<IsRelativistic, 3>& rhs) const = 0;
/// @{
/*!
* Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$, the
* specific internal energy \f$\epsilon\f$ and electron fraction \f$Y_e\f$.
*/
virtual Scalar<double> pressure_from_density_and_energy(
const Scalar<double>& /*rest_mass_density*/,
const Scalar<double>& /*specific_internal_energy*/,
const Scalar<double>& /*electron_fraction*/) const = 0;
virtual Scalar<DataVector> pressure_from_density_and_energy(
const Scalar<DataVector>& /*rest_mass_density*/,
const Scalar<DataVector>& /*specific_internal_energy*/,
const Scalar<DataVector>& /*electron_fraction*/) const = 0;
/// @}

/// @{
/*!
* Computes the pressure \f$p\f$ from the rest mass density \f$\rho\f$, the
* temperature \f$T\f$, and electron fraction \f$Y_e\f$.
*/
virtual Scalar<double> pressure_from_density_and_temperature(
const Scalar<double>& /*rest_mass_density*/,
const Scalar<double>& /*temperature*/,
const Scalar<double>& /*electron_fraction*/) const = 0;
virtual Scalar<DataVector> pressure_from_density_and_temperature(
const Scalar<DataVector>& /*rest_mass_density*/,
const Scalar<DataVector>& /*temperature*/,
const Scalar<DataVector>& /*electron_fraction*/) const = 0;
/// @}

/// @{
/*!
* Computes the temperature \f$T\f$ from the rest mass
* density \f$\rho\f$, the specific internal energy \f$\epsilon\f$,
* and electron fraction \f$Y_e\f$.
*/
virtual Scalar<double> temperature_from_density_and_energy(
const Scalar<double>& /*rest_mass_density*/,
const Scalar<double>& /*specific_internal_energy*/,
const Scalar<double>& /*electron_fraction*/) const = 0;
virtual Scalar<DataVector> temperature_from_density_and_energy(
const Scalar<DataVector>& /*rest_mass_density*/,
const Scalar<DataVector>& /*specific_internal_energy*/,
const Scalar<DataVector>& /*electron_fraction*/) const = 0;
/// @}

/// @{
/*!
* Computes the specific internal energy \f$\epsilon\f$ from the rest mass
* density \f$\rho\f$, the temperature \f$T\f$, and electron fraction
* \f$Y_e\f$.
*/
virtual Scalar<double> specific_internal_energy_from_density_and_temperature(
const Scalar<double>& /*rest_mass_density*/,
const Scalar<double>& /*temperature*/,
const Scalar<double>& /*electron_fraction*/
) const = 0;
virtual Scalar<DataVector>
specific_internal_energy_from_density_and_temperature(
const Scalar<DataVector>& /*rest_mass_density*/,
const Scalar<DataVector>& /*temperature*/,
const Scalar<DataVector>& /*electron_fraction*/
) const = 0;
/// @}

/// @{
/*!
* Computes sound speed squared \f$ c_s^2 \f$.
*
* Interpolate as a function of temperature, rest-mass density and electron
* fraction. Note that this will break thermodynamic consistency with the
* pressure and internal energy interpolated separately. The precise impact of
* this will depend on the EoS and numerical scheme used for the evolution.
*/
virtual Scalar<double> sound_speed_squared_from_density_and_temperature(
const Scalar<double>& /*rest_mass_density*/,
const Scalar<double>& /*temperature*/,
const Scalar<double>& /*electron_fraction*/) const = 0;
virtual Scalar<DataVector> sound_speed_squared_from_density_and_temperature(
const Scalar<DataVector>& /*rest_mass_density*/,
const Scalar<DataVector>& /*temperature*/,
const Scalar<DataVector>& /*electron_fraction*/) const = 0;
/// @}

/// The lower bound of the electron fraction that is valid for this EOS
virtual double electron_fraction_lower_bound() const = 0;

/// The upper bound of the electron fraction that is valid for this EOS
virtual double electron_fraction_upper_bound() const = 0;

/// The lower bound of the rest mass density that is valid for this EOS
virtual double rest_mass_density_lower_bound() const = 0;

/// The upper bound of the rest mass density that is valid for this EOS
virtual double rest_mass_density_upper_bound() const = 0;

/// The lower bound of the specific internal energy that is valid for this EOS
/// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$.
virtual double specific_internal_energy_lower_bound(
const double rest_mass_density, const double electron_fraction) const = 0;

/// The upper bound of the specific internal energy that is valid for this EOS
/// at the given rest mass density \f$\rho\f$ and electron fraction \f$Y_e\f$.
virtual double specific_internal_energy_upper_bound(
const double rest_mass_density, const double electron_fraction) const = 0;

/// The lower bound of the specific enthalpy that is valid for this EOS
virtual double specific_enthalpy_lower_bound() const = 0;

/// The lower bound of the temperature that is valid for this EOS
virtual double temperature_lower_bound() const = 0;

/// The upper bound of the temperature that is valid for this EOS
virtual double temperature_upper_bound() const = 0;
};

/// Compare two equations of state for equality
template <bool IsRelLhs, bool IsRelRhs, size_t ThermoDimLhs,
size_t ThermoDimRhs>
Expand Down Expand Up @@ -425,6 +586,12 @@ bool operator!=(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
chi_from_density_and_energy, \
kappa_times_p_over_rho_squared_from_density_and_energy)

#define EQUATION_OF_STATE_FUNCTIONS_3D \
(pressure_from_density_and_energy, pressure_from_density_and_temperature, \
temperature_from_density_and_energy, \
specific_internal_energy_from_density_and_temperature, \
sound_speed_squared_from_density_and_temperature)

#define EQUATION_OF_STATE_ARGUMENTS_EXPAND(z, n, type) \
BOOST_PP_COMMA_IF(n) const Scalar<type>&

Expand All @@ -442,16 +609,17 @@ bool operator!=(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
* \brief Macro used to generate forward declarations of member functions in
* derived classes
*/
#define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(DERIVED, DIM) \
BOOST_PP_LIST_FOR_EACH( \
EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS_HELPER, DIM, \
BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
BOOST_PP_SUB(DIM, 1), \
(EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D)))) \
\
/* clang-tidy: do not use non-const references */ \
void pup(PUP::er& p) override; /* NOLINT */ \
\
#define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS(DERIVED, DIM) \
BOOST_PP_LIST_FOR_EACH( \
EQUATION_OF_STATE_FORWARD_DECLARE_MEMBERS_HELPER, DIM, \
BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
BOOST_PP_SUB(DIM, 1), \
(EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
EQUATION_OF_STATE_FUNCTIONS_3D)))) \
\
/* clang-tidy: do not use non-const references */ \
void pup(PUP::er& p) override; /* NOLINT */ \
\
explicit DERIVED(CkMigrateMessage* msg);

/// \cond
Expand Down Expand Up @@ -484,7 +652,8 @@ bool operator!=(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
(TEMPLATE, DERIVED, DATA_TYPE, DIM), \
BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
BOOST_PP_SUB(DIM, 1), \
(EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D))))
(EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
EQUATION_OF_STATE_FUNCTIONS_3D))))

/// \cond
#define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER(r, DIM, \
Expand All @@ -494,9 +663,10 @@ bool operator!=(const EquationOfState<IsRelLhs, ThermoDimLhs>& lhs,
DIM, EQUATION_OF_STATE_ARGUMENTS_EXPAND, DataType)) const;
/// \endcond

#define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(DIM) \
BOOST_PP_LIST_FOR_EACH( \
EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER, DIM, \
BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
BOOST_PP_SUB(DIM, 1), \
(EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D))))
#define EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS(DIM) \
BOOST_PP_LIST_FOR_EACH( \
EQUATION_OF_STATE_FORWARD_DECLARE_MEMBER_IMPLS_HELPER, DIM, \
BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TUPLE_ELEM( \
BOOST_PP_SUB(DIM, 1), \
(EQUATION_OF_STATE_FUNCTIONS_1D, EQUATION_OF_STATE_FUNCTIONS_2D, \
EQUATION_OF_STATE_FUNCTIONS_3D))))
Loading

0 comments on commit 4dd2ab1

Please sign in to comment.