Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
riemann_common.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
21
22#include "shambackends/math.hpp"
24#include "shambackends/vec.hpp"
25#include <algorithm>
26#include <array>
27#include <cmath>
28#include <concepts>
29#include <iostream>
30#include <utility>
31namespace shammath {
32
38 template<class T>
39 concept FluidStateSpec = requires(
40 const T &self,
41 typename T::Tcons cons,
42 typename T::Tprim prim,
43 typename T::Tvec n,
44 typename T::Tscal vn) {
45 typename T::Tvec;
46 typename T::Tscal;
47 typename T::Tprim;
48 typename T::Tcons;
49 { self.cons_to_prim(cons) } -> std::convertible_to<typename T::Tprim>;
50 { self.prim_to_cons(prim) } -> std::convertible_to<typename T::Tcons>;
51 { self.sound_speed(prim) } -> std::convertible_to<typename T::Tscal>;
52 { self.vn(prim, n) } -> std::convertible_to<typename T::Tscal>;
53 { self.flux(prim, n) } -> std::convertible_to<typename T::Tcons>;
54 { self.flux(prim, n, vn) } -> std::convertible_to<typename T::Tcons>;
55 };
56
63 template<class T>
64 concept DustFluidStateSpec = requires(
65 const T &self,
66 typename T::Tcons cons,
67 typename T::Tprim prim,
68 typename T::Tvec n,
69 typename T::Tscal vn) {
70 typename T::Tvec;
71 typename T::Tscal;
72 typename T::Tprim;
73 typename T::Tcons;
74 { self.cons_to_prim(cons) } -> std::convertible_to<typename T::Tprim>;
75 { self.prim_to_cons(prim) } -> std::convertible_to<typename T::Tcons>;
76 { self.vn(prim, n) } -> std::convertible_to<typename T::Tscal>;
77 { self.flux(prim, n) } -> std::convertible_to<typename T::Tcons>;
78 { self.flux(prim, n, vn) } -> std::convertible_to<typename T::Tcons>;
79 };
80
81 namespace details {
83 template<class T>
84 concept HasGlobalGamma = requires(const T &self) {
85 { self.gamma() } -> std::convertible_to<typename T::Tscal>;
86 };
87
89 template<class T>
90 concept HasPerStateGamma = requires(const T &self, typename T::Tprim prim) {
91 { self.gamma(prim) } -> std::convertible_to<typename T::Tscal>;
92 };
93 } // namespace details
94
102 template<class T>
105
113 template<FluidStateAdiabaticSpec FSpec>
114 inline constexpr std::pair<typename FSpec::Tscal, typename FSpec::Tscal> get_adiabatic_index_lr(
115 const FSpec &fspec,
116 const typename FSpec::Tprim &prim_l,
117 const typename FSpec::Tprim &prim_r) {
118 if constexpr (details::HasGlobalGamma<FSpec>) {
119 const typename FSpec::Tscal gamma = fspec.gamma();
120 return {gamma, gamma};
121 } else {
122 return {fspec.gamma(prim_l), fspec.gamma(prim_r)};
123 }
124 }
125
126 template<class VecType>
127 struct ConsState {
128 using Tvec = VecType;
129 using Tscal = shambase::VecComponent<Tvec>;
130
131 Tscal rho{}, rhoe{};
132 Tvec rhovel{};
133
134 const ConsState &operator+=(const ConsState &);
135 const ConsState &operator-=(const ConsState &);
136 const ConsState &operator*=(const Tscal);
137 };
138
139 template<class VecType>
140 struct PrimState {
141 using Tvec = VecType;
142 using Tscal = shambase::VecComponent<Tvec>;
143
144 Tscal rho{}, press{};
145 Tvec vel{};
146 };
147
148 template<class Tvec>
149 const ConsState<Tvec> &ConsState<Tvec>::operator+=(const ConsState<Tvec> &cst) {
150 rho += cst.rho;
151 rhoe += cst.rhoe;
152 rhovel += cst.rhovel;
153 return *this;
154 }
155
156 template<class Tvec>
157 const ConsState<Tvec> operator+(const ConsState<Tvec> &lhs, const ConsState<Tvec> &rhs) {
158 return ConsState<Tvec>(lhs) += rhs;
159 }
160
161 template<class Tvec>
162 const ConsState<Tvec> &ConsState<Tvec>::operator-=(const ConsState<Tvec> &cst) {
163 rho -= cst.rho;
164 rhoe -= cst.rhoe;
165 rhovel -= cst.rhovel;
166 return *this;
167 }
168
169 template<class Tvec>
170 const ConsState<Tvec> operator-(const ConsState<Tvec> &lhs, const ConsState<Tvec> &rhs) {
171 return ConsState<Tvec>(lhs) -= rhs;
172 }
173
174 template<class Tvec>
175 const ConsState<Tvec> &ConsState<Tvec>::operator*=(
176 const typename ConsState<Tvec>::Tscal factor) {
177 rho *= factor;
178 rhoe *= factor;
179 rhovel *= factor;
180 return *this;
181 }
182
183 template<class Tvec>
184 const ConsState<Tvec> operator*(
185 const typename ConsState<Tvec>::Tscal factor, const ConsState<Tvec> &rhs) {
186 return ConsState<Tvec>(rhs) *= factor;
187 }
188
189 template<class Tvec>
190 const ConsState<Tvec> operator*(
191 const ConsState<Tvec> &lhs, const typename ConsState<Tvec>::Tscal factor) {
192 return ConsState<Tvec>(lhs) *= factor;
193 }
194
195 template<class VecType>
196 struct Fluxes {
197 using Tvec = VecType;
198 using Tscal = shambase::VecComponent<Tvec>;
199
200 std::array<ConsState<Tvec>, 3> f;
201 };
202
203 template<class Tvec>
204 inline constexpr shambase::VecComponent<Tvec> rhoekin(
205 shambase::VecComponent<Tvec> rho, Tvec v) {
206 using Tscal = shambase::VecComponent<Tvec>;
207 const Tscal v2 = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
208 return 0.5 * rho * v2;
209 }
210
211 template<class Tvec>
212 inline constexpr ConsState<Tvec> prim_to_cons(
213 const PrimState<Tvec> prim, typename PrimState<Tvec>::Tscal gamma) {
214 ConsState<Tvec> cons;
215
216 cons.rho = prim.rho;
217
218 const auto rhoeint = prim.press / (gamma - 1.0);
219 cons.rhoe = rhoeint + rhoekin(prim.rho, prim.vel);
220
221 cons.rhovel[0] = prim.rho * prim.vel[0];
222 cons.rhovel[1] = prim.rho * prim.vel[1];
223 cons.rhovel[2] = prim.rho * prim.vel[2];
224
225 return cons;
226 }
227
228 template<class Tvec>
229 inline constexpr PrimState<Tvec> cons_to_prim(
230 const ConsState<Tvec> cons, typename ConsState<Tvec>::Tscal gamma) {
231 PrimState<Tvec> prim;
232
233 prim.rho = cons.rho;
234
235 prim.vel[0] = cons.rhovel[0] / cons.rho;
236 prim.vel[1] = cons.rhovel[1] / cons.rho;
237 prim.vel[2] = cons.rhovel[2] / cons.rho;
238
239 const auto rhoeint = cons.rhoe - rhoekin(prim.rho, prim.vel);
240 prim.press = (gamma - 1.0) * rhoeint;
241
242 return prim;
243 }
244
251 template<class Tvec>
253 const PrimState<Tvec> prim,
254 Tvec n,
255 typename PrimState<Tvec>::Tscal vn,
256 typename PrimState<Tvec>::Tscal gamma) {
257 ConsState<Tvec> flux;
258
259 const auto rhoeint = prim.press / (gamma - 1.0);
260 const auto rhoe = rhoeint + rhoekin(prim.rho, prim.vel);
261
262 flux.rho = prim.rho * vn;
263
264 flux.rhoe = (rhoe + prim.press) * vn;
265
266 flux.rhovel[0] = prim.rho * vn * prim.vel[0] + prim.press * n[0];
267 flux.rhovel[1] = prim.rho * vn * prim.vel[1] + prim.press * n[1];
268 flux.rhovel[2] = prim.rho * vn * prim.vel[2] + prim.press * n[2];
269
270 return flux;
271 }
272
278 template<class Tvec>
280 const PrimState<Tvec> prim, Tvec n, typename PrimState<Tvec>::Tscal gamma) {
281 const auto vn = n[0] * prim.vel[0] + n[1] * prim.vel[1] + n[2] * prim.vel[2];
282 return hydro_flux_n(prim, n, vn, gamma);
283 }
284
285 template<class Tvec>
286 inline constexpr shambase::VecComponent<Tvec> sound_speed(
287 PrimState<Tvec> prim, shambase::VecComponent<Tvec> gamma) {
288 return sycl::sqrt(gamma * prim.press / prim.rho);
289 }
290
291 template<class Tcons>
292 inline constexpr Tcons y_to_x(const Tcons c) {
293 Tcons cprime;
294 cprime.rho = c.rho;
295 cprime.rhoe = c.rhoe;
296 cprime.rhovel[0] = c.rhovel[1];
297 cprime.rhovel[1] = -c.rhovel[0];
298 cprime.rhovel[2] = c.rhovel[2];
299 return cprime;
300 }
301
302 template<class Tcons>
303 inline constexpr Tcons x_to_y(const Tcons c) {
304 Tcons cprime;
305 cprime.rho = c.rho;
306 cprime.rhoe = c.rhoe;
307 cprime.rhovel[0] = -c.rhovel[1];
308 cprime.rhovel[1] = c.rhovel[0];
309 cprime.rhovel[2] = c.rhovel[2];
310 return cprime;
311 }
312
313 template<class Tcons>
314 inline constexpr Tcons z_to_x(const Tcons c) {
315 Tcons cprime;
316 cprime.rho = c.rho;
317 cprime.rhoe = c.rhoe;
318 cprime.rhovel[0] = c.rhovel[2];
319 cprime.rhovel[1] = c.rhovel[1];
320 cprime.rhovel[2] = -c.rhovel[0];
321 return cprime;
322 }
323
324 template<class Tcons>
325 inline constexpr Tcons x_to_z(const Tcons c) {
326 Tcons cprime;
327 cprime.rho = c.rho;
328 cprime.rhoe = c.rhoe;
329 cprime.rhovel[0] = -c.rhovel[2];
330 cprime.rhovel[1] = c.rhovel[1];
331 cprime.rhovel[2] = c.rhovel[0];
332 return cprime;
333 }
334
335 template<class Tcons>
336 inline constexpr Tcons invert_axis(const Tcons c) {
337 Tcons cprime;
338 cprime.rho = c.rho;
339 cprime.rhoe = c.rhoe;
340 cprime.rhovel[0] = -c.rhovel[0];
341 cprime.rhovel[1] = -c.rhovel[1];
342 cprime.rhovel[2] = -c.rhovel[2];
343 return cprime;
344 }
345
346 // Axis-transform helpers for PrimState, mirroring y_to_x/z_to_x/invert_axis above.
347 // Riemann solvers take primitive states directly (see riemann_hll.hpp etc.), so these
348 // are applied to the inputs; the flux they return is a ConsState and is rotated back
349 // with the untransformed x_to_y/x_to_z/invert_axis.
350 template<class Tprim>
351 inline constexpr Tprim prim_y_to_x(const Tprim p) {
352 Tprim pprime;
353 pprime.rho = p.rho;
354 pprime.press = p.press;
355 pprime.vel[0] = p.vel[1];
356 pprime.vel[1] = -p.vel[0];
357 pprime.vel[2] = p.vel[2];
358 return pprime;
359 }
360
361 template<class Tprim>
362 inline constexpr Tprim prim_z_to_x(const Tprim p) {
363 Tprim pprime;
364 pprime.rho = p.rho;
365 pprime.press = p.press;
366 pprime.vel[0] = p.vel[2];
367 pprime.vel[1] = p.vel[1];
368 pprime.vel[2] = -p.vel[0];
369 return pprime;
370 }
371
372 template<class Tprim>
373 inline constexpr Tprim prim_invert_axis(const Tprim p) {
374 Tprim pprime;
375 pprime.rho = p.rho;
376 pprime.press = p.press;
377 pprime.vel[0] = -p.vel[0];
378 pprime.vel[1] = -p.vel[1];
379 pprime.vel[2] = -p.vel[2];
380 return pprime;
381 }
382
383 template<class VecType>
385 using Tvec = VecType;
386 using Tscal = shambase::VecComponent<Tvec>;
387
388 Tscal rho{};
389 Tvec rhovel{};
390
391 const DustConsState &operator+=(const DustConsState &);
392 const DustConsState &operator-=(const DustConsState &);
393 const DustConsState &operator*=(const Tscal);
394 };
395
396 template<class VecType>
398 using Tvec = VecType;
399 using Tscal = shambase::VecComponent<Tvec>;
400 Tscal rho{};
401 Tvec vel{};
402 };
403
404 template<class Tvec>
405 const DustConsState<Tvec> &DustConsState<Tvec>::operator+=(const DustConsState<Tvec> &d_cst) {
406 rho += d_cst.rho;
407 rhovel += d_cst.rhovel;
408 return *this;
409 }
410
411 template<class Tvec>
412 const DustConsState<Tvec> operator+(
413 const DustConsState<Tvec> &lhs, const DustConsState<Tvec> &rhs) {
414 return DustConsState<Tvec>(lhs) += rhs;
415 }
416
417 template<class Tvec>
418 const DustConsState<Tvec> &DustConsState<Tvec>::operator-=(const DustConsState<Tvec> &d_cst) {
419 rho -= d_cst.rho;
420 rhovel -= d_cst.rhovel;
421 return *this;
422 }
423
424 template<class Tvec>
425 const DustConsState<Tvec> operator-(
426 const DustConsState<Tvec> &lhs, const DustConsState<Tvec> &rhs) {
427 return DustConsState<Tvec>(lhs) -= rhs;
428 }
429
430 template<class Tvec>
431 const DustConsState<Tvec> &DustConsState<Tvec>::operator*=(
432 const typename DustConsState<Tvec>::Tscal factor) {
433 rho *= factor;
434 rhovel *= factor;
435 return *this;
436 }
437
438 template<class Tvec>
439 const DustConsState<Tvec> operator*(
440 const DustConsState<Tvec> &lhs, const typename DustConsState<Tvec>::Tscal factor) {
441 return DustConsState<Tvec>(lhs) *= factor;
442 }
443
444 template<class Tvec>
445 const DustConsState<Tvec> operator*(
446 const typename DustConsState<Tvec>::Tscal factor, const DustConsState<Tvec> &rhs) {
447 return DustConsState<Tvec>(rhs) *= factor;
448 }
449
450 template<class VecType>
451 struct DustFluxes {
452 using Tvec = VecType;
453 using Tscal = shambase::VecComponent<Tvec>;
454 std::array<DustConsState<Tvec>, 3> f;
455 };
456
457 template<class Tvec>
458 inline constexpr DustConsState<Tvec> d_prim_to_cons(const DustPrimState<Tvec> d_prim) {
459 DustConsState<Tvec> d_cons;
460 d_cons.rho = d_prim.rho;
461 d_cons.rhovel = (d_prim.vel * d_prim.rho);
462 return d_cons;
463 }
464
465 template<class Tvec>
466 inline constexpr DustPrimState<Tvec> d_cons_to_prim(const DustConsState<Tvec> d_cons) {
467 DustPrimState<Tvec> d_prim;
468 d_prim.rho = d_cons.rho;
469 d_prim.vel = (d_cons.rhovel * (1 / d_cons.rho));
470 return d_prim;
471 }
472
479 template<class Tvec>
481 const DustPrimState<Tvec> d_prim, Tvec n, typename DustPrimState<Tvec>::Tscal vn) {
482 DustConsState<Tvec> d_flux;
483 d_flux.rho = d_prim.rho * vn;
484 d_flux.rhovel = d_prim.vel * (d_prim.rho * vn);
485 return d_flux;
486 }
487
493 template<class Tvec>
494 inline constexpr DustConsState<Tvec> d_hydro_flux_n(const DustPrimState<Tvec> d_prim, Tvec n) {
495 const auto vn = n[0] * d_prim.vel[0] + n[1] * d_prim.vel[1] + n[2] * d_prim.vel[2];
496 return d_hydro_flux_n(d_prim, n, vn);
497 }
498
499 template<class Tcons>
500 inline constexpr Tcons d_x_to_y(const Tcons c) {
501 Tcons d_cst;
502 d_cst.rho = c.rho;
503 d_cst.rhovel[0] = -c.rhovel[1];
504 d_cst.rhovel[1] = c.rhovel[0];
505 d_cst.rhovel[2] = c.rhovel[2];
506
507 return d_cst;
508 }
509
510 template<class Tcons>
511 inline constexpr Tcons d_y_to_x(const Tcons c) {
512 Tcons d_cst;
513 d_cst.rho = c.rho;
514 d_cst.rhovel[0] = c.rhovel[1];
515 d_cst.rhovel[1] = -c.rhovel[0];
516 d_cst.rhovel[2] = c.rhovel[2];
517 return d_cst;
518 }
519
520 template<class Tcons>
521 inline constexpr Tcons d_x_to_z(const Tcons c) {
522 Tcons d_cst;
523 d_cst.rho = c.rho;
524 d_cst.rhovel[0] = -c.rhovel[2];
525 d_cst.rhovel[1] = c.rhovel[1];
526 d_cst.rhovel[2] = c.rhovel[0];
527 return d_cst;
528 }
529
530 template<class Tcons>
531 inline constexpr Tcons d_z_to_x(const Tcons c) {
532 Tcons d_cst;
533 d_cst.rho = c.rho;
534 d_cst.rhovel[0] = c.rhovel[2];
535 d_cst.rhovel[1] = c.rhovel[1];
536 d_cst.rhovel[2] = -c.rhovel[0];
537 return d_cst;
538 }
539
540 template<class Tcons>
541 inline constexpr Tcons d_invert_axis(const Tcons c) {
542 Tcons d_cst;
543 d_cst.rho = c.rho;
544 d_cst.rhovel = -(c.rhovel);
545 return d_cst;
546 }
547
548 // Axis-transform helpers for DustPrimState, mirroring d_y_to_x/d_z_to_x/d_invert_axis
549 // above. Dust Riemann solvers take primitive states directly, so these are applied to
550 // the inputs; the flux they return is a DustConsState and is rotated back with the
551 // untransformed d_x_to_y/d_x_to_z/d_invert_axis.
552 template<class Tprim>
553 inline constexpr Tprim d_prim_y_to_x(const Tprim p) {
554 Tprim pprime;
555 pprime.rho = p.rho;
556 pprime.vel[0] = p.vel[1];
557 pprime.vel[1] = -p.vel[0];
558 pprime.vel[2] = p.vel[2];
559 return pprime;
560 }
561
562 template<class Tprim>
563 inline constexpr Tprim d_prim_z_to_x(const Tprim p) {
564 Tprim pprime;
565 pprime.rho = p.rho;
566 pprime.vel[0] = p.vel[2];
567 pprime.vel[1] = p.vel[1];
568 pprime.vel[2] = -p.vel[0];
569 return pprime;
570 }
571
572 template<class Tprim>
573 inline constexpr Tprim d_prim_invert_axis(const Tprim p) {
574 Tprim pprime;
575 pprime.rho = p.rho;
576 pprime.vel = -(p.vel);
577 return pprime;
578 }
579
583 template<class VecType>
585 using Tvec = VecType;
586 using Tscal = shambase::VecComponent<Tvec>;
587 using Tprim = PrimState<Tvec>;
588 using Tcons = ConsState<Tvec>;
589
590 Tscal m_gamma; // need a different name than the methods below
591
592 Tprim cons_to_prim(Tcons c) const { return shammath::cons_to_prim(c, m_gamma); }
593 Tcons prim_to_cons(Tprim p) const { return shammath::prim_to_cons(p, m_gamma); }
594 Tscal sound_speed(Tprim p) const { return shammath::sound_speed(p, m_gamma); }
595 Tscal vn(Tprim p, Tvec n) const { return sham::dot(p.vel, n); }
596 Tcons flux(Tprim p, Tvec n, Tscal vn) const {
597 return shammath::hydro_flux_n(p, n, vn, m_gamma);
598 }
599 Tcons flux(Tprim p, Tvec n) const { return shammath::hydro_flux_n(p, n, m_gamma); }
600 Tscal gamma() const { return m_gamma; }
601 };
602
605
612 template<class VecType>
614 using Tvec = VecType;
615 using Tscal = shambase::VecComponent<Tvec>;
616 using Tprim = DustPrimState<Tvec>;
617 using Tcons = DustConsState<Tvec>;
618
619 Tprim cons_to_prim(Tcons c) const { return shammath::d_cons_to_prim(c); }
620 Tcons prim_to_cons(Tprim p) const { return shammath::d_prim_to_cons(p); }
621 Tscal vn(Tprim p, Tvec n) const { return sham::dot(p.vel, n); }
622 Tcons flux(Tprim p, Tvec n, Tscal vn) const { return shammath::d_hydro_flux_n(p, n, vn); }
623 Tcons flux(Tprim p, Tvec n) const { return shammath::d_hydro_flux_n(p, n); }
624 };
625
627
628} // namespace shammath
The flux operations a dust (pressureless) Riemann solver needs, so that solvers (see riemann_dust_hll...
A FluidStateSpec that also exposes the adiabatic index, for solvers (e.g. HLLC) that need gamma direc...
An equation of state paired with the flux/wave-speed operations a Riemann solver needs,...
True if T exposes a single, state-independent adiabatic index via gamma().
True if T exposes a per-primitive-state adiabatic index via gamma(prim).
Namespace for internal details of the logs module.
namespace for math utility
Definition AABB.hpp:26
constexpr DustConsState< Tvec > d_hydro_flux_n(const DustPrimState< Tvec > d_prim, Tvec n, typename DustPrimState< Tvec >::Tscal vn)
Pressureless (dust) flux across a face of normal n, given a precomputed normal velocity vn = dot(d_pr...
constexpr std::pair< typename FSpec::Tscal, typename FSpec::Tscal > get_adiabatic_index_lr(const FSpec &fspec, const typename FSpec::Tprim &prim_l, const typename FSpec::Tprim &prim_r)
Read the left/right adiabatic indices a HLLC-style solver should use for a given L/R pair.
constexpr ConsState< Tvec > hydro_flux_n(const PrimState< Tvec > prim, Tvec n, typename PrimState< Tvec >::Tscal vn, typename PrimState< Tvec >::Tscal gamma)
Euler flux across a face of normal n, given a precomputed normal velocity vn = dot(prim....
FluidStateSpec implementation for an ideal (adiabatic) gas equation of state.
cons_to_prim/prim_to_cons/vn/flux wrapper for a pressureless (dust) fluid. Unlike FluidStateAdiabatic...