diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 5a523daf3f6..2275ddf86c8 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -45,11 +45,11 @@ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations1D) rho = 1.0 rho_v1 = 0.1 rho_v2 = -0.2 - rho_v3 = -0.5 + rho_v3 = -0.5f0 rho_e = 50.0 B1 = 3.0 B2 = -1.2 - B3 = 0.5 + B3 = 0.5f0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) end @@ -92,9 +92,9 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations phi = atan(x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? 1.0 : 1.1691 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) + p = r > 0.5f0 ? 1.0 : 1.245 return prim2cons(SVector(rho, v1, 0.0, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations) end @@ -105,8 +105,8 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) p_over_gamma_minus_one = (rho_e - kin_en - mag_en) p = (equations.gamma - 1) * p_over_gamma_minus_one @@ -150,44 +150,44 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 p_ll = (equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll) + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll) p_rr = (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr # for convenience store v⋅B vel_dot_mag_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_mean = 0.5 * rho_avg / beta_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) # Ignore orientation since it is always "1" in 1D f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean + 0.5 * mag_norm_avg - B1_avg * B1_avg + f2 = f1 * v1_avg + p_mean + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f3 = f1 * v2_avg - B1_avg * B2_avg f4 = f1 * v3_avg - B1_avg * B3_avg f6 = 0.0 f7 = v1_avg * B2_avg - v2_avg * B1_avg f8 = v1_avg * B3_avg - v3_avg * B1_avg # total energy flux is complicated and involves the previous eight components - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + - f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5 * v1_mag_avg + + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg) return SVector(f1, f2, f3, f4, f5, f6, f7, f8) @@ -228,27 +228,27 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on orientation with specific direction averages f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = 0.0 - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -273,8 +273,8 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr = cons2prim(u_rr, equations) # Total pressure, i.e., thermal + magnetic pressures (eq. (12)) - p_tot_ll = p_ll + 0.5 * (B1_ll^2 + B2_ll^2 + B3_ll^2) - p_tot_rr = p_rr + 0.5 * (B1_rr^2 + B2_rr^2 + B3_rr^2) + p_tot_ll = p_ll + 0.5f0 * (B1_ll^2 + B2_ll^2 + B3_ll^2) + p_tot_rr = p_rr + 0.5f0 * (B1_rr^2 + B2_rr^2 + B3_rr^2) # Conserved variables rho_v1_ll = u_ll[2] @@ -502,8 +502,8 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * (rho_e - - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 - + B1 * B1 + B2 * B2 + B3 * B3)) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + + B1 * B1 + B2 * B2 + B3 * B3)) return SVector(rho, v1, v2, v3, p, B1, B2, B3) end @@ -517,11 +517,11 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2)) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2)) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) / (equations.gamma - 1) - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) / (equations.gamma - 1) - 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -541,8 +541,8 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p / (equations.gamma - 1) + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) end @@ -554,17 +554,17 @@ end @inline function pressure(u, equations::IdealGlmMhdEquations1D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2)) + 0.5f0 * (B1^2 + B2^2 + B3^2)) return p end @inline function density_pressure(u, equations::IdealGlmMhdEquations1D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2)) + 0.5f0 * (B1^2 + B2^2 + B3^2)) return rho * p end @@ -576,7 +576,7 @@ end v3 = rho_v3 / rho v_mag = sqrt(v1^2 + v2^2 + v3^2) p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2)) + (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2)) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -584,8 +584,8 @@ end b3 = B3 / sqrt_rho b_square = b1^2 + b2^2 + b3^2 - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) return c_f end @@ -611,7 +611,7 @@ as given by vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 p_ll = (equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll) + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr @@ -619,11 +619,11 @@ as given by vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 p_rr = (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr) + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) @@ -645,19 +645,19 @@ as given by H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equations (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum a_square_roe = ((2.0 - equations.gamma) * X + (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity # Ignore orientation since it is always "1" in 1D c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) return v1_roe, c_f_roe end @@ -691,13 +691,13 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::IdealGlmMhdEquations1D) - return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + return 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] end # Calculate the magnetic energy for a conservative state `cons'. # OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy @inline function energy_magnetic(cons, ::IdealGlmMhdEquations1D) - return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + return 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) end # Calculate internal energy for a conservative state `cons` diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 622410de855..bec1ee50e4e 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -51,11 +51,11 @@ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D) rho = 1.0 rho_v1 = 0.1 rho_v2 = -0.2 - rho_v3 = -0.5 + rho_v3 = -0.5f0 rho_e = 50.0 B1 = 3.0 B2 = -1.2 - B3 = 0.5 + B3 = 0.5f0 psi = 0.0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -102,10 +102,10 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations phi = atan(y_norm, x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) - v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? 1.0 : 1.1691 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) + v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin(phi) + p = r > 0.5f0 ? 1.0 : 1.245 return prim2cons(SVector(rho, v1, v2, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations) end @@ -119,9 +119,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one if orientation == 1 f1 = rho_v1 @@ -158,9 +158,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] @@ -495,36 +495,38 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 p_ll = (equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll - + 0.5f0 * psi_ll^2) p_rr = (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr - + 0.5f0 * psi_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr # for convenience store v⋅B vel_dot_mag_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_mean = 0.5 * rho_avg / beta_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean + 0.5 * mag_norm_avg - B1_avg * B1_avg + f2 = f1 * v1_avg + p_mean + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f3 = f1 * v2_avg - B1_avg * B2_avg f4 = f1 * v3_avg - B1_avg * B3_avg f6 = equations.c_h * psi_avg @@ -532,29 +534,29 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f8 = v1_avg * B3_avg - v3_avg * B1_avg f9 = equations.c_h * B1_avg # total energy flux is complicated and involves the previous eight components - psi_B1_avg = 0.5 * (B1_ll * psi_ll + B1_rr * psi_rr) - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr) + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v1_mag_avg + + 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg - equations.c_h * psi_B1_avg) else f1 = rho_mean * v2_avg f2 = f1 * v1_avg - B1_avg * B2_avg - f3 = f1 * v2_avg + p_mean + 0.5 * mag_norm_avg - B2_avg * B2_avg + f3 = f1 * v2_avg + p_mean + 0.5f0 * mag_norm_avg - B2_avg * B2_avg f4 = f1 * v3_avg - B2_avg * B3_avg f6 = v2_avg * B1_avg - v1_avg * B2_avg f7 = equations.c_h * psi_avg f8 = v2_avg * B3_avg - v3_avg * B2_avg f9 = equations.c_h * B2_avg # total energy flux is complicated and involves the previous eight components - psi_B2_avg = 0.5 * (B2_ll * psi_ll + B2_rr * psi_rr) - v2_mag_avg = 0.5 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr) + v2_mag_avg = 0.5f0 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v2_mag_avg + + 0.5f0 * v2_mag_avg + B2_avg * vel_dot_mag_avg - equations.c_h * psi_B2_avg) end @@ -598,31 +600,31 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = equations.c_h * psi_avg - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) - f9 = equations.c_h * 0.5 * (B1_ll + B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f9 = equations.c_h * 0.5f0 * (B1_ll + B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -633,20 +635,20 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat equations.c_h * (B1_ll * psi_rr + B1_rr * psi_ll))) else # orientation == 2 f1 = rho_mean * v2_avg - f2 = f1 * v1_avg - 0.5 * (B2_ll * B1_rr + B2_rr * B1_ll) + f2 = f1 * v1_avg - 0.5f0 * (B2_ll * B1_rr + B2_rr * B1_ll) f3 = f1 * v2_avg + p_avg + magnetic_square_avg - - 0.5 * (B2_ll * B2_rr + B2_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B2_ll * B3_rr + B2_rr * B3_ll) + 0.5f0 * (B2_ll * B2_rr + B2_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B2_ll * B3_rr + B2_rr * B3_ll) #f5 below - f6 = 0.5 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) + f6 = 0.5f0 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) f7 = equations.c_h * psi_avg - f8 = 0.5 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) - f9 = equations.c_h * 0.5 * (B2_ll + B2_rr) + f8 = 0.5f0 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) + f9 = equations.c_h * 0.5f0 * (B2_ll + B2_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v2_rr + p_rr * v2_ll + 0.5f0 * (+p_ll * v2_rr + p_rr * v2_ll + (v2_ll * B1_ll * B1_rr + v2_rr * B1_rr * B1_ll) + (v2_ll * B3_ll * B3_rr + v2_rr * B3_rr * B3_ll) - @@ -679,41 +681,41 @@ end # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = (f1 * v1_avg + (p_avg + magnetic_square_avg) * normal_direction[1] - - 0.5 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) + 0.5f0 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) f3 = (f1 * v2_avg + (p_avg + magnetic_square_avg) * normal_direction[2] - - 0.5 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) + 0.5f0 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) f4 = (f1 * v3_avg - - 0.5 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) + 0.5f0 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) #f5 below f6 = (equations.c_h * psi_avg * normal_direction[1] + - 0.5 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + v_dot_n_rr * B1_rr - v1_rr * B_dot_n_rr)) f7 = (equations.c_h * psi_avg * normal_direction[2] + - 0.5 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + v_dot_n_rr * B2_rr - v2_rr * B_dot_n_rr)) - f8 = +0.5 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + + f8 = +0.5f0 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + v_dot_n_rr * B3_rr - v3_rr * B_dot_n_rr) - f9 = equations.c_h * 0.5 * (B_dot_n_ll + B_dot_n_rr) + f9 = equations.c_h * 0.5f0 * (B_dot_n_ll + B_dot_n_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + 0.5f0 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + (v_dot_n_ll * B1_ll * B1_rr + v_dot_n_rr * B1_rr * B1_ll) + (v_dot_n_ll * B2_ll * B2_rr + v_dot_n_rr * B2_rr * B2_ll) + (v_dot_n_ll * B3_ll * B3_rr + v_dot_n_rr * B3_rr * B3_ll) @@ -1029,7 +1031,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * (rho_e - - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + B1 * B1 + B2 * B2 + B3 * B3 + psi * psi)) @@ -1045,11 +1047,12 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -1097,8 +1100,8 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p * equations.inv_gamma_minus_one + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5 * psi^2 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) + 0.5f0 * psi^2 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -1110,11 +1113,11 @@ end @inline function pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return p end @@ -1129,16 +1132,16 @@ end v_square = v1^2 + v2^2 + v3^2 return (equations.gamma - 1.0) * - SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) + SVector(0.5f0 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) end @inline function density_pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return rho * p end @@ -1149,9 +1152,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1159,11 +1162,11 @@ end b3 = B3 / sqrt_rho b_square = b1 * b1 + b2 * b2 + b3 * b3 if orientation == 1 # x-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) else - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) end return c_f end @@ -1174,9 +1177,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1188,8 +1191,8 @@ end b_dot_n_squared = (b1 * normal_direction[1] + b2 * normal_direction[2])^2 / norm_squared - c_f = sqrt((0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * + c_f = sqrt((0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * norm_squared) return c_f end @@ -1213,22 +1216,22 @@ as given by v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) @@ -1250,26 +1253,26 @@ as given by H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum a_square_roe = ((2.0 - equations.gamma) * X + (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) if orientation == 1 # x-direction c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v1_roe else # y-direction c_a_roe = B2_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v2_roe end @@ -1285,22 +1288,22 @@ end v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) @@ -1322,13 +1325,13 @@ end H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum a_square_roe = ((2.0 - equations.gamma) * X + (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) @@ -1339,7 +1342,7 @@ end c_a_roe = B_roe_dot_n_squared * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) vel_out_roe = (v1_roe * normal_direction[1] + v2_roe * normal_direction[2]) @@ -1378,13 +1381,13 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::IdealGlmMhdEquations2D) - return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + return 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] end # Calculate the magnetic energy for a conservative state `cons'. # OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy @inline function energy_magnetic(cons, ::IdealGlmMhdEquations2D) - return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + return 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) end # Calculate internal energy for a conservative state `cons` diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 321e501b051..679afa62b24 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -49,11 +49,11 @@ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) rho = 1.0 rho_v1 = 0.1 rho_v2 = -0.2 - rho_v3 = -0.5 + rho_v3 = -0.5f0 rho_e = 50.0 B1 = 3.0 B2 = -1.2 - B3 = 0.5 + B3 = 0.5f0 psi = 0.0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -77,7 +77,7 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation ny = r / sqrt(r^2 + 1) sqr = 1 Va = omega / (ny * sqr) - phi_alv = omega / ny * (nx * (x[1] - 0.5 * r) + ny * (x[2] - 0.5 * r)) - Va * t + phi_alv = omega / ny * (nx * (x[1] - 0.5f0 * r) + ny * (x[2] - 0.5f0 * r)) - Va * t rho = 1.0 v1 = -e * ny * cos(phi_alv) / rho @@ -112,11 +112,11 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations theta = iszero(r) ? 0.0 : acos(z_norm / r) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) * sin(theta) - v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) * sin(theta) - v3 = r > 0.5 ? 0.0 : 0.1882 * cos(theta) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? 1.0 : 1.1691 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) * sin(theta) + v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin(phi) * sin(theta) + v3 = r > 0.5f0 ? 0.0 : 0.1882 * cos(theta) + p = r > 0.5f0 ? 1.0 : 1.245 return prim2cons(SVector(rho, v1, v2, v3, p, 1.0, 1.0, 1.0, 0.0), equations) end @@ -130,9 +130,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one if orientation == 1 f1 = rho_v1 @@ -180,9 +180,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + @@ -342,33 +342,33 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr # for convenience store v⋅B vel_dot_mag_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_mean = 0.5 * rho_avg / beta_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean + 0.5 * mag_norm_avg - B1_avg * B1_avg + f2 = f1 * v1_avg + p_mean + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f3 = f1 * v2_avg - B1_avg * B2_avg f4 = f1 * v3_avg - B1_avg * B3_avg f6 = equations.c_h * psi_avg @@ -376,46 +376,46 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f8 = v1_avg * B3_avg - v3_avg * B1_avg f9 = equations.c_h * B1_avg # total energy flux is complicated and involves the previous eight components - psi_B1_avg = 0.5 * (B1_ll * psi_ll + B1_rr * psi_rr) - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr) + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v1_mag_avg + + 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg - equations.c_h * psi_B1_avg) elseif orientation == 2 f1 = rho_mean * v2_avg f2 = f1 * v1_avg - B2_avg * B1_avg - f3 = f1 * v2_avg + p_mean + 0.5 * mag_norm_avg - B2_avg * B2_avg + f3 = f1 * v2_avg + p_mean + 0.5f0 * mag_norm_avg - B2_avg * B2_avg f4 = f1 * v3_avg - B2_avg * B3_avg f6 = v2_avg * B1_avg - v1_avg * B2_avg f7 = equations.c_h * psi_avg f8 = v2_avg * B3_avg - v3_avg * B2_avg f9 = equations.c_h * B2_avg # total energy flux is complicated and involves the previous eight components - psi_B2_avg = 0.5 * (B2_ll * psi_ll + B2_rr * psi_rr) - v2_mag_avg = 0.5 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr) + v2_mag_avg = 0.5f0 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v2_mag_avg + + 0.5f0 * v2_mag_avg + B2_avg * vel_dot_mag_avg - equations.c_h * psi_B2_avg) else f1 = rho_mean * v3_avg f2 = f1 * v1_avg - B3_avg * B1_avg f3 = f1 * v2_avg - B3_avg * B2_avg - f4 = f1 * v3_avg + p_mean + 0.5 * mag_norm_avg - B3_avg * B3_avg + f4 = f1 * v3_avg + p_mean + 0.5f0 * mag_norm_avg - B3_avg * B3_avg f6 = v3_avg * B1_avg - v1_avg * B3_avg f7 = v3_avg * B2_avg - v2_avg * B3_avg f8 = equations.c_h * psi_avg f9 = equations.c_h * B3_avg # total energy flux is complicated and involves the previous eight components - psi_B3_avg = 0.5 * (B3_ll * psi_ll + B3_rr * psi_rr) - v3_mag_avg = 0.5 * (v3_ll * mag_norm_ll + v3_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B3_avg = 0.5f0 * (B3_ll * psi_ll + B3_rr * psi_rr) + v3_mag_avg = 0.5f0 * (v3_ll * mag_norm_ll + v3_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v3_mag_avg + + 0.5f0 * v3_mag_avg + B3_avg * vel_dot_mag_avg - equations.c_h * psi_B3_avg) end @@ -459,31 +459,31 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = equations.c_h * psi_avg - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) - f9 = equations.c_h * 0.5 * (B1_ll + B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f9 = equations.c_h * 0.5f0 * (B1_ll + B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -494,20 +494,20 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat equations.c_h * (B1_ll * psi_rr + B1_rr * psi_ll))) elseif orientation == 2 f1 = rho_mean * v2_avg - f2 = f1 * v1_avg - 0.5 * (B2_ll * B1_rr + B2_rr * B1_ll) + f2 = f1 * v1_avg - 0.5f0 * (B2_ll * B1_rr + B2_rr * B1_ll) f3 = f1 * v2_avg + p_avg + magnetic_square_avg - - 0.5 * (B2_ll * B2_rr + B2_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B2_ll * B3_rr + B2_rr * B3_ll) + 0.5f0 * (B2_ll * B2_rr + B2_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B2_ll * B3_rr + B2_rr * B3_ll) #f5 below - f6 = 0.5 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) + f6 = 0.5f0 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) f7 = equations.c_h * psi_avg - f8 = 0.5 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) - f9 = equations.c_h * 0.5 * (B2_ll + B2_rr) + f8 = 0.5f0 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) + f9 = equations.c_h * 0.5f0 * (B2_ll + B2_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v2_rr + p_rr * v2_ll + 0.5f0 * (+p_ll * v2_rr + p_rr * v2_ll + (v2_ll * B1_ll * B1_rr + v2_rr * B1_rr * B1_ll) + (v2_ll * B3_ll * B3_rr + v2_rr * B3_rr * B3_ll) - @@ -518,20 +518,20 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat equations.c_h * (B2_ll * psi_rr + B2_rr * psi_ll))) else # orientation == 3 f1 = rho_mean * v3_avg - f2 = f1 * v1_avg - 0.5 * (B3_ll * B1_rr + B3_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B3_ll * B2_rr + B3_rr * B2_ll) + f2 = f1 * v1_avg - 0.5f0 * (B3_ll * B1_rr + B3_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B3_ll * B2_rr + B3_rr * B2_ll) f4 = f1 * v3_avg + p_avg + magnetic_square_avg - - 0.5 * (B3_ll * B3_rr + B3_rr * B3_ll) + 0.5f0 * (B3_ll * B3_rr + B3_rr * B3_ll) #f5 below - f6 = 0.5 * (v3_ll * B1_ll - v1_ll * B3_ll + v3_rr * B1_rr - v1_rr * B3_rr) - f7 = 0.5 * (v3_ll * B2_ll - v2_ll * B3_ll + v3_rr * B2_rr - v2_rr * B3_rr) + f6 = 0.5f0 * (v3_ll * B1_ll - v1_ll * B3_ll + v3_rr * B1_rr - v1_rr * B3_rr) + f7 = 0.5f0 * (v3_ll * B2_ll - v2_ll * B3_ll + v3_rr * B2_rr - v2_rr * B3_rr) f8 = equations.c_h * psi_avg - f9 = equations.c_h * 0.5 * (B3_ll + B3_rr) + f9 = equations.c_h * 0.5f0 * (B3_ll + B3_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v3_rr + p_rr * v3_ll + 0.5f0 * (+p_ll * v3_rr + p_rr * v3_ll + (v3_ll * B1_ll * B1_rr + v3_rr * B1_rr * B1_ll) + (v3_ll * B2_ll * B2_rr + v3_rr * B2_rr * B2_ll) - @@ -568,43 +568,43 @@ end # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = (f1 * v1_avg + (p_avg + magnetic_square_avg) * normal_direction[1] - - 0.5 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) + 0.5f0 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) f3 = (f1 * v2_avg + (p_avg + magnetic_square_avg) * normal_direction[2] - - 0.5 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) + 0.5f0 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) f4 = (f1 * v3_avg + (p_avg + magnetic_square_avg) * normal_direction[3] - - 0.5 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) + 0.5f0 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) #f5 below f6 = (equations.c_h * psi_avg * normal_direction[1] + - 0.5 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + v_dot_n_rr * B1_rr - v1_rr * B_dot_n_rr)) f7 = (equations.c_h * psi_avg * normal_direction[2] + - 0.5 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + v_dot_n_rr * B2_rr - v2_rr * B_dot_n_rr)) f8 = (equations.c_h * psi_avg * normal_direction[3] + - 0.5 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + v_dot_n_rr * B3_rr - v3_rr * B_dot_n_rr)) - f9 = equations.c_h * 0.5 * (B_dot_n_ll + B_dot_n_rr) + f9 = equations.c_h * 0.5f0 * (B_dot_n_ll + B_dot_n_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + 0.5f0 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + (v_dot_n_ll * B1_ll * B1_rr + v_dot_n_rr * B1_rr * B1_ll) + (v_dot_n_ll * B2_ll * B2_rr + v_dot_n_rr * B2_rr * B2_ll) + (v_dot_n_ll * B3_ll * B3_rr + v_dot_n_rr * B3_rr * B3_ll) @@ -951,7 +951,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * (rho_e - - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + B1 * B1 + B2 * B2 + B3 * B3 + psi * psi)) @@ -966,11 +966,12 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -991,8 +992,8 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p * equations.inv_gamma_minus_one + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5 * psi^2 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) + 0.5f0 * psi^2 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -1003,21 +1004,21 @@ end @inline function pressure(u, equations::IdealGlmMhdEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return p end @inline function density_pressure(u, equations::IdealGlmMhdEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return rho * p end @@ -1028,9 +1029,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1038,14 +1039,14 @@ end b3 = B3 / sqrt_rho b_square = b1 * b1 + b2 * b2 + b3 * b3 if orientation == 1 # x-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) elseif orientation == 2 # y-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) else # z-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b3^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b3^2)) end return c_f end @@ -1056,9 +1057,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1072,8 +1073,8 @@ end b2 * normal_direction[2] + b3 * normal_direction[3])^2 / norm_squared - c_f = sqrt((0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * + c_f = sqrt((0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * norm_squared) return c_f end @@ -1096,22 +1097,22 @@ Compute the fast magnetoacoustic wave speed using Roe averages as given by v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) @@ -1133,32 +1134,32 @@ Compute the fast magnetoacoustic wave speed using Roe averages as given by H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum a_square_roe = ((2.0 - equations.gamma) * X + (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) if orientation == 1 # x-direction c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v1_roe elseif orientation == 2 # y-direction c_a_roe = B2_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v2_roe else # z-direction c_a_roe = B3_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v3_roe end @@ -1174,22 +1175,22 @@ end v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) @@ -1211,13 +1212,13 @@ end H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum a_square_roe = ((2.0 - equations.gamma) * X + (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) @@ -1230,7 +1231,7 @@ end c_a_roe = B_roe_dot_n_squared * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) vel_out_roe = (v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + v3_roe * normal_direction[3]) @@ -1270,13 +1271,13 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::IdealGlmMhdEquations3D) - return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + return 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] end # Calculate the magnetic energy for a conservative state `cons'. # OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy @inline function energy_magnetic(cons, ::IdealGlmMhdEquations3D) - return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + return 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) end # Calculate internal energy for a conservative state `cons`