Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
mhd.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
20#include "shambase/numeric_limits.hpp"
21#include "shambackends/math.hpp"
22#include "shamcomm/logs.hpp"
23#include "shammodels/sph/Solver.hpp"
24#include "shammodels/sph/SolverConfig.hpp"
25#include "shammodels/sph/math/forces.hpp"
26#include "shammodels/sph/math/q_ab.hpp"
27#include "shamphys/mhd.hpp"
28#include "shamunits/Constants.hpp"
29#include <tuple>
30
31namespace shamrock::sph::mhd {
32
33 enum MHDType { Ideal = 0, NonIdeal = 1 };
34
35 // mag tension form the Tricco 2023 formula
36 template<class Tvec, class Tscal>
37 inline Tvec B_dot_grad_W(
38 Tscal m_b,
39 Tscal rho_a_sq,
40 Tscal rho_b_sq,
41 Tvec B_a,
42 Tvec B_b,
43 Tscal omega_a,
44 Tscal omega_b,
45 Tvec nabla_Wab_ha,
46 Tvec nabla_Wab_hb,
47 Tscal mu_0) {
48
49 Tscal sub_fact_a = rho_a_sq * omega_a;
50 Tscal sub_fact_b = rho_b_sq * omega_b;
51
52 Tscal B_dot_grad_W_a = sycl::dot(B_a, nabla_Wab_ha);
53 Tscal B_dot_grad_W_b = sycl::dot(B_b, nabla_Wab_hb);
54
55 Tvec acc_a = ((B_dot_grad_W_a) *sham::inv_sat_zero(sub_fact_a)) * B_a / mu_0;
56 Tvec acc_b = ((B_dot_grad_W_b) *sham::inv_sat_zero(sub_fact_b)) * B_b / mu_0;
57
58 return -m_b * (acc_a + acc_b);
59 }
60
61 // from the Phantom paper formula
62 template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
63 inline Tvec mag_tension(
64 Tscal m_b,
65 Tvec B_a,
66 Tvec B_b,
67 Tvec nabla_Wab_ha,
68 Tvec nabla_Wab_hb,
69 Tscal sub_fact_a,
70 Tscal sub_fact_b,
71 Tscal mu_0) {
72 Tvec magnetic_tension_term;
73
74 for (int i = 0; i < 3; ++i) {
75 for (int j = 0; j < 3; ++j) {
76
77 Tscal mag_tension_a = -(1. / mu_0) * B_a[i] * B_a[j];
78 Tscal mag_tension_b = -(1. / mu_0) * B_b[i] * B_b[j];
79 Tscal acc_mag_tension_a = mag_tension_a * sham::inv_sat_zero(sub_fact_a);
80 Tscal acc_mag_tension_b = mag_tension_b * sham::inv_sat_zero(sub_fact_b);
81
82 magnetic_tension_term[i] += -m_b
83 * (acc_mag_tension_a * nabla_Wab_ha[j]
84 + acc_mag_tension_b * nabla_Wab_hb[j]);
85 }
86 }
87
88 return magnetic_tension_term;
89 }
90
91 template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
92 inline Tvec fdivB(
93 Tscal m_b,
94 Tvec B_a,
95 Tvec B_b,
96 Tvec nabla_Wab_ha,
97 Tvec nabla_Wab_hb,
98 Tscal sub_fact_a,
99 Tscal sub_fact_b,
100 Tscal mu_0) {
101
102 Tscal acc_fdivB_a = sycl::dot(B_a, nabla_Wab_ha) * sham::inv_sat_zero(sub_fact_a);
103 Tscal acc_fdivB_b = sycl::dot(B_b, nabla_Wab_hb) * sham::inv_sat_zero(sub_fact_b);
104 return -0.5 * B_a * m_b * (acc_fdivB_a + acc_fdivB_b)
105 / mu_0; // tested, this is what works best
106 }
107
108 template<class Tvec, class Tscal>
109 inline Tscal lambda_artes(
110 Tscal m_b,
111 Tscal rho_a_sq,
112 Tscal rho_b_sq,
113 Tscal vsigb,
114 Tvec B_a,
115 Tvec B_b,
116 Tscal omega_a,
117 Tscal omega_b,
118 Tscal Fab_a,
119 Tscal Fab_b) {
120
121 Tscal B_ab_sq = sycl::dot(B_a - B_b, B_a - B_b);
122
123 Tscal sub_fact_a = rho_a_sq * omega_a;
124 Tscal sub_fact_b = rho_b_sq * omega_b;
125
126 Tscal acc_a = Fab_a * sham::inv_sat_zero(sub_fact_a);
127 Tscal acc_b = Fab_b * sham::inv_sat_zero(sub_fact_b);
128
129 Tscal artres = -0.25 * m_b * vsigb * (acc_a + acc_b) * B_ab_sq;
130 return artres;
131 }
132
133 template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
134 inline Tscal dB_on_rho_induction_term(
135 Tscal m_b, Tscal rho_a_sq, Tvec B_a, Tscal omega_a, Tvec nabla_Wab_ha) {
136
137 Tscal sub_fact_a = rho_a_sq * omega_a;
138 Tscal induction_term_no_vab
139 = -sham::inv_sat_zero(sub_fact_a) * m_b * sycl::dot(B_a, nabla_Wab_ha);
140
141 return induction_term_no_vab;
142 }
143
144 template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
145 inline Tvec dB_on_rho_psi_term(
146 Tscal m_b,
147 Tscal rho_a_sq,
148 Tscal rho_b_sq,
149 Tscal psi_a,
150 Tscal psi_b,
151 Tscal omega_a,
152 Tscal omega_b,
153 Tvec nabla_Wab_ha,
154 Tvec nabla_Wab_hb) {
155
156 Tscal sub_fact_a = rho_a_sq * omega_a;
157 Tscal sub_fact_b = rho_b_sq * omega_b;
158
159 Tvec psisubterm_a = ((psi_a) *sham::inv_sat_zero(sub_fact_a)) * nabla_Wab_ha;
160 Tvec psisubterm_b = ((psi_b) *sham::inv_sat_zero(sub_fact_b)) * nabla_Wab_hb;
161
162 Tvec psiterm = -m_b * (psisubterm_a + psisubterm_a);
163
164 return psiterm;
165 }
166
167 template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
168 inline Tscal dpsi_on_ch_parabolic_propag(
169 Tscal m_b, Tscal rho_a, Tvec B_a, Tvec B_b, Tscal omega_a, Tvec nabla_Wab_ha, Tscal ch_a) {
170
171 Tscal sub_fact_a = rho_a * omega_a;
172 Tvec B_ab = (B_a - B_b);
173
174 Tscal divB_a = -(1. * sham::inv_sat_zero(sub_fact_a)) * m_b * sycl::dot(B_ab, nabla_Wab_ha);
175
176 Tscal parabolic_propag = m_b * (ch_a * sham::inv_sat_zero(sub_fact_a))
177 * sycl::dot(B_ab, nabla_Wab_ha); //-ch_a * divB_a;
178
179 return parabolic_propag;
180 }
181
182 template<class Tvec, class Tscal, MHDType MHD_mode = Ideal>
183 inline Tscal dpsi_on_ch_parabolic_diff(
184 Tscal m_b,
185 Tscal rho_a,
186 Tvec v_ab,
187 Tscal psi_a,
188 Tscal omega_a,
189 Tvec nabla_Wab_ha,
190 Tscal ch_a) {
191
192 Tscal sub_fact_a = 2. * rho_a * omega_a * ch_a;
193
194 Tscal parabolic_diff
195 = m_b * (psi_a * sham::inv_sat_zero(sub_fact_a)) * sycl::dot(v_ab, nabla_Wab_ha);
196
197 return parabolic_diff;
198 }
199
200 template<class Kernel, class Tvec, class Tscal, MHDType MHD_mode = Ideal>
201 inline void add_to_derivs_spmhd(
202 Tscal pmass,
203 Tvec dr,
204 Tscal rab,
205 Tscal rho_a,
206 Tscal rho_a_sq,
207 Tscal omega_a_rho_a_inv,
208 Tscal rho_a_inv,
209 Tscal rho_b,
210 Tscal omega_a,
211 Tscal omega_b,
212 Tscal Fab_a,
213 Tscal Fab_b,
214 Tvec vxyz_a,
215 Tvec vxyz_b,
216 Tscal u_a,
217 Tscal u_b,
218 Tscal P_a,
219 Tscal P_b,
220 Tscal cs_a,
221 Tscal cs_b,
222 Tscal h_a,
223 Tscal h_b,
224
225 Tscal alpha_u,
226
227 Tvec B_a,
228 Tvec B_b,
229
230 Tscal psi_a,
231 Tscal psi_b,
232
233 Tscal mu_0,
234 Tscal sigma_mhd,
235
236 Tvec &dv_dt,
237 Tscal &du_dt,
238 Tvec &dB_on_rho_dt,
239 Tscal &dpsi_on_ch_dt,
240 Tscal &drho_dt,
241
242 Tvec &mag_pressure,
243 Tvec &mag_tension,
244 Tvec &gas_pressure,
245 Tvec &tensile_corr,
246
247 Tscal &psi_propag,
248 Tscal &psi_diff,
249 Tscal &psi_cons,
250
251 Tscal &u_pressure_viscous_heating) {
252
253 using namespace shamrock::sph;
254
255 Tvec v_ab = vxyz_a - vxyz_b;
256 Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
257
258 Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
259 Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
260
261 Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
262 Tscal vsig_a = shamphys::MHD_physics<Tvec, Tscal>::vsig_MHD(
263 v_ab, r_ab_unit, cs_a, B_a, rho_a, mu_0, 1., 1.);
264 Tscal vsig_b = shamphys::MHD_physics<Tvec, Tscal>::vsig_MHD(
265 v_ab, r_ab_unit, cs_a, B_b, rho_b, mu_0, 1., 1.);
266
267 Tscal dWab_a = Fab_a;
268 Tscal dWab_b = Fab_b;
269
270 Tscal sub_fact_a = rho_a_sq * omega_a;
271 Tscal sub_fact_b = rho_b * rho_b * omega_b;
272
273 Tscal v_shock_a = shamphys::MHD_physics<Tvec, Tscal>::v_shock(cs_a, B_a, rho_a, mu_0);
274 Tscal v_shock_b = shamphys::MHD_physics<Tvec, Tscal>::v_shock(cs_b, B_b, rho_b, mu_0);
275 Tscal vsig_B = shamphys::MHD_physics<Tvec, Tscal>::vsigB(v_ab, r_ab_unit);
276
277 Tscal qa_ab = shamrock::sph::q_av(rho_a, vsig_a, v_ab_r_ab);
278 Tscal qb_ab = shamrock::sph::q_av(rho_b, vsig_b, v_ab_r_ab);
279
280 Tscal AV_P_a = P_a + qa_ab;
281 Tscal AV_P_b = P_b + qb_ab;
282
283 Tvec sum_mag_tension, sum_fdivB = {0., 0., 0.};
284 Tscal sum_psi_propag, sum_psi_diff = 0.;
285
286 // dv/dt = gas_pressure_pishock + magnetic_pressure_term + sum_mag_tension + sum_fdivB;
287 // du/dt = pressure_term + viscous_heating + shock_conductivity + artificial_resistivity;
288 // d(B/rho)/dt = induction_term + psi_term + dissipation_term
289 // d(psi/ch)/dt = psi_propag + psi_diff + psi_cons (not in the SPH sum, hence added in
290 // update_derivs)
291
292 // dv/dt terms
293 sum_fdivB += fdivB(
294 pmass, B_a, B_b, r_ab_unit * dWab_a, r_ab_unit * dWab_b, sub_fact_a, sub_fact_b, mu_0);
295
296 Tvec gas_pressure_pishock = sph::sph_pressure_symetric(
297 pmass,
298 rho_a_sq,
299 rho_b * rho_b,
300 AV_P_a,
301 AV_P_b,
302 omega_a,
303 omega_b,
304 r_ab_unit * dWab_a,
305 r_ab_unit * dWab_b);
306
307 sum_mag_tension += -B_dot_grad_W(
308 pmass,
309 rho_a_sq,
310 rho_b * rho_b,
311 B_a,
312 B_b,
313 omega_a,
314 omega_b,
315 r_ab_unit * dWab_a,
316 r_ab_unit * dWab_b,
317 mu_0);
318
319 Tvec magnetic_pressure_term = (1. / (2 * mu_0))
320 * sph::sph_pressure_symetric(
321 pmass,
322 rho_a_sq,
323 rho_b * rho_b,
324 sycl::dot(B_a, B_a),
325 sycl::dot(B_b, B_b),
326 omega_a,
327 omega_b,
328 r_ab_unit * dWab_a,
329 r_ab_unit * dWab_b);
330
331 dv_dt += gas_pressure_pishock + magnetic_pressure_term + sum_mag_tension + sum_fdivB;
332 // end dv/dt terms
333
334 // get terms for debugging
335 mag_pressure += magnetic_pressure_term;
336 gas_pressure += gas_pressure_pishock;
337 mag_tension += sum_mag_tension;
338 tensile_corr += sum_fdivB;
339 // end get terms for debugging
340
341 // du/dt terms
342 u_pressure_viscous_heating = sph::duint_dt_pressure(
343 pmass, AV_P_a, omega_a_rho_a_inv * rho_a_inv, v_ab, r_ab_unit * dWab_a);
344
345 du_dt += u_pressure_viscous_heating;
346
347 du_dt += sph::lambda_shock_conductivity(
348 pmass,
349 alpha_u,
350 vsig_u,
351 u_a - u_b,
352 dWab_a * omega_a_rho_a_inv,
353 dWab_b / (rho_b * omega_b));
354
355 du_dt += lambda_artes(
356 pmass, rho_a_sq, rho_b * rho_b, vsig_B, B_a, B_b, omega_a, omega_b, Fab_a, Fab_b);
357
358 // end du/dt terms
359
360 // d(B/rho)/dt terms
361 Tscal rho_diss_term_a = Fab_a * sham::inv_sat_zero(sub_fact_a);
362 Tscal rho_diss_term_b = Fab_b * sham::inv_sat_zero(sub_fact_b);
363
364 Tvec dB_on_rho_dissipation_term
365 = 0.5 * pmass * (rho_diss_term_a + rho_diss_term_b) * (B_a - B_b) * vsig_B;
366
367 dB_on_rho_dt
368 += v_ab * dB_on_rho_induction_term(pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_b);
369
370 dB_on_rho_dt += dB_on_rho_psi_term(
371 pmass,
372 rho_a_sq,
373 rho_b * rho_b,
374 psi_a,
375 psi_b,
376 omega_a,
377 omega_b,
378 r_ab_unit * dWab_a,
379 r_ab_unit * dWab_b);
380
381 dB_on_rho_dt += dB_on_rho_dissipation_term;
382
383 // end d(B/rho)/dt terms
384
385 // d(psi/ch)/dt terms
386 sum_psi_propag = dpsi_on_ch_parabolic_propag(
387 pmass, rho_a, B_a, B_b, omega_a, r_ab_unit * dWab_a, v_shock_a);
388
389 sum_psi_diff = dpsi_on_ch_parabolic_diff(
390 pmass, rho_a, v_ab, psi_a, omega_a, r_ab_unit * dWab_a, v_shock_a);
391
392 dpsi_on_ch_dt += sum_psi_propag;
393 dpsi_on_ch_dt += sum_psi_diff;
394
395 // end d(psi/ch)/dt terms
396
397 // get debugging terms
398 psi_propag += sum_psi_propag;
399 psi_diff += sum_psi_diff;
400 // end get debugging terms
401
402 // for conservative checks
403 drho_dt += (1. / omega_a) * pmass * sycl::dot(v_ab, r_ab_unit * dWab_a);
404 }
405
406} // namespace shamrock::sph::mhd
sham::inv_sat_positive
T inv_sat_positive(T v, T minvsat=T{1e-9}, T satval=T{0.}) noexcept
inverse saturated (positive numbers only)
Definition math.hpp:841
sham::inv_sat_zero
T inv_sat_zero(T v, T satval=T{0.}) noexcept
inverse saturated (zero version)
Definition math.hpp:870
numeric_limits.hpp
q_ab.hpp
shamrock::sph::q_av
constexpr Tscal q_av(const Tscal &rho, const Tscal &vsig, const Tscal &v_scal_rhat)
eq.40
Definition q_ab.hpp:37
forces.hpp
file containing formulas for sph forces
shamphys::MHD_physics
Definition mhd.hpp:24
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8