Skip to content

Commit e0745b5

Browse files
AstroKrielpre-commit-ci[bot]
AstroKriel
and
pre-commit-ci[bot]
authored
HLLD solver: Fix pressure scaling for low-Mach flows (#1372)
### Description This PR updates the HLLD solver to apply the same low-Mach pressure correction already implemented in the HLLC solver. ### Related issues Fixes the remaining issues detailed in #528 ### Checklist - [x] I have added a description (see above). - [x] I have added a link to any related issues (if applicable; see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [x] I have added tests for any new physics that this PR adds to the code. - [ ] *(For quokka-astro org members)* I have manually triggered the GPU tests with the magic comment `/azp run`. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
1 parent e61113d commit e0745b5

4 files changed

Lines changed: 410 additions & 37 deletions

File tree

‎inputs/MHDBalsaraVortex.in‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
# *****************************************************************
2+
# Problem size and geometry
3+
# *****************************************************************
4+
geometry.prob_lo = -5.0 -5.0 -5.0
5+
geometry.prob_hi = 5.0 5.0 5.0
6+
geometry.is_periodic = 1 1 1
7+
8+
# *****************************************************************
9+
# VERBOSITY
10+
# *****************************************************************
11+
amr.v = 1
12+
13+
# *****************************************************************
14+
# Resolution and refinement
15+
# *****************************************************************
16+
amr.n_cell = 32 32 8
17+
amr.max_level = 0
18+
amr.max_grid_size = 32
19+
amr.blocking_factor_x = 32
20+
amr.blocking_factor_y = 32
21+
amr.blocking_factor_z = 8
22+
23+
do_reflux = 0
24+
do_subcycle = 0
25+
26+
plotfile_interval = 10
27+
plotfile_prefix = balsara_vortex_plt
28+
checkpoint_interval = -1
29+
cfl = 0.3
30+
max_timesteps = 10000
31+
32+
hydro.rk_integrator_order = 2
33+
hydro.reconstruction_order = 5
34+
hydro.use_dual_energy = 0
35+
36+
setup.vortex_Mach = 0.1
37+
setup.vortex_b_magn = 0.1
38+
setup.advection = 1
39+
setup.num_orbits = 3

‎src/hydro/HLLD.hpp‎

Lines changed: 67 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,8 @@ namespace quokka::Riemann
1515
{
1616
constexpr double DELTA = 1.0e-4;
1717

18-
// HLLD solver following Miyoshi and Kusano (2005), hereafter MK5.
18+
// HLLD solver based on Miyoshi & Kusano (2005), hereafter MK5.
19+
// Corrections based on Minoshima & Miyoshi (2021), hereafter MM21.
1920
template <typename problem_t, int N_scalars, int N_mscalars, int fluxdim>
2021
AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_mscalars> const &sL, quokka::HydroState<N_scalars, N_mscalars> const &sR,
2122
const double gamma, const double bx, const double perp_v_jump)
@@ -40,16 +41,20 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
4041
ConsHydro1D<N_scalars> u_star_R{};
4142

4243
// frequently used term
43-
double const bx_sq = SQUARE(bx);
44+
const double vel_magn_sq_L = SQUARE(sL.u) + (SQUARE(sL.v) + SQUARE(sL.w));
45+
const double vel_magn_sq_R = SQUARE(sR.u) + (SQUARE(sR.v) + SQUARE(sR.w));
46+
const double bx_sq = SQUARE(bx);
47+
const double b_magn_sq_L = bx_sq + (SQUARE(sL.by) + SQUARE(sL.bz));
48+
const double b_magn_sq_R = bx_sq + (SQUARE(sR.by) + SQUARE(sR.bz));
4449

4550
// compute L/R states for select conserved variables
4651
// (group transverse vector components for floating-point associativity symmetry)
4752
// magnetic pressure
48-
double const pb_L = 0.5 * (bx_sq + (SQUARE(sL.by) + SQUARE(sL.bz)));
49-
double const pb_R = 0.5 * (bx_sq + (SQUARE(sR.by) + SQUARE(sR.bz)));
53+
const double pb_L = 0.5 * b_magn_sq_L;
54+
const double pb_R = 0.5 * b_magn_sq_R;
5055
// kinetic energy
51-
double const ke_L = 0.5 * sL.rho * (SQUARE(sL.u) + (SQUARE(sL.v) + SQUARE(sL.w)));
52-
double const ke_R = 0.5 * sR.rho * (SQUARE(sR.u) + (SQUARE(sR.v) + SQUARE(sR.w)));
56+
const double ke_L = 0.5 * (sL.rho * vel_magn_sq_L);
57+
const double ke_R = 0.5 * (sR.rho * vel_magn_sq_R);
5358
// set left conserved states
5459
u_L.rho = sL.rho;
5560
u_L.mx = sL.u * u_L.rho;
@@ -76,18 +81,21 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
7681

7782
//--- Step 2. Compute L & R wave speeds according to MK5, eqn. (67)
7883

79-
const double fspd_L = FastMagnetoSonicSpeed(gamma, sL, bx);
80-
const double fspd_R = FastMagnetoSonicSpeed(gamma, sR, bx);
81-
spds[0] = std::min(sL.u - fspd_L, sR.u - fspd_R);
82-
spds[4] = std::max(sL.u + fspd_L, sR.u + fspd_R);
83-
const double fspd_m = -std::min(0.0, spds[0]);
84-
const double fspd_p = std::max(0.0, spds[4]);
84+
const double spd_fms_L = FastMagnetoSonicSpeed(gamma, sL, bx);
85+
const double spd_fms_R = FastMagnetoSonicSpeed(gamma, sR, bx);
86+
const double spd_fms_max = std::max(spd_fms_L, spd_fms_R);
87+
const double u_min = std::min(sL.u, sR.u);
88+
const double u_max = std::max(sL.u, sR.u);
89+
spds[0] = std::min(0.0, u_min - spd_fms_max);
90+
spds[4] = std::max(0.0, u_max + spd_fms_max);
91+
const double fspd_m = -spds[0];
92+
const double fspd_p = spds[4];
8593

8694
//--- Step 3. Compute L/R fluxes
8795

8896
// total pressure
89-
double ptot_L = sL.P + pb_L;
90-
double ptot_R = sR.P + pb_R;
97+
const double ptot_L = sL.P + pb_L;
98+
const double ptot_R = sR.P + pb_R;
9199
// fluxes on the left side of the interface
92100
f_L.rho = u_L.mx;
93101
f_L.mx = u_L.mx * sL.u + ptot_L - bx_sq;
@@ -116,25 +124,23 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
116124
//--- Step 4. Compute middle and Alfven wave speeds
117125

118126
// MK5: S_i - u_i (for i=L or R)
119-
double siui_L = spds[0] - sL.u;
120-
double siui_R = spds[4] - sR.u;
127+
const double siui_L = spds[0] - sL.u;
128+
const double siui_R = spds[4] - sR.u;
121129
// carbuncle detector
122-
const double max_spd = std::max(fspd_L, fspd_R);
123-
const double para_v_jump = sL.u - sR.u; // negative -> compression
124-
double theta = 1.0;
130+
const double para_v_jump = sR.u - sL.u; // negative -> compression
125131
// tp := shock anisotropy, clamped to [0, 1], with theta = tp^4
126-
const double denom_tp = std::max(1e-14, max_spd - std::min(perp_v_jump, 0.0));
127-
double tp = (max_spd - std::min(para_v_jump, 0.0)) / denom_tp;
132+
const double denom_tp = std::max(1e-14, spd_fms_max - std::min(perp_v_jump, 0.0));
133+
double tp = (spd_fms_max - std::min(para_v_jump, 0.0)) / denom_tp;
128134
tp = amrex::Clamp(tp, 0.0, 1.0);
129-
theta = SQUARE(SQUARE(tp));
130-
// modified middle speed S_M from MK5 eqn (38)
135+
const double theta = SQUARE(SQUARE(tp));
136+
// modified middle speed S_M from MK5 eqn 38 with theta from MM21 eqn 9
131137
const double sm_denom = (siui_R * u_R.rho - siui_L * u_L.rho);
132138
spds[2] = (siui_R * u_R.mx - siui_L * u_L.mx + theta * (ptot_L - ptot_R)) / sm_denom;
133139
// S_i - S_M (for i=L or R)
134-
double sism_L = spds[0] - spds[2];
135-
double sism_R = spds[4] - spds[2];
136-
double sism_inv_L = 1.0 / sism_L;
137-
double sism_inv_R = 1.0 / sism_R;
140+
const double sism_L = spds[0] - spds[2];
141+
const double sism_R = spds[4] - spds[2];
142+
const double sism_inv_L = 1.0 / sism_L;
143+
const double sism_inv_R = 1.0 / sism_R;
138144
// MK5: rho_i from eqn (43)
139145
u_star_L.rho = u_L.rho * siui_L * sism_inv_L;
140146
u_star_R.rho = u_R.rho * siui_R * sism_inv_R;
@@ -145,21 +151,45 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
145151
u_star_R.scalar[n] = u_R.scalar[n] * siui_R * sism_inv_R;
146152
}
147153

