10#include <deal.II/base/tensor.h>
18 namespace EulerInitialStates
32 template <
typename Description,
int dim,
typename Number>
47 const std::string subsection)
50 , hyperbolic_system_(hyperbolic_system)
53 if constexpr (!View::have_gamma) {
54 this->add_parameter(
"gamma", gamma_,
"The ratio of specific heats");
57 primitive_left_[0] = 1.4;
58 primitive_left_[1] = 0.0;
59 primitive_left_[2] = 1.0;
61 "primitive state left",
63 "Initial 1d primitive state (rho, u, p) on the left");
65 primitive_right_[0] = 1.4;
66 primitive_right_[1] = 0.0;
67 primitive_right_[2] = 1.0;
69 "primitive state right",
71 "Initial 1d primitive state (rho, u, p) on the right");
74 const auto prepare_riemann_data = [&]() {
75 const auto view = hyperbolic_system_.template view<dim, Number>();
76 if constexpr (View::have_gamma) {
77 gamma_ = view.gamma();
80 const Number p_L = primitive_left_[2];
81 const Number p_R = primitive_right_[2];
83 p_star_ = compute_pstar(p_L, p_R, primitive_left_, primitive_right_);
85 const Number u_L = primitive_left_[1];
86 u_star_ = u_L - fZofP(p_star_, primitive_left_);
89 const Number u_R = primitive_right_[1];
90 std::cout <<
"left data = " << primitive_left_
91 <<
"\nright data = " << primitive_right_
92 <<
"\np_star = " << p_star_
93 <<
"\nu_star = " << u_star_
94 <<
"\nVerifying u_star = "
95 << u_R + fZofP(p_star_, primitive_right_) << std::endl;
98 lambda_left_minus_ = lambda(p_star_, primitive_left_, -1.);
100 lambda_intermediate(p_star_, primitive_left_, -1.);
101 lambda_right_minus_ =
102 lambda_intermediate(p_star_, primitive_right_, 1.);
103 lambda_right_plus_ = lambda(p_star_, primitive_right_, 1.);
107 std::cout <<
"lambda_left_minus = " << lambda_left_minus_
108 <<
"\nlambda_left_plus = " << lambda_left_plus_
109 <<
"\nlambda_right_minus = " << lambda_right_minus_
110 <<
"\nlambda_right_plus = " << lambda_right_plus_
115 this->parse_parameters_call_back.connect(prepare_riemann_data);
116 prepare_riemann_data();
122 const auto view = hyperbolic_system_.template view<dim, Number>();
124 const double &x = point[0];
126 const Number xi = x / t;
128 dealii::Tensor<1, 3, Number> primitive_state;
130 if (t < 1.e-14 && x < 0.) {
131 primitive_state = primitive_left_;
133 std::cout <<
"Left primitive state: " << primitive_state << std::endl;
136 }
else if (t < 1.e-14 && x > 0.) {
137 primitive_state = primitive_right_;
139 std::cout <<
"Right primitive state: " << primitive_state
143 }
else if (xi < lambda_left_minus_) {
145 primitive_state = primitive_left_;
147 std::cout <<
"Left primitive state: " << primitive_state << std::endl;
150 }
else if (xi < lambda_left_plus_) {
152 expansion_solution(p_star_, xi, primitive_left_, -1.);
153 primitive_state = c_LL;
155 std::cout <<
"Left expansion state: " << primitive_state << std::endl;
158 }
else if (xi < u_star_) {
159 primitive_state = cstar_solution(p_star_, u_star_, primitive_left_);
161 const Number p_L = primitive_left_[2];
163 primitive_state = expansion_solution(
164 p_star_, lambda_left_plus_, primitive_left_, -1.);
166 std::cout <<
"Left cstar state: " << primitive_state << std::endl;
169 }
else if (xi < lambda_right_minus_) {
170 primitive_state = cstar_solution(p_star_, u_star_, primitive_right_);
172 const Number p_R = primitive_right_[2];
174 primitive_state = expansion_solution(
175 p_star_, lambda_right_minus_, primitive_right_, 1.);
177 std::cout <<
"Right cstar state: " << primitive_state << std::endl;
180 }
else if (xi < lambda_right_plus_) {
182 expansion_solution(p_star_, xi, primitive_right_, 1.);
184 std::cout <<
"Right expansion state: " << primitive_state
190 primitive_state = primitive_right_;
192 std::cout <<
"Right primitive state: " << primitive_state
197 return view.from_initial_state(primitive_state);
209 dealii::Tensor<1, 3, Number> primitive_left_;
210 dealii::Tensor<1, 3, Number> primitive_right_;
222 Number lambda_left_minus_;
223 Number lambda_left_plus_;
224 Number lambda_right_minus_;
225 Number lambda_right_plus_;
233 Number fZofP(
const Number &p_in,
234 const dealii::Tensor<1, 3, Number> &data_in)
const
237 const Number rho_Z = data_in[0];
238 const Number p_Z = data_in[2];
240 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
242 const Number A_Z = 2. / (gamma_ + 1.) / rho_Z;
243 const Number B_Z = (gamma_ - 1.) / (gamma_ + 1.) * p_Z;
245 const Number exp = 0.5 * (gamma_ - 1.) / gamma_;
246 Number left_brach = 2. * c_Z / (gamma_ - 1.);
247 left_brach *= (
std::pow(p_in / p_Z, exp) - 1.);
249 Number f_of_p = (p_in - p_Z) * std::sqrt(A_Z / (p_in + B_Z));
258 Number dfZofP(
const Number &p_in,
259 const dealii::Tensor<1, 3, Number> &data_in)
const
262 const Number rho_Z = data_in[0];
263 const Number p_Z = data_in[2];
265 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
267 const Number A_Z = 2. / (gamma_ + 1.) / rho_Z;
268 const Number B_Z = (gamma_ - 1.) / (gamma_ + 1.) * p_Z;
270 Number exp = 0.5 * (gamma_ - 1.) / gamma_;
271 Number left_brach = 2. * c_Z / (gamma_ - 1.) * exp;
274 left_brach *=
std::pow(p_in / p_Z, exp - 1.) / p_Z;
276 Number right_branch =
std::pow(A_Z / (p_in + B_Z), 1.5);
277 right_branch *= (2. * B_Z + p_in + p_Z) / (2. * A_Z);
279 Number df_of_p = right_branch;
282 df_of_p = left_brach;
288 Number dphi(
const Number &p_in,
289 const dealii::Tensor<1, 3, Number> &data_left,
290 const dealii::Tensor<1, 3, Number> &data_right)
const
292 return dfZofP(p_in, data_left) + dfZofP(p_in, data_right);
296 Number phi(
const Number &p_in,
297 const dealii::Tensor<1, 3, Number> &data_left,
298 const dealii::Tensor<1, 3, Number> &data_right)
const
300 const Number u_L = data_left[1];
301 const Number u_R = data_right[1];
303 return fZofP(p_in, data_right) + fZofP(p_in, data_left) + u_R - u_L;
307 Number lambda(
const Number &p_in,
308 const dealii::Tensor<1, 3, Number> &data_in,
309 const Number &sign)
const
312 const Number rho_Z = data_in[0];
313 const Number u_Z = data_in[1];
314 const Number p_Z = data_in[2];
316 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
318 const Number radicand =
319 1. + 0.5 * (gamma_ + 1.) / gamma_ * std::max(p_in / p_Z - 1., 0.);
321 return u_Z + sign * c_Z * std::sqrt(radicand);
325 Number lambda_intermediate(
const Number &p_in,
326 const dealii::Tensor<1, 3, Number> &data_in,
327 const Number &sign)
const
329 const Number rho_Z = data_in[0];
330 const Number u_Z = data_in[1];
331 const Number p_Z = data_in[2];
333 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
335 const auto lambda_value = lambda(p_in, data_in, sign);
337 const Number f_of_p = fZofP(p_in, data_in);
339 const Number exp = 0.5 * (gamma_ - 1.) / gamma_;
340 const Number expansion_speed =
341 u_Z + sign * (f_of_p + c_Z *
std::pow(p_in / p_Z, exp));
343 Number result = lambda_value;
345 result = expansion_speed;
351 dealii::Tensor<1, 3, Number>
352 cstar_solution(
const Number &p_star,
353 const Number &u_star,
354 const dealii::Tensor<1, 3, Number> &data_in)
const
356 const Number rho_Z = data_in[0];
357 const Number p_Z = data_in[2];
360 const Number p_ratio = p_star / p_Z;
361 const Number gamma_ratio = (gamma_ - 1.) / (gamma_ + 1.);
363 const Number numerator = rho_Z * (p_ratio + gamma_ratio);
364 const Number denominator = gamma_ratio * p_ratio + 1.;
366 Number rho_star = numerator / denominator;
368 auto result = data_in;
369 result[0] = rho_star;
377 dealii::Tensor<1, 3, Number>
378 expansion_solution(
const Number & ,
380 const dealii::Tensor<1, 3, Number> &data_in,
381 const Number &sign)
const
383 const Number rho_Z = data_in[0];
384 const Number u_Z = data_in[1];
385 const Number p_Z = data_in[2];
387 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
390 const Number gamma_ratio = (gamma_ - 1.) / (gamma_ + 1.);
392 const Number first = 2. / (gamma_ + 1.);
393 const Number second = gamma_ratio / c_Z * (u_Z - xi);
394 const Number exp = 2. / (gamma_ - 1.);
396 Number rho_expansion = rho_Z *
std::pow(first - sign * second, exp);
399 const Number factor = p_Z /
std::pow(rho_Z, gamma_);
400 const Number p_expansion = factor *
std::pow(rho_expansion, gamma_);
403 const Number u_expansion = u_Z + sign * fZofP(p_expansion, data_in);
405 auto result = data_in;
406 result[0] = rho_expansion;
407 result[1] = u_expansion;
408 result[2] = p_expansion;
417 double compute_pstar(
double p_1,
419 dealii::Tensor<1, 3, Number> data_1,
420 dealii::Tensor<1, 3, Number> data_2)
422 constexpr Number eps = std::numeric_limits<Number>::epsilon();
428 std::swap(data_1, data_2);
432 const double phi_1 = phi(p_1, data_1, data_2);
433 const double phi_2 = phi(p_2, data_1, data_2);
434 Assert(phi_1 * phi_2 <= 0.,
436 "Euler::ExactRiemannSolver: failed to compute p_star."));
446 std::cout <<
"Computing p_star with a bisection method." << std::endl;
449 unsigned int iter = 0;
450 for (; iter < 200; ++iter) {
453 if (std::abs(p_2 - p_1) < 10. * eps * std::max(p_1, p_2)) {
457 const double phi_2 = phi(p_2, data_1, data_2);
460 const double phi_1 = phi(p_1, data_1, data_2);
462 std::cout <<
"\niter: " << iter <<
"\n";
463 std::cout <<
"p_1: " << p_1 <<
"\n";
464 std::cout <<
"p_2: " << p_2 <<
"\n";
465 std::cout <<
"phi_1: " << phi_1 <<
"\n";
466 std::cout <<
"phi_2: " << phi_2 <<
"\n";
469 const auto p_m = 0.5 * (p_2 + p_1);
470 const double phi_m = phi(p_m, data_1, data_2);
472 if (phi_m * phi_2 >= 0.) {
480 const double phi_2 = phi(p_2, data_1, data_2);
481 std::cout <<
"After " << iter <<
" iterations:"
482 <<
"\np_star = " << p_2 <<
"\nphi(p_star) = " << phi_2
483 <<
"\n|p_2 - p_1| = " << std::abs(p_2 - p_1) << std::endl;
typename View::state_type state_type
typename Description::template HyperbolicSystemView< dim, Number > View
typename Description::HyperbolicSystem HyperbolicSystem
ExactRiemannSolution(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
typename View::ScalarNumber ScalarNumber
state_type compute(const dealii::Point< dim > &point, Number t) final
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem