4237779830

Build Type

push

github

Committed by GitHub

Commit Message

Merge pull request #4755 from nilsvu/test_bbh_domain

Run Details

63865 of 66590 relevant lines covered (95.91%)

411287.72 hits per line

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

98.62

/src/PointwiseFunctions/AnalyticSolutions/NewtonianEuler/RiemannProblem.cpp

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

#include "PointwiseFunctions/AnalyticSolutions/NewtonianEuler/RiemannProblem.hpp"

#include <cmath>
#include <cstddef>
#include <limits>
#include <pup.h>

#include "DataStructures/DataVector.hpp"
#include "DataStructures/Tensor/Tensor.hpp"
#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
#include "Utilities/ContainerHelpers.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeWithValue.hpp"
#include "Utilities/Math.hpp"

namespace {

// Any of the two functions of the pressure and the initial data
// in Eqns. (4.6) and (4.7) of Toro.
template <size_t Dim>
struct FunctionOfPressureAndData {
  FunctionOfPressureAndData(const typename NewtonianEuler::Solutions::
                                RiemannProblem<Dim>::InitialData& data,
                            const double adiabatic_index)
      : state_pressure_(data.pressure_),
        constant_a_(data.constant_a_),
        constant_b_(data.constant_b_) {
    prefactor_ = 2.0 * data.sound_speed_ / (adiabatic_index - 1.0);
    exponent_ = 0.5 * (adiabatic_index - 1.0) / adiabatic_index;
  }

  double operator()(const double pressure) const {
    // Value depends on whether the initial state is a shock
    // (pressure > pressure of initial state) or a rarefaction wave
    // (pressure <= pressure of initial state)
    return pressure > state_pressure_
               ? (pressure - state_pressure_) *
                     sqrt(constant_a_ / (pressure + constant_b_))
               : prefactor_ *
                     (pow(pressure / state_pressure_, exponent_) - 1.0);
  }

 private:
  double state_pressure_ = std::numeric_limits<double>::signaling_NaN();
  double constant_a_ = std::numeric_limits<double>::signaling_NaN();
  double constant_b_ = std::numeric_limits<double>::signaling_NaN();

  // Auxiliary variables for computing the rarefaction wave
  double prefactor_ = std::numeric_limits<double>::signaling_NaN();
  double exponent_ = std::numeric_limits<double>::signaling_NaN();
};

}  // namespace

namespace NewtonianEuler::Solutions {

template <size_t Dim>
RiemannProblem<Dim>::RiemannProblem(
    const double adiabatic_index, const double initial_position,
    const double left_mass_density,
    const std::array<double, Dim>& left_velocity, const double left_pressure,
    const double right_mass_density,
    const std::array<double, Dim>& right_velocity, const double right_pressure,
    const double pressure_star_tol)
    : adiabatic_index_(adiabatic_index),
      initial_position_(initial_position),
      left_initial_data_(left_mass_density, left_velocity, left_pressure,
                         adiabatic_index, propagation_axis_),
      right_initial_data_(right_mass_density, right_velocity, right_pressure,
                          adiabatic_index, propagation_axis_),
      pressure_star_tol_(pressure_star_tol),
      equation_of_state_(adiabatic_index) {
  const double delta_u = right_initial_data_.normal_velocity_ -
                         left_initial_data_.normal_velocity_;

  // The pressure positivity condition must be met (Eqn. 4.40 of Toro).
  ASSERT(2.0 *
                     (left_initial_data_.sound_speed_ +
                      right_initial_data_.sound_speed_) /
                     (adiabatic_index_ - 1.0) -
                 delta_u >
             0.0,
         "The pressure positivity condition must be met. Initial data do not "
         "satisfy this criterion.");

  // Before evaluating the solution at any (x, t), the type of solution
  // (shock or rarefaction) on each side of the contact discontinuity
  // must be sorted out. To do so, the first step is to compute the
  // state variables in the star region by solving a transcendental equation,
  // which is what all the math in this constructor is about.

  // Compute bracket for root finder according to value of the function whose
  // root we want (Eqn. 4.39 of Toro.)
  const FunctionOfPressureAndData<Dim> f_of_p_left(left_initial_data_,
                                                   adiabatic_index_);
  const FunctionOfPressureAndData<Dim> f_of_p_right(right_initial_data_,
                                                    adiabatic_index_);
  const auto p_minmax =
      std::minmax(left_initial_data_.pressure_, right_initial_data_.pressure_);
  const double f_min =
      f_of_p_left(p_minmax.first) + f_of_p_right(p_minmax.first) + delta_u;
  const double f_max =
      f_of_p_left(p_minmax.second) + f_of_p_right(p_minmax.second) + delta_u;

  if (f_min > 0.0 and f_max > 0.0) {
    const double exponent = 0.5 * (adiabatic_index_ - 1.0) / adiabatic_index_;
    pressure_star_ = std::pow(
        (left_initial_data_.sound_speed_ + right_initial_data_.sound_speed_ -
         0.5 * (adiabatic_index_ - 1.0) * delta_u) /
            (left_initial_data_.sound_speed_ *
                 std::pow(left_initial_data_.pressure_, -exponent) +
             right_initial_data_.sound_speed_ *
                 std::pow(right_initial_data_.pressure_, -exponent)),
        1.0 / exponent);
    ASSERT(std::abs(f_of_p_left(pressure_star_) + f_of_p_right(pressure_star_) +
                    delta_u) < 1.0e-8,
           "Failed to analytically solve correctly.");
  } else {
    double pressure_lower = std::numeric_limits<double>::signaling_NaN();
    double pressure_upper = std::numeric_limits<double>::signaling_NaN();
    if (f_min < 0.0 and f_max > 0.0) {
      pressure_lower = p_minmax.first;
      pressure_upper = p_minmax.second;
    } else {
      pressure_lower = p_minmax.second;
      pressure_upper = 10.0 * pressure_lower;  // Arbitrary upper bound < \infty
    }

    // Now get pressure by solving transcendental equation.
    const auto f_of_p = [&f_of_p_left, &f_of_p_right,
                         &delta_u](const double pressure) {
      // Function of pressure in Eqn. (4.5) of Toro.
      return f_of_p_left(pressure) + f_of_p_right(pressure) + delta_u;
    };
    try {
      pressure_star_ = RootFinder::toms748(
          f_of_p, pressure_lower, pressure_upper, pressure_star_tol_, 1.0e-15);
    } catch (std::exception& exception) {
      ERROR(
          "Failed to find p_* with Newton-Raphson root finder. Got "
          "exception message:\n"
          << exception.what()
          << "\nIf the residual is small you can change the tolerance for the "
             "root finder in the input file.");
    }
  }

  // Calculated p_*, u_* is obtained from Eqn. (4.9) in Toro.
  velocity_star_ =
      0.5 * (left_initial_data_.normal_velocity_ +
             right_initial_data_.normal_velocity_ -
             f_of_p_left(pressure_star_) + f_of_p_right(pressure_star_));
}

template <size_t Dim>
void RiemannProblem<Dim>::pup(PUP::er& p) {
  p | adiabatic_index_;
  p | initial_position_;
  p | left_initial_data_;
  p | right_initial_data_;
  p | pressure_star_tol_;
  p | pressure_star_;
  p | velocity_star_;
  p | equation_of_state_;
}

template <size_t Dim>
RiemannProblem<Dim>::InitialData::InitialData(
    const double mass_density, const std::array<double, Dim>& velocity,
    const double pressure, const double adiabatic_index,
    const size_t propagation_axis)
    : mass_density_(mass_density), velocity_(velocity), pressure_(pressure) {
  ASSERT(mass_density_ > 0.0,
         "The mass density must be positive. Value given: " << mass_density);
  ASSERT(pressure_ > 0.0,
         "The pressure must be positive. Value given: " << pressure);

  sound_speed_ = sqrt(adiabatic_index * pressure / mass_density);
  normal_velocity_ = gsl::at(velocity, propagation_axis);
  constant_a_ = 2.0 / (adiabatic_index + 1.0) / mass_density;
  constant_b_ = (adiabatic_index - 1.0) * pressure / (adiabatic_index + 1.0);
}

template <size_t Dim>
void RiemannProblem<Dim>::InitialData::pup(PUP::er& p) {
  p | mass_density_;
  p | velocity_;
  p | pressure_;
  p | sound_speed_;
  p | normal_velocity_;
  p | constant_a_;
  p | constant_b_;
}

template <size_t Dim>
RiemannProblem<Dim>::Wave::Wave(const InitialData& data,
                                const double pressure_star,
                                const double velocity_star,
                                const double adiabatic_index, const Side& side)
    : pressure_ratio_(pressure_star / data.pressure_),
      is_shock_(pressure_ratio_ > 1.0),
      data_(data),
      shock_(data, pressure_ratio_, adiabatic_index, side),
      rarefaction_(data, pressure_ratio_, velocity_star, adiabatic_index,
                   side) {}

template <size_t Dim>
double RiemannProblem<Dim>::Wave::mass_density(const double x_shifted,
                                               const double t) const {
  return (is_shock_ ? shock_.mass_density(x_shifted, t, data_)
                    : rarefaction_.mass_density(x_shifted, t, data_));
}

template <size_t Dim>
double RiemannProblem<Dim>::Wave::normal_velocity(
    const double x_shifted, const double t, const double velocity_star) const {
  return (is_shock_ ? shock_.normal_velocity(x_shifted, t, data_, velocity_star)
                    : rarefaction_.normal_velocity(x_shifted, t, data_,
                                                   velocity_star));
}

template <size_t Dim>
double RiemannProblem<Dim>::Wave::pressure(const double x_shifted,
                                           const double t,
                                           const double pressure_star) const {
  return (is_shock_
              ? shock_.pressure(x_shifted, t, data_, pressure_star)
              : rarefaction_.pressure(x_shifted, t, data_, pressure_star));
}

template <size_t Dim>
RiemannProblem<Dim>::Shock::Shock(const InitialData& data,
                                  const double pressure_ratio,
                                  const double adiabatic_index,
                                  const Side& side)
    : direction_(side == Side::Left ? -1.0 : 1.0) {
  const double gamma_mm = adiabatic_index - 1.0;
  const double gamma_pp = adiabatic_index + 1.0;
  const double gamma_mm_over_gamma_pp = gamma_mm / gamma_pp;

  mass_density_star_ = data.mass_density_ *
                       (pressure_ratio + gamma_mm_over_gamma_pp) /
                       (pressure_ratio * gamma_mm_over_gamma_pp + 1.0);

  shock_speed_ =
      data.normal_velocity_ +
      direction_ * data.sound_speed_ *
          sqrt(0.5 * (gamma_pp * pressure_ratio + gamma_mm) / adiabatic_index);
}

template <size_t Dim>
double RiemannProblem<Dim>::Shock::mass_density(const double x_shifted,
                                                const double t,
                                                const InitialData& data) const {
  return mass_density_star_ +
         (data.mass_density_ - mass_density_star_) *
             step_function(direction_ * (x_shifted - shock_speed_ * t));
}

template <size_t Dim>
double RiemannProblem<Dim>::Shock::normal_velocity(
    const double x_shifted, const double t, const InitialData& data,
    const double velocity_star) const {
  return velocity_star +
         (data.normal_velocity_ - velocity_star) *
             step_function(direction_ * (x_shifted - shock_speed_ * t));
}

template <size_t Dim>
double RiemannProblem<Dim>::Shock::pressure(const double x_shifted,
                                            const double t,
                                            const InitialData& data,
                                            const double pressure_star) const {
  return pressure_star +
         (data.pressure_ - pressure_star) *
             step_function(direction_ * (x_shifted - shock_speed_ * t));
}

template <size_t Dim>
RiemannProblem<Dim>::Rarefaction::Rarefaction(const InitialData& data,
                                              const double pressure_ratio,
                                              const double velocity_star,
                                              const double adiabatic_index,
                                              const Side& side)
    : direction_(side == Side::Left ? -1.0 : 1.0) {
  gamma_mm_ = adiabatic_index - 1.0;
  gamma_pp_ = adiabatic_index + 1.0;
  mass_density_star_ =
      data.mass_density_ * pow(pressure_ratio, 1.0 / adiabatic_index);
  sound_speed_star_ =
      data.sound_speed_ *
      pow(pressure_ratio, 0.5 * (adiabatic_index - 1.0) / adiabatic_index);
  head_speed_ = data.normal_velocity_ + direction_ * data.sound_speed_;
  tail_speed_ = velocity_star + direction_ * sound_speed_star_;
}

template <size_t Dim>
double RiemannProblem<Dim>::Rarefaction::mass_density(
    const double x_shifted, const double t, const InitialData& data) const {
  const double s = (t > 0.0 ? (x_shifted / t) : 0.0);
  return direction_ * (x_shifted - tail_speed_ * t) < 0.0
             ? mass_density_star_
             : (direction_ * (x_shifted - tail_speed_ * t) >= 0.0 and
                        direction_ * (x_shifted - head_speed_ * t) < 0.0
                    ? (data.mass_density_ *
                       pow((2.0 - direction_ * gamma_mm_ *
                                      (data.normal_velocity_ - s) /
                                      data.sound_speed_) /
                               gamma_pp_,
                           2.0 / gamma_mm_))
                    : data.mass_density_);
}

template <size_t Dim>
double RiemannProblem<Dim>::Rarefaction::normal_velocity(
    const double x_shifted, const double t, const InitialData& data,
    const double velocity_star) const {
  const double s = (t > 0.0 ? (x_shifted / t) : 0.0);
  return direction_ * (x_shifted - tail_speed_ * t) < 0.0
             ? velocity_star
             : (direction_ * (x_shifted - tail_speed_ * t) >= 0.0 and
                        direction_ * (x_shifted - head_speed_ * t) < 0.0
                    ? (2.0 *
                       (0.5 * gamma_mm_ * data.normal_velocity_ + s -
                        direction_ * data.sound_speed_) /
                       gamma_pp_)
                    : data.normal_velocity_);
}

template <size_t Dim>
double RiemannProblem<Dim>::Rarefaction::pressure(
    const double x_shifted, const double t, const InitialData& data,
    const double pressure_star) const {
  const double s = (t > 0.0 ? (x_shifted / t) : 0.0);
  return direction_ * (x_shifted - tail_speed_ * t) < 0.0
             ? pressure_star
             : (direction_ * (x_shifted - tail_speed_ * t) >= 0.0 and
                        direction_ * (x_shifted - head_speed_ * t) < 0.0
                    ? (data.pressure_ *
                       pow((2.0 - direction_ * gamma_mm_ *
                                      (data.normal_velocity_ - s) /
                                      data.sound_speed_) /
                               gamma_pp_,
                           2.0 * (gamma_mm_ + 1.0) / gamma_mm_))
                    : data.pressure_);
}

template <size_t Dim>
template <typename DataType>
tuples::TaggedTuple<Tags::MassDensity<DataType>> RiemannProblem<Dim>::variables(
    const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
    tmpl::list<Tags::MassDensity<DataType>> /*meta*/, const Wave& left,
    const Wave& right) const {
  auto mass_density = make_with_value<Scalar<DataType>>(x_shifted, 0.0);
  const double u_star_times_t = velocity_star_ * t;
  for (size_t s = 0; s < get_size(get<0>(x_shifted)); ++s) {
    const double x_shifted_s = get_element(x_shifted.get(propagation_axis_), s);
    get_element(get(mass_density), s) =
        (x_shifted_s < u_star_times_t ? left.mass_density(x_shifted_s, t)
                                      : right.mass_density(x_shifted_s, t));
  }
  return mass_density;
}

template <size_t Dim>
template <typename DataType>
tuples::TaggedTuple<Tags::Velocity<DataType, Dim>>
RiemannProblem<Dim>::variables(
    const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
    tmpl::list<Tags::Velocity<DataType, Dim>> /*meta*/, const Wave& left,
    const Wave& right) const {
  auto velocity = make_with_value<tnsr::I<DataType, Dim, Frame::Inertial>>(
      get<0>(x_shifted), 0.0);

  const double u_star_times_t = velocity_star_ * t;
  for (size_t s = 0; s < get_size(get<0>(x_shifted)); ++s) {
    const double x_shifted_s = get_element(x_shifted.get(propagation_axis_), s);

    size_t index = propagation_axis_ % Dim;
    get_element(velocity.get(index), s) =
        (x_shifted_s < u_star_times_t
             ? left.normal_velocity(x_shifted_s, t, velocity_star_)
             : right.normal_velocity(x_shifted_s, t, velocity_star_));

    if (Dim > 1) {
      index = (propagation_axis_ + 1) % Dim;
      get_element(velocity.get(index), s) =
          (x_shifted_s < u_star_times_t
               ? gsl::at(left_initial_data_.velocity_, index)
               : gsl::at(right_initial_data_.velocity_, index));
    }

    if (Dim > 2) {
      index = (propagation_axis_ + 2) % Dim;
      get_element(velocity.get(index), s) =
          (x_shifted_s < u_star_times_t
               ? gsl::at(left_initial_data_.velocity_, index)
               : gsl::at(right_initial_data_.velocity_, index));
    }
  }
  return velocity;
}

template <size_t Dim>
template <typename DataType>
tuples::TaggedTuple<Tags::Pressure<DataType>> RiemannProblem<Dim>::variables(
    const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
    tmpl::list<Tags::Pressure<DataType>> /*meta*/, const Wave& left,
    const Wave& right) const {
  auto pressure = make_with_value<Scalar<DataType>>(x_shifted, 0.0);
  const double u_star_times_t = velocity_star_ * t;
  for (size_t s = 0; s < get_size(get<0>(x_shifted)); ++s) {
    const double x_shifted_s = get_element(x_shifted.get(propagation_axis_), s);
    get_element(get(pressure), s) =
        (x_shifted_s < u_star_times_t
             ? left.pressure(x_shifted_s, t, pressure_star_)
             : right.pressure(x_shifted_s, t, pressure_star_));
  }
  return pressure;
}

template <size_t Dim>
template <typename DataType>
tuples::TaggedTuple<Tags::SpecificInternalEnergy<DataType>>
RiemannProblem<Dim>::variables(
    const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
    tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/,
    const Wave& left, const Wave& right) const {
  return equation_of_state_.specific_internal_energy_from_density_and_pressure(
      get<Tags::MassDensity<DataType>>(
          variables(x_shifted, t, tmpl::list<Tags::MassDensity<DataType>>{},
                    left, right)),
      get<Tags::Pressure<DataType>>(variables(
          x_shifted, t, tmpl::list<Tags::Pressure<DataType>>{}, left, right)));
}

template <size_t Dim>
bool operator==(const RiemannProblem<Dim>& lhs,
                const RiemannProblem<Dim>& rhs) {
  return lhs.adiabatic_index_ == rhs.adiabatic_index_ and
         lhs.initial_position_ == rhs.initial_position_ and
         lhs.left_initial_data_ == rhs.left_initial_data_ and
         lhs.right_initial_data_ == rhs.right_initial_data_ and
         lhs.pressure_star_tol_ == rhs.pressure_star_tol_ and
         lhs.pressure_star_ == rhs.pressure_star_ and
         lhs.velocity_star_ == rhs.velocity_star_;
}

template <size_t Dim>
bool operator!=(const RiemannProblem<Dim>& lhs,
                const RiemannProblem<Dim>& rhs) {
  return not(lhs == rhs);
}

#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(1, data)
#define TAG(data) BOOST_PP_TUPLE_ELEM(2, data)

#define INSTANTIATE_CLASS(_, data)                            \
  template class RiemannProblem<DIM(data)>;                   \
  template bool operator==(const RiemannProblem<DIM(data)>&,  \
                           const RiemannProblem<DIM(data)>&); \
  template bool operator!=(const RiemannProblem<DIM(data)>&,  \
                           const RiemannProblem<DIM(data)>&);

GENERATE_INSTANTIATIONS(INSTANTIATE_CLASS, (1, 2, 3))

#define INSTANTIATE_SCALARS(_, data)                                         \
  template tuples::TaggedTuple<TAG(data) < DTYPE(data)> >                    \
      RiemannProblem<DIM(data)>::variables(                                  \
          const tnsr::I<DTYPE(data), DIM(data), Frame::Inertial>& x_shifted, \
          const double t, tmpl::list<TAG(data) < DTYPE(data)> >,             \
          const Wave& left, const Wave& right) const;

GENERATE_INSTANTIATIONS(INSTANTIATE_SCALARS, (1, 2, 3), (double, DataVector),
                        (Tags::MassDensity, Tags::Pressure,
                         Tags::SpecificInternalEnergy))

#define INSTANTIATE_VELOCITY(_, data)                                        \
  template tuples::TaggedTuple<TAG(data) < DTYPE(data), DIM(data)> >         \
      RiemannProblem<DIM(data)>::variables(                                  \
          const tnsr::I<DTYPE(data), DIM(data), Frame::Inertial>& x_shifted, \
          const double t, tmpl::list<TAG(data) < DTYPE(data), DIM(data)> >,  \
          const Wave& left, const Wave& right) const;

GENERATE_INSTANTIATIONS(INSTANTIATE_VELOCITY, (1, 2, 3), (double, DataVector),
                        (Tags::Velocity))

#undef DIM
#undef DTYPE
#undef TAG
#undef INSTANTIATE_CLASS
#undef INSTANTIATE_SCALARS
#undef INSTANTIATE_VELOCITY

}  // namespace NewtonianEuler::Solutions