148-
double u_star_rho_inv_L = 1.0 / u_star_L.rho;
149-
double u_star_rho_inv_R = 1.0 / u_star_R.rho;
150-
double rho_sqrt_L = std::sqrt(u_star_L.rho);
151-
double rho_sqrt_R = std::sqrt(u_star_R.rho);
154+
const double u_star_rho_inv_L = 1.0 / u_star_L.rho;
155+
const double u_star_rho_inv_R = 1.0 / u_star_R.rho;
156+
const double rho_sqrt_L = std::sqrt(u_star_L.rho);
157+
const double rho_sqrt_R = std::sqrt(u_star_R.rho);
152158
// MK5: eqn (51)
153159
spds[1] = spds[2] - std::abs(bx) / rho_sqrt_L;
154160
spds[3] = spds[2] + std::abs(bx) / rho_sqrt_R;
155161

156162
//--- Step 5. Compute intermediate states
157163

158-
// compute total pressure
159-
// MK5: eqn (41) can be calculated (more explicitly) via eqn (23)
160-
double ptot_star_L = ptot_L + u_L.rho * siui_L * (spds[2] - sL.u);
161-
double ptot_star_R = ptot_R + u_R.rho * siui_R * (spds[2] - sR.u);
162-
double const ptot_star = 0.5 * (ptot_star_L + ptot_star_R);
164+
// // total pressure (no correction)
165+
// // MK5: eqn (41) can be calculated (more explicitly) via eqn (23)
166+
// double ptot_star_L = ptot_L + u_L.rho * siui_L * (spds[2] - sL.u);
167+
// double ptot_star_R = ptot_R + u_R.rho * siui_R * (spds[2] - sR.u);
168+
// double const ptot_star = 0.5 * (ptot_star_L + ptot_star_R);
169+
170+
// total pressure (w/ MM21 low-Mach correction)
171+
// MM21 eqn 8
172+
const double spd_ca_sq_L = b_magn_sq_L / sL.rho;
173+
const double spd_ca_sq_R = b_magn_sq_R / sR.rho;
174+
const double spd_cax_sq_L = bx_sq / sL.rho;
175+
const double spd_cax_sq_R = bx_sq / sR.rho;
176+
// cu := modified fast magnetosonic speed (MM21 eqn 14)
177+
const double cu_sqrt_arg_L = std::max(0.0, SQUARE(spd_ca_sq_L + vel_magn_sq_L) - 4.0 * vel_magn_sq_L * spd_cax_sq_L);
178+
const double cu_sqrt_arg_R = std::max(0.0, SQUARE(spd_ca_sq_R + vel_magn_sq_R) - 4.0 * vel_magn_sq_R * spd_cax_sq_R);
179+
const double spd_cu_sq_L = 0.5 * ((spd_ca_sq_L + vel_magn_sq_L) + std::sqrt(cu_sqrt_arg_L));
180+
const double spd_cu_sq_R = 0.5 * ((spd_ca_sq_R + vel_magn_sq_R) + std::sqrt(cu_sqrt_arg_R));
181+
const double spd_cu_L = std::sqrt(std::max(0.0, spd_cu_sq_L));
182+
const double spd_cu_R = std::sqrt(std::max(0.0, spd_cu_sq_R));
183+
const double spd_cu_max = std::max(spd_cu_L, spd_cu_R);
184+
// limiter (MM21 eqn 16)
185+
const double chi = std::min(1.0, spd_cu_max / std::max(1e-14, spd_fms_max));
186+
const double phi = chi * (2.0 - chi);
187+
// corrected total star pressure (MM21 eqn 15)
188+
const double ptot_star_numer_term1 = (siui_R * u_R.rho) * ptot_L;
189+
const double ptot_star_numer_term2 = (siui_L * u_L.rho) * ptot_R;
190+
const double ptot_star_numer_term3 = phi * (u_L.rho * u_R.rho) * (siui_R * siui_L) * (sR.u - sL.u);
191+
const double ptot_star_numer = ptot_star_numer_term1 - ptot_star_numer_term2 + ptot_star_numer_term3;
192+
const double ptot_star = ptot_star_numer / sm_denom;
163193

