Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
riemann_dust.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
19#include "shambackends/math.hpp"
20#include "shambackends/typeAliasVec.hpp"
21#include "shambackends/vec.hpp"
22#include <algorithm>
23#include <array>
24#include <cmath>
25#include <iostream>
26
27namespace shammath {
28 template<class Tvec_>
29 struct DustConsState {
30 using Tvec = Tvec_;
31 using Tscal = shambase::VecComponent<Tvec>;
32
33 Tscal rho{};
34 Tvec rhovel{};
35
36 const DustConsState &operator+=(const DustConsState &);
37 const DustConsState &operator-=(const DustConsState &);
38 const DustConsState &operator*=(const Tscal);
39 };
40
41 template<class Tvec_>
42 struct DustPrimState {
43 using Tvec = Tvec_;
44 using Tscal = shambase::VecComponent<Tvec>;
45 Tscal rho{};
46 Tvec vel{};
47 };
48
49 template<class Tvec>
50 const DustConsState<Tvec> &DustConsState<Tvec>::operator+=(const DustConsState<Tvec> &d_cst) {
51 rho += d_cst.rho;
52 rhovel += d_cst.rhovel;
53 return *this;
54 }
55
56 template<class Tvec>
57 const DustConsState<Tvec> operator+(
58 const DustConsState<Tvec> &lhs, const DustConsState<Tvec> &rhs) {
59 return DustConsState<Tvec>(lhs) += rhs;
60 }
61
62 template<class Tvec>
63 const DustConsState<Tvec> &DustConsState<Tvec>::operator-=(const DustConsState<Tvec> &d_cst) {
64 rho -= d_cst.rho;
65 rhovel -= d_cst.rhovel;
66 return *this;
67 }
68
69 template<class Tvec>
70 const DustConsState<Tvec> operator-(
71 const DustConsState<Tvec> &lhs, const DustConsState<Tvec> &rhs) {
72 return DustConsState<Tvec>(lhs) -= rhs;
73 }
74
75 template<class Tvec>
76 const DustConsState<Tvec> &DustConsState<Tvec>::operator*=(
77 const typename DustConsState<Tvec>::Tscal factor) {
78 rho *= factor;
79 rhovel *= factor;
80 return *this;
81 }
82
83 template<class Tvec>
84 const DustConsState<Tvec> operator*(
85 const DustConsState<Tvec> &lhs, const typename DustConsState<Tvec>::Tscal factor) {
86 return DustConsState<Tvec>(lhs) *= factor;
87 }
88
89 template<class Tvec>
90 const DustConsState<Tvec> operator*(
91 const typename DustConsState<Tvec>::Tscal factor, const DustConsState<Tvec> &rhs) {
92 return DustConsState<Tvec>(rhs) *= factor;
93 }
94
95 template<class Tvec_>
96 struct DustFluxes {
97 using Tvec = Tvec_;
98 using Tscal = shambase::VecComponent<Tvec>;
99 std::array<DustConsState<Tvec>, 3> F;
100 };
101
102 template<class Tvec>
103 inline constexpr DustConsState<Tvec> d_prim_to_cons(const DustPrimState<Tvec> d_prim) {
104 DustConsState<Tvec> d_cons;
105 d_cons.rho = d_prim.rho;
106 d_cons.rhovel = (d_prim.vel * d_prim.rho);
107 return d_cons;
108 }
109
110 template<class Tvec>
111 inline constexpr DustPrimState<Tvec> d_cons_to_prim(const DustConsState<Tvec> d_cons) {
112 DustPrimState<Tvec> d_prim;
113 d_prim.rho = d_cons.rho;
114 d_prim.vel = (d_cons.rhovel * (1 / d_cons.rho));
115 return d_prim;
116 }
117
118 template<class Tvec>
119 inline constexpr DustConsState<Tvec> d_hydro_flux_x(const DustConsState<Tvec> d_cons) {
120 DustConsState<Tvec> d_flux;
121 const DustPrimState<Tvec> d_prim = d_cons_to_prim<Tvec>(d_cons);
122 const typename DustConsState<Tvec>::Tscal x_vel{d_prim.vel[0]};
123 d_flux.rho = d_cons.rhovel[0];
124 d_flux.rhovel = d_prim.vel * (d_cons.rho * x_vel);
125 return d_flux;
126 }
127
128 template<class Tcons>
129 inline constexpr Tcons d_x_to_y(const Tcons c) {
130 Tcons d_cst;
131 d_cst.rho = c.rho;
132 d_cst.rhovel[0] = -c.rhovel[1];
133 d_cst.rhovel[1] = c.rhovel[0];
134 d_cst.rhovel[2] = c.rhovel[2];
135
136 return d_cst;
137 }
138
139 template<class Tcons>
140 inline constexpr Tcons d_y_to_x(const Tcons c) {
141 Tcons d_cst;
142 d_cst.rho = c.rho;
143 d_cst.rhovel[0] = c.rhovel[1];
144 d_cst.rhovel[1] = -c.rhovel[0];
145 d_cst.rhovel[2] = c.rhovel[2];
146 return d_cst;
147 }
148
149 template<class Tcons>
150 inline constexpr Tcons d_x_to_z(const Tcons c) {
151 Tcons d_cst;
152 d_cst.rho = c.rho;
153 d_cst.rhovel[0] = -c.rhovel[2];
154 d_cst.rhovel[1] = c.rhovel[1];
155 d_cst.rhovel[2] = c.rhovel[0];
156 return d_cst;
157 }
158
159 template<class Tcons>
160 inline constexpr Tcons d_z_to_x(const Tcons c) {
161 Tcons d_cst;
162 d_cst.rho = c.rho;
163 d_cst.rhovel[0] = c.rhovel[2];
164 d_cst.rhovel[1] = c.rhovel[1];
165 d_cst.rhovel[2] = -c.rhovel[0];
166 return d_cst;
167 }
168
169 template<class Tcons>
170 inline constexpr Tcons d_invert_axis(const Tcons c) {
171 Tcons d_cst;
172 d_cst.rho = c.rho;
173 d_cst.rhovel = -(c.rhovel);
174 return d_cst;
175 }
176
177 // Krapp et al. 2024, A Fast second-order solver for stiff multifluid dust and gas hydrodynamics
178 // Appendice E
179 template<class Tcons>
180 inline constexpr auto d_hll_flux_x(Tcons cL, Tcons cR) {
181 Tcons d_flux;
182 const auto d_primL = d_cons_to_prim(cL);
183 const auto d_primR = d_cons_to_prim(cR);
184
185 const auto S = sham::max(sham::abs(d_primL.vel[0]), sham::abs(d_primR.vel[0]));
186
187 const auto fL = d_hydro_flux_x(cL);
188 const auto fR = d_hydro_flux_x(cR);
189
190 return 0.5 * ((fL + fR) - S * (cR - cL));
191 }
192
193 // Huang & Bai, 2022 ,A Multifluid Dust Module in Athena++: Algorithms and Numerical Tests
194 // Equation (32)
195 template<class Tcons>
196 inline constexpr auto huang_bai_flux_x(Tcons cL, Tcons cR) {
197 Tcons d_flux;
198 const auto d_primL = d_cons_to_prim(cL);
199 const auto d_primR = d_cons_to_prim(cR);
200
201 const auto fL = d_hydro_flux_x(cL);
202 const auto fR = d_hydro_flux_x(cR);
203
204 if (d_primL.vel[0] > 0 && d_primR.vel[0] > 0)
205 d_flux = fL;
206 else if (d_primL.vel[0] < 0 && d_primR.vel[0] < 0)
207 d_flux = fR;
208 else if (d_primL.vel[0] < 0 && d_primR.vel[0] > 0)
209 d_flux *= 0;
210 else if (d_primL.vel[0] > 0 && d_primR.vel[0] < 0)
211 d_flux = (fL + fR);
212
213 return d_flux;
214 }
215
216 template<class Tcons>
217 inline constexpr Tcons d_hll_flux_y(Tcons cL, Tcons cR) {
218 return d_x_to_y(d_hll_flux_x(d_y_to_x(cL), d_y_to_x(cR)));
219 }
220
221 template<class Tcons>
222 inline constexpr Tcons d_hll_flux_z(Tcons cL, Tcons cR) {
223 return d_x_to_z(d_hll_flux_x(d_z_to_x(cL), d_z_to_x(cR)));
224 }
225
226 template<class Tcons>
227 inline constexpr Tcons d_hll_flux_mx(Tcons cL, Tcons cR) {
228 return d_invert_axis(d_hll_flux_x(d_invert_axis(cL), d_invert_axis(cR)));
229 }
230
231 template<class Tcons>
232 inline constexpr Tcons d_hll_flux_my(Tcons cL, Tcons cR) {
233 return d_invert_axis(d_hll_flux_y(d_invert_axis(cL), d_invert_axis(cR)));
234 }
235
236 template<class Tcons>
237 inline constexpr Tcons d_hll_flux_mz(Tcons cL, Tcons cR) {
238 return d_invert_axis(d_hll_flux_z(d_invert_axis(cL), d_invert_axis(cR)));
239 }
240
241 template<class Tcons>
242 inline constexpr Tcons huang_bai_flux_y(Tcons cL, Tcons cR) {
243 return d_x_to_y(huang_bai_flux_x(d_y_to_x(cL), d_y_to_x(cR)));
244 }
245
246 template<class Tcons>
247 inline constexpr Tcons huang_bai_flux_z(Tcons cL, Tcons cR) {
248 return d_x_to_z(huang_bai_flux_x(d_z_to_x(cL), d_z_to_x(cR)));
249 }
250
251 template<class Tcons>
252 inline constexpr Tcons huang_bai_flux_mx(Tcons cL, Tcons cR) {
253 return d_invert_axis(huang_bai_flux_x(d_invert_axis(cL), d_invert_axis(cR)));
254 }
255
256 template<class Tcons>
257 inline constexpr Tcons huang_bai_flux_my(Tcons cL, Tcons cR) {
258 return d_invert_axis(huang_bai_flux_y(d_invert_axis(cL), d_invert_axis(cR)));
259 }
260
261 template<class Tcons>
262 inline constexpr Tcons huang_bai_flux_mz(Tcons cL, Tcons cR) {
263 return d_invert_axis(huang_bai_flux_z(d_invert_axis(cL), d_invert_axis(cR)));
264 }
265} // namespace shammath
shammath
namespace for math utility
Definition AABB.hpp:26
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8