ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
initial_state_exact_riemann_solution.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2023 - 2024 by the ryujin authors
4// Copyright (C) 2025 by Triad National Security, LLC
5//
6
7#pragma once
8
9#include <cmath>
10#include <deal.II/base/tensor.h>
11
13
14// #define DEBUG_SOLUTION
15
16namespace ryujin
17{
18 namespace EulerInitialStates
19 {
32 template <typename Description, int dim, typename Number>
33 class ExactRiemannSolution : public InitialState<Description, dim, Number>
34 {
35 public:
37
39 using View =
40 typename Description::template HyperbolicSystemView<dim, Number>;
41 using state_type = typename View::state_type;
42
43 using ScalarNumber = typename View::ScalarNumber;
44
45
46 ExactRiemannSolution(const HyperbolicSystem &hyperbolic_system,
47 const std::string subsection)
48 : InitialState<Description, dim, Number>("exact riemann solution",
49 subsection)
50 , hyperbolic_system_(hyperbolic_system)
51 {
52 gamma_ = 1.4;
53 if constexpr (!View::have_gamma) {
54 this->add_parameter("gamma", gamma_, "The ratio of specific heats");
55 }
56
57 primitive_left_[0] = 1.4;
58 primitive_left_[1] = 0.0;
59 primitive_left_[2] = 1.0;
60 this->add_parameter(
61 "primitive state left",
62 primitive_left_,
63 "Initial 1d primitive state (rho, u, p) on the left");
64
65 primitive_right_[0] = 1.4;
66 primitive_right_[1] = 0.0;
67 primitive_right_[2] = 1.0;
68 this->add_parameter(
69 "primitive state right",
70 primitive_right_,
71 "Initial 1d primitive state (rho, u, p) on the right");
72
73 // Convert the primitive states to conserved states
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();
78 }
79
80 const Number p_L = primitive_left_[2];
81 const Number p_R = primitive_right_[2];
82
83 p_star_ = compute_pstar(p_L, p_R, primitive_left_, primitive_right_);
84
85 const Number u_L = primitive_left_[1];
86 u_star_ = u_L - fZofP(p_star_, primitive_left_);
87
88#ifdef DEBUG_SOLUTION
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;
96#endif
97
98 lambda_left_minus_ = lambda(p_star_, primitive_left_, -1.);
99 lambda_left_plus_ =
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.);
104
105
106#ifdef DEBUG_SOLUTION
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_
111 << std::endl;
112#endif
113 };
114
115 this->parse_parameters_call_back.connect(prepare_riemann_data);
116 prepare_riemann_data();
117 }
118
119
120 state_type compute(const dealii::Point<dim> &point, Number t) final
121 {
122 const auto view = hyperbolic_system_.template view<dim, Number>();
123
124 const double &x = point[0];
125
126 const Number xi = x / t;
127
128 dealii::Tensor<1, 3, Number> primitive_state;
129
130 if (t < 1.e-14 && x < 0.) {
131 primitive_state = primitive_left_;
132#ifdef DEBUG_SOLUTION
133 std::cout << "Left primitive state: " << primitive_state << std::endl;
134#endif
135
136 } else if (t < 1.e-14 && x > 0.) {
137 primitive_state = primitive_right_;
138#ifdef DEBUG_SOLUTION
139 std::cout << "Right primitive state: " << primitive_state
140 << std::endl;
141#endif
142
143 } else if (xi < lambda_left_minus_) {
144 /* Left state: */
145 primitive_state = primitive_left_;
146#ifdef DEBUG_SOLUTION
147 std::cout << "Left primitive state: " << primitive_state << std::endl;
148#endif
149
150 } else if (xi < lambda_left_plus_) {
151 const auto c_LL =
152 expansion_solution(p_star_, xi, primitive_left_, -1.);
153 primitive_state = c_LL;
154#ifdef DEBUG_SOLUTION
155 std::cout << "Left expansion state: " << primitive_state << std::endl;
156#endif
157
158 } else if (xi < u_star_) {
159 primitive_state = cstar_solution(p_star_, u_star_, primitive_left_);
160
161 const Number p_L = primitive_left_[2];
162 if (p_star_ < p_L)
163 primitive_state = expansion_solution(
164 p_star_, lambda_left_plus_, primitive_left_, -1.);
165#ifdef DEBUG_SOLUTION
166 std::cout << "Left cstar state: " << primitive_state << std::endl;
167#endif
168
169 } else if (xi < lambda_right_minus_) {
170 primitive_state = cstar_solution(p_star_, u_star_, primitive_right_);
171
172 const Number p_R = primitive_right_[2];
173 if (p_star_ < p_R)
174 primitive_state = expansion_solution(
175 p_star_, lambda_right_minus_, primitive_right_, 1.);
176#ifdef DEBUG_SOLUTION
177 std::cout << "Right cstar state: " << primitive_state << std::endl;
178#endif
179
180 } else if (xi < lambda_right_plus_) {
181 primitive_state =
182 expansion_solution(p_star_, xi, primitive_right_, 1.);
183#ifdef DEBUG_SOLUTION
184 std::cout << "Right expansion state: " << primitive_state
185 << std::endl;
186#endif
187
188 } else {
189 /* Right state: */
190 primitive_state = primitive_right_;
191#ifdef DEBUG_SOLUTION
192 std::cout << "Right primitive state: " << primitive_state
193 << std::endl;
194#endif
195 }
196
197 return view.from_initial_state(primitive_state);
198 }
199
200 private:
202
206
207 Number gamma_;
208
209 dealii::Tensor<1, 3, Number> primitive_left_;
210 dealii::Tensor<1, 3, Number> primitive_right_;
211
213
217
218 const HyperbolicSystem &hyperbolic_system_;
219
220 Number p_star_;
221 Number u_star_;
222 Number lambda_left_minus_;
223 Number lambda_left_plus_;
224 Number lambda_right_minus_;
225 Number lambda_right_plus_;
226
228
232
233 Number fZofP(const Number &p_in,
234 const dealii::Tensor<1, 3, Number> &data_in) const
235 {
236 // Get left/right data
237 const Number rho_Z = data_in[0];
238 const Number p_Z = data_in[2];
239
240 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
241
242 const Number A_Z = 2. / (gamma_ + 1.) / rho_Z;
243 const Number B_Z = (gamma_ - 1.) / (gamma_ + 1.) * p_Z;
244
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.);
248
249 Number f_of_p = (p_in - p_Z) * std::sqrt(A_Z / (p_in + B_Z));
250
251 if (p_in <= p_Z)
252 f_of_p = left_brach;
253
254 return f_of_p;
255 }
256
257
258 Number dfZofP(const Number &p_in,
259 const dealii::Tensor<1, 3, Number> &data_in) const
260 {
261 // Get left/right data
262 const Number rho_Z = data_in[0];
263 const Number p_Z = data_in[2];
264
265 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
266
267 const Number A_Z = 2. / (gamma_ + 1.) / rho_Z;
268 const Number B_Z = (gamma_ - 1.) / (gamma_ + 1.) * p_Z;
269
270 Number exp = 0.5 * (gamma_ - 1.) / gamma_;
271 Number left_brach = 2. * c_Z / (gamma_ - 1.) * exp;
272 exp -= 1.;
273
274 left_brach *= std::pow(p_in / p_Z, exp - 1.) / p_Z;
275
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);
278
279 Number df_of_p = right_branch;
280
281 if (p_in <= p_Z)
282 df_of_p = left_brach;
283
284 return df_of_p;
285 }
286
287
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
291 {
292 return dfZofP(p_in, data_left) + dfZofP(p_in, data_right);
293 }
294
295
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
299 {
300 const Number u_L = data_left[1];
301 const Number u_R = data_right[1];
302
303 return fZofP(p_in, data_right) + fZofP(p_in, data_left) + u_R - u_L;
304 }
305
306
307 Number lambda(const Number &p_in,
308 const dealii::Tensor<1, 3, Number> &data_in,
309 const Number &sign) const
310 {
311 // Get left/right data
312 const Number rho_Z = data_in[0];
313 const Number u_Z = data_in[1];
314 const Number p_Z = data_in[2];
315
316 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
317
318 const Number radicand =
319 1. + 0.5 * (gamma_ + 1.) / gamma_ * std::max(p_in / p_Z - 1., 0.);
320
321 return u_Z + sign * c_Z * std::sqrt(radicand);
322 }
323
324
325 Number lambda_intermediate(const Number &p_in,
326 const dealii::Tensor<1, 3, Number> &data_in,
327 const Number &sign) const
328 {
329 const Number rho_Z = data_in[0];
330 const Number u_Z = data_in[1];
331 const Number p_Z = data_in[2];
332
333 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
334
335 const auto lambda_value = lambda(p_in, data_in, sign);
336
337 const Number f_of_p = fZofP(p_in, data_in);
338
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));
342
343 Number result = lambda_value;
344 if (p_in < p_Z)
345 result = expansion_speed;
346
347 return result;
348 }
349
350
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
355 {
356 const Number rho_Z = data_in[0];
357 const Number p_Z = data_in[2];
358
359 // Define rho_star
360 const Number p_ratio = p_star / p_Z;
361 const Number gamma_ratio = (gamma_ - 1.) / (gamma_ + 1.);
362
363 const Number numerator = rho_Z * (p_ratio + gamma_ratio);
364 const Number denominator = gamma_ratio * p_ratio + 1.;
365
366 Number rho_star = numerator / denominator;
367
368 auto result = data_in;
369 result[0] = rho_star;
370 result[1] = u_star;
371 result[2] = p_star;
372
373 return result;
374 }
375
376
377 dealii::Tensor<1, 3, Number>
378 expansion_solution(const Number & /*p_star*/,
379 const Number &xi,
380 const dealii::Tensor<1, 3, Number> &data_in,
381 const Number &sign) const
382 {
383 const Number rho_Z = data_in[0];
384 const Number u_Z = data_in[1];
385 const Number p_Z = data_in[2];
386
387 const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z);
388
389 // Define rho_expansion
390 const Number gamma_ratio = (gamma_ - 1.) / (gamma_ + 1.);
391
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.);
395
396 Number rho_expansion = rho_Z * std::pow(first - sign * second, exp);
397
398 // Define p_expansion
399 const Number factor = p_Z / std::pow(rho_Z, gamma_);
400 const Number p_expansion = factor * std::pow(rho_expansion, gamma_);
401
402 // Define u_expansion
403 const Number u_expansion = u_Z + sign * fZofP(p_expansion, data_in);
404
405 auto result = data_in;
406 result[0] = rho_expansion;
407 result[1] = u_expansion;
408 result[2] = p_expansion;
409
410 return result;
411 }
412
413
417 double compute_pstar(double p_1,
418 double p_2,
419 dealii::Tensor<1, 3, Number> data_1,
420 dealii::Tensor<1, 3, Number> data_2)
421 {
422 constexpr Number eps = std::numeric_limits<Number>::epsilon();
423
424 // Ensure that p_1 <= p_2
425
426 if (p_1 > p_2) {
427 std::swap(p_1, p_2);
428 std::swap(data_1, data_2);
429 }
430
431#ifdef DEBUG
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.,
435 dealii::ExcMessage(
436 "Euler::ExactRiemannSolver: failed to compute p_star."));
437#endif
438
439 //
440 // We simply compute the root of phi with a bisection method down
441 // to machine precision. This is not terribly efficient but luckily
442 // happens only once during initialization.
443 //
444
445#ifdef DEBUG_SOLUTION
446 std::cout << "Computing p_star with a bisection method." << std::endl;
447#endif
448
449 unsigned int iter = 0;
450 for (; iter < 200; ++iter) {
451
452 // Check for convergence:
453 if (std::abs(p_2 - p_1) < 10. * eps * std::max(p_1, p_2)) {
454 break;
455 }
456
457 const double phi_2 = phi(p_2, data_1, data_2);
458
459#ifdef DEBUG_SOLUTION
460 const double phi_1 = phi(p_1, data_1, data_2);
461
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";
467#endif
468
469 const auto p_m = 0.5 * (p_2 + p_1);
470 const double phi_m = phi(p_m, data_1, data_2);
471
472 if (phi_m * phi_2 >= 0.) {
473 p_2 = p_m;
474 } else {
475 p_1 = p_m;
476 }
477 }
478
479#ifdef DEBUG_SOLUTION
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;
484#endif
485
486 return p_2;
487 }
488
490 };
491 } // namespace EulerInitialStates
492} // namespace ryujin
typename Description::template HyperbolicSystemView< dim, Number > View
ExactRiemannSolution(const HyperbolicSystem &hyperbolic_system, const std::string subsection)
state_type compute(const dealii::Point< dim > &point, Number t) final
T pow(const T x, const T b)
Euler::HyperbolicSystem HyperbolicSystem
Definition: description.h:32