164194
// MK5: u_L^(star, dstar) from, eqn (39)
165195
u_star_L.mx = u_star_L.rho * spds[2];
@@ -205,7 +235,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
205235
}
206236
// vec(v_R^star) dot vec(b_R^star)
207237
// group transverse momenta-components for floating-point associativity
208-
double vb_star_R = (u_star_R.mx * bx + (u_star_R.my * u_star_R.by + u_star_R.mz * u_star_R.bz)) * u_star_rho_inv_R;
238+
const double vb_star_R = (u_star_R.mx * bx + (u_star_R.my * u_star_R.by + u_star_R.mz * u_star_R.bz)) * u_star_rho_inv_R;
209239
// MK5: eqn (48)
210240
u_star_R.E = (siui_R * u_R.E - ptot_R * sR.u + ptot_star * spds[2] + bx * (sR.u * bx + (sR.v * u_R.by + sR.w * u_R.bz) - vb_star_R)) * sism_inv_R;
211241

‎src/problems/MHDBalsaraVortex/CMakeLists.txt‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
if(AMReX_SPACEDIM EQUAL 3)
2+
set(JOB_NAME MHDBalsaraVortex)
3+
4+
add_executable(${JOB_NAME} test${JOB_NAME}.cpp ${QuokkaObjSources})
5+
6+
if(AMReX_GPU_BACKEND MATCHES "CUDA")
7+
setup_target_for_cuda_compilation(${JOB_NAME})
8+
endif()
9+
10+
add_test(
11+
NAME ${JOB_NAME}
12+
COMMAND ${JOB_NAME} ../inputs/${JOB_NAME}.in
13+
WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
14+
endif()

0 commit comments

Comments
 (0)