1	// Distributed under the MIT License.
2	// See LICENSE.txt for details.
3
4	#include "PointwiseFunctions/AnalyticSolutions/NewtonianEuler/RiemannProblem.hpp"
5
6	#include <cmath>
7	#include <cstddef>
8	#include <limits>
9	#include <pup.h>
10
11	#include "DataStructures/DataVector.hpp"
12	#include "DataStructures/Tensor/Tensor.hpp"
13	#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
14	#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
15	#include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
16	#include "Utilities/ContainerHelpers.hpp"
17	#include "Utilities/GenerateInstantiations.hpp"
18	#include "Utilities/Gsl.hpp"
19	#include "Utilities/MakeWithValue.hpp"
20	#include "Utilities/Math.hpp"
21
22	namespace {
23
24	// Any of the two functions of the pressure and the initial data
25	// in Eqns. (4.6) and (4.7) of Toro.
26	template <size_t Dim>
27	struct FunctionOfPressureAndData {
28	FunctionOfPressureAndData(const typename NewtonianEuler::Solutions::	46✔
29	RiemannProblem<Dim>::InitialData& data,
30	const double adiabatic_index)
31	: state_pressure_(data.pressure_),	46✔
32	constant_a_(data.constant_a_),	46✔
33	constant_b_(data.constant_b_) {	46✔
34	prefactor_ = 2.0 * data.sound_speed_ / (adiabatic_index - 1.0);	46✔
35	exponent_ = 0.5 * (adiabatic_index - 1.0) / adiabatic_index;	46✔
36	}	46✔
37
38	double operator()(const double pressure) const {	578✔
39	// Value depends on whether the initial state is a shock
40	// (pressure > pressure of initial state) or a rarefaction wave
41	// (pressure <= pressure of initial state)
42	return pressure > state_pressure_	578✔
43	? (pressure - state_pressure_) *	250✔
44	sqrt(constant_a_ / (pressure + constant_b_))	250✔
45	: prefactor_ *	328✔
46	(pow(pressure / state_pressure_, exponent_) - 1.0);	578✔
47	}
48
49	private:
50	double state_pressure_ = std::numeric_limits<double>::signaling_NaN();
51	double constant_a_ = std::numeric_limits<double>::signaling_NaN();
52	double constant_b_ = std::numeric_limits<double>::signaling_NaN();
53
54	// Auxiliary variables for computing the rarefaction wave
55	double prefactor_ = std::numeric_limits<double>::signaling_NaN();
56	double exponent_ = std::numeric_limits<double>::signaling_NaN();
57	};
58
59	} // namespace
60
61	namespace NewtonianEuler::Solutions {
62
63	template <size_t Dim>
64	RiemannProblem<Dim>::RiemannProblem(	26✔
65	const double adiabatic_index, const double initial_position,
66	const double left_mass_density,
67	const std::array<double, Dim>& left_velocity, const double left_pressure,
68	const double right_mass_density,
69	const std::array<double, Dim>& right_velocity, const double right_pressure,
70	const double pressure_star_tol)
71	: adiabatic_index_(adiabatic_index),
72	initial_position_(initial_position),
73	left_initial_data_(left_mass_density, left_velocity, left_pressure,
74	adiabatic_index, propagation_axis_),
75	right_initial_data_(right_mass_density, right_velocity, right_pressure,
76	adiabatic_index, propagation_axis_),
77	pressure_star_tol_(pressure_star_tol),
78	equation_of_state_(adiabatic_index) {	27✔
79	const double delta_u = right_initial_data_.normal_velocity_ -	24✔
80	left_initial_data_.normal_velocity_;	24✔
81
82	// The pressure positivity condition must be met (Eqn. 4.40 of Toro).
83	ASSERT(2.0 *	25✔
84	(left_initial_data_.sound_speed_ +
85	right_initial_data_.sound_speed_) /
86	(adiabatic_index_ - 1.0) -
87	delta_u >
88	0.0,
89	"The pressure positivity condition must be met. Initial data do not "
90	"satisfy this criterion.");
91
92	// Before evaluating the solution at any (x, t), the type of solution
93	// (shock or rarefaction) on each side of the contact discontinuity
94	// must be sorted out. To do so, the first step is to compute the
95	// state variables in the star region by solving a transcendental equation,
96	// which is what all the math in this constructor is about.
97
98	// Compute bracket for root finder according to value of the function whose
99	// root we want (Eqn. 4.39 of Toro.)
100	const FunctionOfPressureAndData<Dim> f_of_p_left(left_initial_data_,	23✔
101	adiabatic_index_);
102	const FunctionOfPressureAndData<Dim> f_of_p_right(right_initial_data_,	23✔
103	adiabatic_index_);
104	const auto p_minmax =
105	std::minmax(left_initial_data_.pressure_, right_initial_data_.pressure_);	23✔
106	const double f_min =	23✔
107	f_of_p_left(p_minmax.first) + f_of_p_right(p_minmax.first) + delta_u;	23✔
108	const double f_max =	23✔
109	f_of_p_left(p_minmax.second) + f_of_p_right(p_minmax.second) + delta_u;	23✔
110
111	if (f_min > 0.0 and f_max > 0.0) {	23✔
112	const double exponent = 0.5 * (adiabatic_index_ - 1.0) / adiabatic_index_;	1✔
113	pressure_star_ = std::pow(	1✔
114	(left_initial_data_.sound_speed_ + right_initial_data_.sound_speed_ -	1✔
115	0.5 * (adiabatic_index_ - 1.0) * delta_u) /	1✔
116	(left_initial_data_.sound_speed_ *	1✔
117	std::pow(left_initial_data_.pressure_, -exponent) +	1✔
118	right_initial_data_.sound_speed_ *	1✔
119	std::pow(right_initial_data_.pressure_, -exponent)),	1✔
120	1.0 / exponent);
121	ASSERT(std::abs(f_of_p_left(pressure_star_) + f_of_p_right(pressure_star_) +	1✔
122	delta_u) < 1.0e-8,
123	"Failed to analytically solve correctly.");	1✔
124	} else {
125	double pressure_lower = std::numeric_limits<double>::signaling_NaN();	22✔
126	double pressure_upper = std::numeric_limits<double>::signaling_NaN();	22✔
127	if (f_min < 0.0 and f_max > 0.0) {	22✔
128	pressure_lower = p_minmax.first;	21✔
129	pressure_upper = p_minmax.second;	21✔
130	} else {
131	pressure_lower = p_minmax.second;	1✔
132	pressure_upper = 10.0 * pressure_lower; // Arbitrary upper bound < \infty	1✔
133	}
134
135	// Now get pressure by solving transcendental equation.
136	const auto f_of_p = [&f_of_p_left, &f_of_p_right,	679✔
137	&delta_u](const double pressure) {
138	// Function of pressure in Eqn. (4.5) of Toro.
139	return f_of_p_left(pressure) + f_of_p_right(pressure) + delta_u;	219✔
140	};
141	try {
142	pressure_star_ = RootFinder::toms748(	22✔
143	f_of_p, pressure_lower, pressure_upper, pressure_star_tol_, 1.0e-15);
144	} catch (std::exception& exception) {	×
145	ERROR(	×
146	"Failed to find p_* with Newton-Raphson root finder. Got "
147	"exception message:\n"
148	<< exception.what()
149	<< "\nIf the residual is small you can change the tolerance for the "
150	"root finder in the input file.");
151	}
152	}
153
154	// Calculated p_, u_ is obtained from Eqn. (4.9) in Toro.
155	velocity_star_ =	23✔
156	0.5 * (left_initial_data_.normal_velocity_ +	46✔
157	right_initial_data_.normal_velocity_ -	46✔
158	f_of_p_left(pressure_star_) + f_of_p_right(pressure_star_));	23✔
159	}	23✔
160
161	template <size_t Dim>
162	void RiemannProblem<Dim>::pup(PUP::er& p) {	18✔
163	p \| adiabatic_index_;	18✔
164	p \| initial_position_;	18✔
165	p \| left_initial_data_;	18✔
166	p \| right_initial_data_;	18✔
167	p \| pressure_star_tol_;	18✔
168	p \| pressure_star_;	18✔
169	p \| velocity_star_;	18✔
170	p \| equation_of_state_;	18✔
171	}	18✔
172
173	template <size_t Dim>
174	RiemannProblem<Dim>::InitialData::InitialData(	50✔
175	const double mass_density, const std::array<double, Dim>& velocity,
176	const double pressure, const double adiabatic_index,
177	const size_t propagation_axis)
178	: mass_density_(mass_density), velocity_(velocity), pressure_(pressure) {	50✔
179	ASSERT(mass_density_ > 0.0,	51✔
180	"The mass density must be positive. Value given: " << mass_density);
181	ASSERT(pressure_ > 0.0,	50✔
182	"The pressure must be positive. Value given: " << pressure);
183
184	sound_speed_ = sqrt(adiabatic_index * pressure / mass_density);	48✔
185	normal_velocity_ = gsl::at(velocity, propagation_axis);	48✔
186	constant_a_ = 2.0 / (adiabatic_index + 1.0) / mass_density;	48✔
187	constant_b_ = (adiabatic_index - 1.0) * pressure / (adiabatic_index + 1.0);	48✔
188	}	48✔
189
190	template <size_t Dim>
191	void RiemannProblem<Dim>::InitialData::pup(PUP::er& p) {	36✔
192	p \| mass_density_;	36✔
193	p \| velocity_;	36✔
194	p \| pressure_;	36✔
195	p \| sound_speed_;	36✔
196	p \| normal_velocity_;	36✔
197	p \| constant_a_;	36✔
198	p \| constant_b_;	36✔
199	}	36✔
200
201	template <size_t Dim>
202	RiemannProblem<Dim>::Wave::Wave(const InitialData& data,	48✔
203	const double pressure_star,
204	const double velocity_star,
205	const double adiabatic_index, const Side& side)
206	: pressure_ratio_(pressure_star / data.pressure_),	48✔
207	is_shock_(pressure_ratio_ > 1.0),	48✔
208	data_(data),
209	shock_(data, pressure_ratio_, adiabatic_index, side),
210	rarefaction_(data, pressure_ratio_, velocity_star, adiabatic_index,
211	side) {}	48✔
212
213	template <size_t Dim>
214	double RiemannProblem<Dim>::Wave::mass_density(const double x_shifted,	144✔
215	const double t) const {
216	return (is_shock_ ? shock_.mass_density(x_shifted, t, data_)	280✔
217	: rarefaction_.mass_density(x_shifted, t, data_));	280✔
218	}
219
220	template <size_t Dim>
221	double RiemannProblem<Dim>::Wave::normal_velocity(	72✔
222	const double x_shifted, const double t, const double velocity_star) const {
223	return (is_shock_ ? shock_.normal_velocity(x_shifted, t, data_, velocity_star)	140✔
224	: rarefaction_.normal_velocity(x_shifted, t, data_,	68✔
225	velocity_star));	72✔
226	}
227
228	template <size_t Dim>
229	double RiemannProblem<Dim>::Wave::pressure(const double x_shifted,	144✔
230	const double t,
231	const double pressure_star) const {
232	return (is_shock_	144✔
233	? shock_.pressure(x_shifted, t, data_, pressure_star)	280✔
234	: rarefaction_.pressure(x_shifted, t, data_, pressure_star));	280✔
235	}
236
237	template <size_t Dim>
238	RiemannProblem<Dim>::Shock::Shock(const InitialData& data,	48✔
239	const double pressure_ratio,
240	const double adiabatic_index,
241	const Side& side)
242	: direction_(side == Side::Left ? -1.0 : 1.0) {	48✔
243	const double gamma_mm = adiabatic_index - 1.0;	48✔
244	const double gamma_pp = adiabatic_index + 1.0;	48✔
245	const double gamma_mm_over_gamma_pp = gamma_mm / gamma_pp;	48✔
246
247	mass_density_star_ = data.mass_density_ *	48✔
248	(pressure_ratio + gamma_mm_over_gamma_pp) /	48✔
249	(pressure_ratio * gamma_mm_over_gamma_pp + 1.0);	48✔
250
251	shock_speed_ =	48✔
252	data.normal_velocity_ +	48✔
253	direction_ * data.sound_speed_ *	48✔
254	sqrt(0.5 * (gamma_pp * pressure_ratio + gamma_mm) / adiabatic_index);	48✔
255	}	48✔
256
257	template <size_t Dim>
258	double RiemannProblem<Dim>::Shock::mass_density(const double x_shifted,	8✔
259	const double t,
260	const InitialData& data) const {
261	return mass_density_star_ +	8✔
262	(data.mass_density_ - mass_density_star_) *	16✔
263	step_function(direction_ * (x_shifted - shock_speed_ * t));	8✔
264	}
265
266	template <size_t Dim>
267	double RiemannProblem<Dim>::Shock::normal_velocity(	4✔
268	const double x_shifted, const double t, const InitialData& data,
269	const double velocity_star) const {
270	return velocity_star +
271	(data.normal_velocity_ - velocity_star) *	8✔
272	step_function(direction_ * (x_shifted - shock_speed_ * t));	4✔
273	}
274
275	template <size_t Dim>
276	double RiemannProblem<Dim>::Shock::pressure(const double x_shifted,	8✔
277	const double t,
278	const InitialData& data,
279	const double pressure_star) const {
280	return pressure_star +
281	(data.pressure_ - pressure_star) *	16✔
282	step_function(direction_ * (x_shifted - shock_speed_ * t));	8✔
283	}
284
285	template <size_t Dim>
286	RiemannProblem<Dim>::Rarefaction::Rarefaction(const InitialData& data,	48✔
287	const double pressure_ratio,
288	const double velocity_star,
289	const double adiabatic_index,
290	const Side& side)
291	: direction_(side == Side::Left ? -1.0 : 1.0) {	48✔
292	gamma_mm_ = adiabatic_index - 1.0;	48✔
293	gamma_pp_ = adiabatic_index + 1.0;	48✔
294	mass_density_star_ =	48✔
295	data.mass_density_ * pow(pressure_ratio, 1.0 / adiabatic_index);	48✔
296	sound_speed_star_ =	48✔
297	data.sound_speed_ *	48✔
298	pow(pressure_ratio, 0.5 * (adiabatic_index - 1.0) / adiabatic_index);	48✔
299	head_speed_ = data.normal_velocity_ + direction_ * data.sound_speed_;	48✔
300	tail_speed_ = velocity_star + direction_ * sound_speed_star_;	48✔
301	}	48✔
302
303	template <size_t Dim>
304	double RiemannProblem<Dim>::Rarefaction::mass_density(	136✔
305	const double x_shifted, const double t, const InitialData& data) const {
306	const double s = (t > 0.0 ? (x_shifted / t) : 0.0);	136✔
307	return direction_ * (x_shifted - tail_speed_ * t) < 0.0	136✔
308	? mass_density_star_	200✔
309	: (direction_ * (x_shifted - tail_speed_ * t) >= 0.0 and	128✔
310	direction_ * (x_shifted - head_speed_ * t) < 0.0	64✔
311	? (data.mass_density_ *	64✔
312	pow((2.0 - direction_ * gamma_mm_ *	64✔
313	(data.normal_velocity_ - s) /	64✔
314	data.sound_speed_) /	64✔
315	gamma_pp_,	64✔
316	2.0 / gamma_mm_))	64✔
317	: data.mass_density_);	136✔
318	}
319
320	template <size_t Dim>
321	double RiemannProblem<Dim>::Rarefaction::normal_velocity(	68✔
322	const double x_shifted, const double t, const InitialData& data,
323	const double velocity_star) const {
324	const double s = (t > 0.0 ? (x_shifted / t) : 0.0);	68✔
325	return direction_ * (x_shifted - tail_speed_ * t) < 0.0	68✔
326	? velocity_star	100✔
327	: (direction_ * (x_shifted - tail_speed_ * t) >= 0.0 and	64✔
328	direction_ * (x_shifted - head_speed_ * t) < 0.0	32✔
329	? (2.0 *	32✔
330	(0.5 * gamma_mm_ * data.normal_velocity_ + s -	32✔
331	direction_ * data.sound_speed_) /	32✔
332	gamma_pp_)	32✔
333	: data.normal_velocity_);	68✔
334	}
335
336	template <size_t Dim>
337	double RiemannProblem<Dim>::Rarefaction::pressure(	136✔
338	const double x_shifted, const double t, const InitialData& data,
339	const double pressure_star) const {
340	const double s = (t > 0.0 ? (x_shifted / t) : 0.0);	136✔
341	return direction_ * (x_shifted - tail_speed_ * t) < 0.0	136✔
342	? pressure_star	200✔
343	: (direction_ * (x_shifted - tail_speed_ * t) >= 0.0 and	128✔
344	direction_ * (x_shifted - head_speed_ * t) < 0.0	64✔
345	? (data.pressure_ *	64✔
346	pow((2.0 - direction_ * gamma_mm_ *	64✔
347	(data.normal_velocity_ - s) /	64✔
348	data.sound_speed_) /	64✔
349	gamma_pp_,	64✔
350	2.0 * (gamma_mm_ + 1.0) / gamma_mm_))	64✔
351	: data.pressure_);	136✔
352	}
353
354	template <size_t Dim>
355	template <typename DataType>
356	tuples::TaggedTuple<Tags::MassDensity<DataType>> RiemannProblem<Dim>::variables(	48✔
357	const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
358	tmpl::list<Tags::MassDensity<DataType>> /meta/, const Wave& left,
359	const Wave& right) const {
360	auto mass_density = make_with_value<Scalar<DataType>>(x_shifted, 0.0);	72✔
361	const double u_star_times_t = velocity_star_ * t;	48✔
362	for (size_t s = 0; s < get_size(get<0>(x_shifted)); ++s) {	432✔
363	const double x_shifted_s = get_element(x_shifted.get(propagation_axis_), s);	144✔
364	get_element(get(mass_density), s) =	432✔
365	(x_shifted_s < u_star_times_t ? left.mass_density(x_shifted_s, t)	144✔
366	: right.mass_density(x_shifted_s, t));
367	}
368	return mass_density;	72✔
369	}
370
371	template <size_t Dim>
372	template <typename DataType>
373	tuples::TaggedTuple<Tags::Velocity<DataType, Dim>>
374	RiemannProblem<Dim>::variables(	24✔
375	const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
376	tmpl::list<Tags::Velocity<DataType, Dim>> /meta/, const Wave& left,
377	const Wave& right) const {
378	auto velocity = make_with_value<tnsr::I<DataType, Dim, Frame::Inertial>>(	20✔
379	get<0>(x_shifted), 0.0);	48✔
380
381	const double u_star_times_t = velocity_star_ * t;	24✔
382	for (size_t s = 0; s < get_size(get<0>(x_shifted)); ++s) {	264✔
383	const double x_shifted_s = get_element(x_shifted.get(propagation_axis_), s);	72✔
384
385	size_t index = propagation_axis_ % Dim;	72✔
386	get_element(velocity.get(index), s) =	96✔
387	(x_shifted_s < u_star_times_t
388	? left.normal_velocity(x_shifted_s, t, velocity_star_)	76✔
389	: right.normal_velocity(x_shifted_s, t, velocity_star_));	4✔
390
391	if (Dim > 1) {
392	index = (propagation_axis_ + 1) % Dim;	48✔
393	get_element(velocity.get(index), s) =	72✔
394	(x_shifted_s < u_star_times_t
395	? gsl::at(left_initial_data_.velocity_, index)	92✔
396	: gsl::at(right_initial_data_.velocity_, index));	8✔
397	}
398
399	if (Dim > 2) {
400	index = (propagation_axis_ + 2) % Dim;	24✔
401	get_element(velocity.get(index), s) =	48✔
402	(x_shifted_s < u_star_times_t
403	? gsl::at(left_initial_data_.velocity_, index)	48✔
404	: gsl::at(right_initial_data_.velocity_, index));	×
405	}
406	}
407	return velocity;	40✔
408	}
409
410	template <size_t Dim>
411	template <typename DataType>
412	tuples::TaggedTuple<Tags::Pressure<DataType>> RiemannProblem<Dim>::variables(	48✔
413	const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
414	tmpl::list<Tags::Pressure<DataType>> /meta/, const Wave& left,
415	const Wave& right) const {
416	auto pressure = make_with_value<Scalar<DataType>>(x_shifted, 0.0);	72✔
417	const double u_star_times_t = velocity_star_ * t;	48✔
418	for (size_t s = 0; s < get_size(get<0>(x_shifted)); ++s) {	432✔
419	const double x_shifted_s = get_element(x_shifted.get(propagation_axis_), s);	144✔
420	get_element(get(pressure), s) =	432✔
421	(x_shifted_s < u_star_times_t
422	? left.pressure(x_shifted_s, t, pressure_star_)	152✔
423	: right.pressure(x_shifted_s, t, pressure_star_));	8✔
424	}
425	return pressure;	72✔
426	}
427
428	template <size_t Dim>
429	template <typename DataType>
430	tuples::TaggedTuple<Tags::SpecificInternalEnergy<DataType>>
431	RiemannProblem<Dim>::variables(	24✔
432	const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, const double t,
433	tmpl::list<Tags::SpecificInternalEnergy<DataType>> /meta/,
434	const Wave& left, const Wave& right) const {
435	return equation_of_state_.specific_internal_energy_from_density_and_pressure(
436	get<Tags::MassDensity<DataType>>(	24✔
437	variables(x_shifted, t, tmpl::list<Tags::MassDensity<DataType>>{},
438	left, right)),
439	get<Tags::Pressure<DataType>>(variables(	36✔
440	x_shifted, t, tmpl::list<Tags::Pressure<DataType>>{}, left, right)));	60✔
441	}
442
443	template <size_t Dim>
444	bool operator==(const RiemannProblem<Dim>& lhs,	48✔
445	const RiemannProblem<Dim>& rhs) {
446	return lhs.adiabatic_index_ == rhs.adiabatic_index_ and	96✔
447	lhs.initial_position_ == rhs.initial_position_ and	48✔
448	lhs.left_initial_data_ == rhs.left_initial_data_ and	48✔
449	lhs.right_initial_data_ == rhs.right_initial_data_ and	48✔
450	lhs.pressure_star_tol_ == rhs.pressure_star_tol_ and	48✔
451	lhs.pressure_star_ == rhs.pressure_star_ and	144✔
452	lhs.velocity_star_ == rhs.velocity_star_;	96✔
453	}
454
455	template <size_t Dim>
456	bool operator!=(const RiemannProblem<Dim>& lhs,	6✔
457	const RiemannProblem<Dim>& rhs) {
458	return not(lhs == rhs);	6✔
459	}
460
461	#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
462	#define DTYPE(data) BOOST_PP_TUPLE_ELEM(1, data)
463	#define TAG(data) BOOST_PP_TUPLE_ELEM(2, data)
464
465	#define INSTANTIATE_CLASS(_, data) \
466	template class RiemannProblem<DIM(data)>; \
467	template bool operator==(const RiemannProblem<DIM(data)>&, \
468	const RiemannProblem<DIM(data)>&); \
469	template bool operator!=(const RiemannProblem<DIM(data)>&, \
470	const RiemannProblem<DIM(data)>&);
471
472	GENERATE_INSTANTIATIONS(INSTANTIATE_CLASS, (1, 2, 3))
473
474	#define INSTANTIATE_SCALARS(_, data) \
475	template tuples::TaggedTuple<TAG(data) < DTYPE(data)> > \
476	RiemannProblem<DIM(data)>::variables( \
477	const tnsr::I<DTYPE(data), DIM(data), Frame::Inertial>& x_shifted, \
478	const double t, tmpl::list<TAG(data) < DTYPE(data)> >, \
479	const Wave& left, const Wave& right) const;
480
481	GENERATE_INSTANTIATIONS(INSTANTIATE_SCALARS, (1, 2, 3), (double, DataVector),
482	(Tags::MassDensity, Tags::Pressure,
483	Tags::SpecificInternalEnergy))
484
485	#define INSTANTIATE_VELOCITY(_, data) \
486	template tuples::TaggedTuple<TAG(data) < DTYPE(data), DIM(data)> > \
487	RiemannProblem<DIM(data)>::variables( \
488	const tnsr::I<DTYPE(data), DIM(data), Frame::Inertial>& x_shifted, \
489	const double t, tmpl::list<TAG(data) < DTYPE(data), DIM(data)> >, \
490	const Wave& left, const Wave& right) const;
491
492	GENERATE_INSTANTIATIONS(INSTANTIATE_VELOCITY, (1, 2, 3), (double, DataVector),
493	(Tags::Velocity))
494
495	#undef DIM
496	#undef DTYPE
497	#undef TAG
498	#undef INSTANTIATE_CLASS
499	#undef INSTANTIATE_SCALARS
500	#undef INSTANTIATE_VELOCITY
501
502	} // namespace NewtonianEuler::Solutions

sxs-collaboration / spectre / 4237779830

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