From a9ee0501236dda2a77fe2cd4794530abf6496c2c Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 16 Jun 2024 20:06:14 -1000 Subject: [PATCH 1/9] start --- src/equations/ideal_glm_mhd_1d.jl | 126 +++++++------- src/equations/ideal_glm_mhd_2d.jl | 245 ++++++++++++++------------- src/equations/ideal_glm_mhd_3d.jl | 271 +++++++++++++++--------------- 3 files changed, 323 insertions(+), 319 deletions(-) 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` From 41c90468398f1a404d6c187a14d7ed1aea4fb8d7 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 27 Jun 2024 19:17:15 -1000 Subject: [PATCH 2/9] continue --- src/equations/ideal_glm_mhd_1d.jl | 59 ++++++++++++++-------------- src/equations/ideal_glm_mhd_2d.jl | 64 ++++++++++++++++--------------- 2 files changed, 64 insertions(+), 59 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 2275ddf86c8..3253e4410f8 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -42,13 +42,14 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations1D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = -convert(RealT, 0.2) rho_v3 = -0.5f0 - rho_e = 50.0 - B1 = 3.0 - B2 = -1.2 + rho_e = 50 + B1 = 3 + B2 = -convert(RealT, 1.2) B3 = 0.5f0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) end @@ -61,14 +62,15 @@ An Alfvén wave as smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations1D) # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1], γ = 5/3 - rho = 1.0 - v1 = 0.0 + RealT = eltype(x) + rho = 1 + v1 = 0 # TODO: sincospi - si, co = sincos(2 * pi * x[1]) - v2 = 0.1 * si - v3 = 0.1 * co - p = 0.1 - B1 = 1.0 + si, co = sincos(2 * convert(RealT, pi) * x[1]) + v2 = convert(RealT, 0.1) * si + v3 = convert(RealT, 0.1) * co + p = convert(RealT, 0.1) + B1 = 1 B2 = v2 B3 = v3 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) @@ -86,17 +88,18 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0,) x_norm = x[1] - inicenter[1] r = sqrt(x_norm^2) phi = atan(x_norm) # Calculate primitive variables - 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 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, 0.0, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations) + return prim2cons(SVector(rho, v1, 0, 0, p, 1, 1, 1, 0), equations) end # Calculate 1D flux in for a single point @@ -117,7 +120,7 @@ end f4 = rho_v1 * v3 - B1 * B3 f5 = (kin_en + equations.gamma * p_over_gamma_minus_one + 2 * mag_en) * v1 - B1 * (v1 * B1 + v2 * B2 + v3 * B3) - f6 = 0.0 + f6 = 0 f7 = v1 * B2 - v2 * B1 f8 = v1 * B3 - v3 * B1 @@ -180,7 +183,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, 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 + f6 = 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 @@ -242,7 +245,7 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat 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 + f6 = 0 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 @@ -585,7 +588,7 @@ end b_square = b1^2 + b2^2 + b3^2 c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) return c_f end @@ -628,8 +631,8 @@ as given by # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -649,14 +652,14 @@ as given by 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) * + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * (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) + a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4 * a_square_roe * c_a_roe) c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) return v1_roe, c_f_roe @@ -666,9 +669,9 @@ end @inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations1D) # Pressure p = (equations.gamma - 1) * - (cons[5] - 1 / 2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + (cons[5] - 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] - - 1 / 2 * (cons[6]^2 + cons[7]^2 + cons[8]^2)) + 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2)) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index bec1ee50e4e..73be5c93481 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -48,15 +48,16 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = -convert(RealT, 0.2) rho_v3 = -0.5f0 - rho_e = 50.0 - B1 = 3.0 - B2 = -1.2 + rho_e = 50 + B1 = 3 + B2 = -convert(RealT, 1.2) B3 = 0.5f0 - psi = 0.0 + psi = 0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -70,16 +71,16 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation # domain must be set to [0, 1/cos(α)] x [0, 1/sin(α)], γ = 5/3 alpha = 0.25 * pi x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = 0.1 * sin(2.0 * pi * x_perp) - rho = 1.0 + B_perp = 0.1 * sin(2 * pi * x_perp) + rho = 1 v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = 0.1 * cos(2.0 * pi * x_perp) + v3 = 0.1 * cos(2 * pi * x_perp) p = 0.1 B1 = cos(alpha) + v1 B2 = sin(alpha) + v2 B3 = v3 - psi = 0.0 + psi = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end @@ -95,6 +96,7 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] @@ -102,12 +104,12 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations phi = atan(y_norm, x_norm) # Calculate primitive variables - 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 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, v2, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations) + return prim2cons(SVector(rho, v1, v2, 0, p, 1, 1, 1, 0), equations) end # Pre-defined source terms should be implemented as @@ -1131,8 +1133,8 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 - return (equations.gamma - 1.0) * - SVector(0.5f0 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) + return (equations.gamma - 1) * + SVector(0.5f0 * v_square, -v1, -v2, -v3, 1, -B1, -B2, -B3, -psi) end @inline function density_pressure(u, equations::IdealGlmMhdEquations2D) @@ -1163,10 +1165,10 @@ end b_square = b1 * b1 + b2 * b2 + b3 * b3 if orientation == 1 # x-direction c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) else c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)) end return c_f end @@ -1236,8 +1238,8 @@ as given by # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1257,21 +1259,21 @@ as given by 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) * + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * (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) + 4 * a_square_roe * c_a_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) + 4 * a_square_roe * c_a_roe) c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v2_roe end @@ -1308,8 +1310,8 @@ end # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1329,8 +1331,8 @@ end 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) * + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed @@ -1404,7 +1406,7 @@ end # State validation for Newton-bisection method of subcell IDP limiting @inline function Base.isvalid(u, equations::IdealGlmMhdEquations2D) p = pressure(u, equations) - if u[1] <= 0.0 || p <= 0.0 + if u[1] <= 0 || p <= 0 return false end return true From 14bfdaf083b6a2267d1281a4696956c0ea265508 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 28 Jun 2024 07:36:30 -1000 Subject: [PATCH 3/9] complete equations --- src/equations/ideal_glm_mhd_2d.jl | 15 ++++--- src/equations/ideal_glm_mhd_3d.jl | 73 ++++++++++++++++--------------- 2 files changed, 46 insertions(+), 42 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 73be5c93481..ab2a4b066a1 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -69,14 +69,15 @@ An Alfvén wave as smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations2D) # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1/cos(α)] x [0, 1/sin(α)], γ = 5/3 - alpha = 0.25 * pi + RealT = eltype(x) + alpha = 0.25f0 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = 0.1 * sin(2 * pi * x_perp) + B_perp = convert(RealT, 0.1) * sinpi(2 * x_perp) rho = 1 v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = 0.1 * cos(2 * pi * x_perp) - p = 0.1 + v3 = convert(RealT, 0.1) * cospi(2 * x_perp) + p = convert(RealT, 0.1) B1 = cos(alpha) + v1 B2 = sin(alpha) + v2 B3 = v3 @@ -1355,11 +1356,11 @@ end @inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations2D) # Pressure p = (equations.gamma - 1) * - (cons[5] - 1 / 2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + (cons[5] - 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] - - 1 / 2 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) - - 1 / 2 * cons[9]^2) + 0.5f0 * cons[9]^2) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 679afa62b24..fc43b81759e 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -46,15 +46,16 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = -convert(RealT, 0.2) rho_v3 = -0.5f0 - rho_e = 50.0 - B1 = 3.0 - B2 = -1.2 + rho_e = 50 + B1 = 3 + B2 = -convert(RealT, 1.2) B3 = 0.5f0 - psi = 0.0 + psi = 0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -67,19 +68,20 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation # Alfvén wave in three space dimensions # Altmann thesis http://dx.doi.org/10.18419/opus-3895 # domain must be set to [-1, 1]^3, γ = 5/3 + RealT = eltype(x) p = 1 - omega = 2 * pi # may be multiplied by frequency + omega = 2 * convert(RealT, pi) # may be multiplied by frequency # r: length-variable = length of computational domain - r = 2 + r = 2.0f0 # e: epsilon = 0.2 - e = 0.2 + e = convert(RealT, 0.2) nx = 1 / sqrt(r^2 + 1) ny = r / sqrt(r^2 + 1) sqr = 1 Va = omega / (ny * sqr) phi_alv = omega / ny * (nx * (x[1] - 0.5f0 * r) + ny * (x[2] - 0.5f0 * r)) - Va * t - rho = 1.0 + rho = 1 v1 = -e * ny * cos(phi_alv) / rho v2 = e * nx * cos(phi_alv) / rho v3 = e * sin(phi_alv) / rho @@ -103,22 +105,23 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0, 0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) phi = atan(y_norm, x_norm) - theta = iszero(r) ? 0.0 : acos(z_norm / r) + theta = iszero(r) ? zero(RealT) : acos(z_norm / r) # Calculate primitive variables - 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 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) * sin(theta) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) * sin(theta) + v3 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(theta) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, v2, v3, p, 1.0, 1.0, 1.0, 0.0), equations) + return prim2cons(SVector(rho, v1, v2, v3, p, 1, 1, 1, 0), equations) end # Pre-defined source terms should be implemented as @@ -1040,13 +1043,13 @@ end b_square = b1 * b1 + b2 * b2 + b3 * b3 if orientation == 1 # x-direction c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) elseif orientation == 2 # y-direction c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)) else # z-direction c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b3^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b3^2)) end return c_f end @@ -1117,8 +1120,8 @@ Compute the fast magnetoacoustic wave speed using Roe averages as given by # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1138,27 +1141,27 @@ Compute the fast magnetoacoustic wave speed using Roe averages as given by 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) * + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * (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) + 4 * a_square_roe * c_a_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) + 4 * a_square_roe * c_a_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) + 4 * a_square_roe * c_a_roe) c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v3_roe end @@ -1195,8 +1198,8 @@ end # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1216,8 +1219,8 @@ end 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) * + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed @@ -1243,11 +1246,11 @@ end @inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations3D) # Pressure p = (equations.gamma - 1) * - (cons[5] - 1 / 2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + (cons[5] - 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] - - 1 / 2 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) - - 1 / 2 * cons[9]^2) + 0.5f0 * cons[9]^2) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) From b90b4d4f2932f55f315ddbfb798c50b579760c83 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 28 Jun 2024 11:31:00 -1000 Subject: [PATCH 4/9] unit test 1D --- test/test_type.jl | 88 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 79 insertions(+), 9 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index c65d9b62679..b54b59b92dd 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -109,10 +109,9 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT)) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT)) orientation = 1 directions = [1, 2] - cons = SVector(one(RealT), one(RealT), one(RealT)) surface_flux_function = flux_lax_friedrichs @@ -204,12 +203,11 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT), zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), - one(RealT)) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) orientations = [1, 2] directions = [1, 2, 3, 4] normal_direction = SVector(one(RealT), zero(RealT)) - cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) surface_flux_function = flux_lax_friedrichs flux_lmars = FluxLMARS(340) @@ -380,13 +378,11 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT), zero(RealT), zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), - one(RealT), one(RealT)) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) orientations = [1, 2, 3] directions = [1, 2, 3, 4, 5, 6] normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) - cons = SVector(one(RealT), one(RealT), one(RealT), - one(RealT), one(RealT)) surface_flux_function = flux_lax_friedrichs flux_lmars = FluxLMARS(340) @@ -1028,6 +1024,80 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Ideal Glm Mhd 1D" begin + for RealT in (Float32, Float64) + equations = @inferred IdealGlmMhdEquations1D(RealT(1.4)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT)) + orientation = 1 + directions = [1, 2] + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == + RealT + @test typeof(@inferred Trixi.max_abs_speeds(u, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == + RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + + for direction in directions + @test typeof(Trixi.calc_fast_wavespeed(cons, direction, equations)) == RealT + @test eltype(Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, direction, + equations)) == RealT + end + + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_magnetic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred cross_helicity(cons, equations)) == RealT + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) From df560f1c824cfa7a57b9cba17f0c19e4b9afa7cd Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sat, 29 Jun 2024 19:09:19 -1000 Subject: [PATCH 5/9] unit test 2D --- test/test_type.jl | 149 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 144 insertions(+), 5 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index b54b59b92dd..a861c0266b0 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1030,11 +1030,12 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT), - one(RealT), - one(RealT), - one(RealT), - one(RealT)) + u = u_ll = u_rr = cons = SVector(one(RealT), zero(RealT), zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT)) orientation = 1 directions = [1, 2] @@ -1098,6 +1099,144 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Ideal Glm Mhd 2D" begin + for RealT in (Float32, Float64) + equations = @inferred IdealGlmMhdEquations2D(RealT(1.4)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), zero(RealT), zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), + zero(RealT)) + nonconservative_type_local = Trixi.NonConservativeLocal() + nonconservative_type_symmetric = Trixi.NonConservativeSymmetric() + nonconservative_terms = [1, 2] + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + normal_direction, + equations)) == RealT + else + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, + u_rr, + orientation, + equations)) == + RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == + RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, orientation, + equations)) == RealT + end + for nonconservative_term in nonconservative_terms + @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, + orientation, + equations, + nonconservative_type_local, + nonconservative_term)) == + RealT + @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, + u_rr, + orientation, + equations, + nonconservative_type_symmetric, + nonconservative_term)) == + RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test eltype(@inferred Trixi.gradient_conservative(pressure, u, equations)) == + RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, normal_direction, + equations)) == RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + normal_direction, + equations)) == RealT + for orientation in orientations + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, orientation, + equations)) == + RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + orientation, + equations)) == RealT + end + + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_magnetic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred cross_helicity(cons, equations)) == RealT + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) From 28f45cc0db8580804ad5afadd00bd5d7a8ecd2a5 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sat, 29 Jun 2024 19:29:28 -1000 Subject: [PATCH 6/9] unit test 3D --- test/test_type.jl | 114 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index a861c0266b0..fc9a9561e6c 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1237,6 +1237,120 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Ideal Glm Mhd 3D" begin + for RealT in (Float32, Float64) + equations = @inferred IdealGlmMhdEquations3D(RealT(1.4)) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), zero(RealT), zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT)) + orientations = [1, 2, 3] + directions = [1, 2, 3, 4, 5, 6] + normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), + zero(RealT), + zero(RealT)) + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + normal_direction, + equations)) == RealT + else + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, + equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == + RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, normal_direction, + equations)) == RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + normal_direction, + equations)) == RealT + for orientation in orientations + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, orientation, + equations)) == + RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + orientation, + equations)) == RealT + end + + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_magnetic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred cross_helicity(cons, equations)) == RealT + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) From cbbea3ca6f8f3008843eba8b0521adf90053318f Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 1 Jul 2024 08:09:43 -1000 Subject: [PATCH 7/9] check errors --- src/equations/ideal_glm_mhd_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index ab2a4b066a1..01a825bbd9c 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -72,11 +72,11 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation RealT = eltype(x) alpha = 0.25f0 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = convert(RealT, 0.1) * sinpi(2 * x_perp) + B_perp = convert(RealT, 0.1) * sin(2 * convert(RealT, pi) * x_perp) rho = 1 v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = convert(RealT, 0.1) * cospi(2 * x_perp) + v3 = convert(RealT, 0.1) * cos(2 * convert(RealT, pi) * x_perp) p = convert(RealT, 0.1) B1 = cos(alpha) + v1 B2 = sin(alpha) + v2 From 5b29c15d94231c93248baa882ec4de42450a8d0b Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 1 Jul 2024 13:48:06 -1000 Subject: [PATCH 8/9] check again --- src/equations/ideal_glm_mhd_3d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index fc43b81759e..2ffaa575243 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -72,7 +72,7 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation p = 1 omega = 2 * convert(RealT, pi) # may be multiplied by frequency # r: length-variable = length of computational domain - r = 2.0f0 + r = convert(RealT, 2) # e: epsilon = 0.2 e = convert(RealT, 0.2) nx = 1 / sqrt(r^2 + 1) From 2cfe92bcc95ca12f5233a60a4e1b73ad360fddb0 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 1 Jul 2024 14:25:10 -1000 Subject: [PATCH 9/9] fixed --- src/equations/ideal_glm_mhd_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 01a825bbd9c..ab2a4b066a1 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -72,11 +72,11 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation RealT = eltype(x) alpha = 0.25f0 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = convert(RealT, 0.1) * sin(2 * convert(RealT, pi) * x_perp) + B_perp = convert(RealT, 0.1) * sinpi(2 * x_perp) rho = 1 v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = convert(RealT, 0.1) * cos(2 * convert(RealT, pi) * x_perp) + v3 = convert(RealT, 0.1) * cospi(2 * x_perp) p = convert(RealT, 0.1) B1 = cos(alpha) + v1 B2 = sin(alpha) + v2