diff --git a/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl new file mode 100644 index 00000000000..80a42697d51 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl @@ -0,0 +1,613 @@ + +using OrdinaryDiffEq +using Trixi +using NLsolve: nlsolve +using Infiltrator + +# Warm moist bubble test case from +# Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) + +@inline function getqvs(eps, p, t) + # Equation for saturation vapor pressure + # Bolton (1980, MWR, p. 1047) + es = 611.2 * exp(17.67 * (t - 273.15) / (t - 29.65)) + return eps * es / (p - es) +end + +# Compute a base state in hydrostatic equilibrium +# The hydrostatic ODE is integrated upwards from surface (z=0) to `height` +# A one dimensional distribution (in z) with `nz` points is stored +# Adapted from https://www2.mmm.ucar.edu/people/bryan/Code/mbm.F +function calculate_hydrostatic_base_state(g, c_pd, R_d, c_pv, R_v, c_pl, L_00, + p_ref, theta_ref, r_t, nz, height) + eps = R_d / R_v + tolerance = 0.0001 # tolerance for iterative scheme + relaxation = 0.3 # relaxation factor for iterative scheme + + theta_0 = Vector{Float64}(undef, nz) # potential temperature + theta_rho_0 = Vector{Float64}(undef, nz) # density potential temperature + p_0 = Vector{Float64}(undef, nz) # pressure + q_v0 = Vector{Float64}(undef, nz) # vapor mixing ratio + q_l0 = Vector{Float64}(undef, nz) # liquid water mixing ratio + + # Start at the surface + p_0[1] = p_ref # ==> exner = 1 + q_v0[1] = getqvs(eps, p_ref, theta_ref) # T = theta_ref * exner = theta_ref + q_l0[1] = r_t - q_v0[1] + theta_0[1] = theta_ref + theta_rho_0[1] = theta_ref * (1.0 + inv(eps) * q_v0[1]) / (1.0 + q_v0[1] + q_l0[1]) + + T_last = theta_0[1] # T = theta_ref * exner = theta_ref + exner_last = 1 + + # Integrate upwards + dz = height / (nz - 1) + for k in 2:nz + theta_new = theta_0[k - 1] + theta_tmp = theta_0[k - 1] + theta_rho_tmp = theta_rho_0[k - 1] + + exner_tmp = 1.0 + p_tmp = 1.0 + T_tmp = 1.0 + q_v_tmp = 1.0 + q_l_tmp = 1.0 + steps = 0 + converged = false + while !converged + steps = steps + 1 + + # Trapezoidal rule? Only applied to theta_rho? + exner_tmp = exner_last - + dz * g / (c_pd * 0.5 * (theta_rho_0[k - 1] + theta_rho_tmp)) + p_tmp = p_ref * (exner_tmp^(c_pd / R_d)) + T_tmp = theta_new * exner_tmp + q_v_tmp = getqvs(eps, p_tmp, T_tmp) + q_l_tmp = r_t - q_v_tmp + theta_rho_tmp = theta_tmp * (1.0 + inv(eps) * q_v_tmp) / + (1.0 + q_v_tmp + q_l_tmp) + + # Trapezoidal rule? + T_bar = 0.5 * (T_last + T_tmp) + q_vbar = 0.5 * (q_v0[k - 1] + q_v_tmp) + q_lbar = 0.5 * (q_l0[k - 1] + q_l_tmp) + + lhv = L_00 - (c_pl - c_pv) * T_bar # latent heat + R_m = R_d + R_v * q_vbar # gas constant for moist air + c_pm = c_pd + c_pv * q_vbar + c_pl * q_lbar # c_p for moist air + + # An exact thermodynamic equation (e.g., for CM1) + theta_tmp = theta_0[k - 1] * + exp(-lhv * (q_v_tmp - q_v0[k - 1]) / (c_pm * T_bar) + + (R_m / c_pm - R_d / c_pd) * log(p_tmp / p_0[k - 1])) + + if abs(theta_tmp - theta_new) > tolerance + theta_new = theta_new + relaxation * (theta_tmp - theta_new) + if steps > 100 + error("calculate_hydrostatic_base_state: not converging!") + end + else + converged = true + if q_l_tmp < 0.0 + error("calculate_hydrostatic_base_state: liquid ratio negative!") + end + end + end + + theta_0[k] = theta_tmp + theta_rho_0[k] = theta_tmp * (1.0 + inv(eps) * q_v_tmp) / (1.0 + q_v_tmp + q_l_tmp) + p_0[k] = p_tmp + q_v0[k] = q_v_tmp + q_l0[k] = q_l_tmp + + T_last = T_tmp + exner_last = exner_tmp + end + + return theta_0, theta_rho_0, p_0, q_v0, q_l0 +end + +struct WarmBubbleSetup + g::Float64 # gravity of earth + c_pd::Float64 # heat capacity for constant pressure (dry air) + c_vd::Float64 # heat capacity for constant volume (dry air) + R_d::Float64 # gas constant (dry air) + c_pv::Float64 # heat capacity for constant pressure (water vapor) + c_vv::Float64 # heat capacity for constant volume (water vapor) + R_v::Float64 # gas constant (water vapor) + c_pl::Float64 # heat capacity for constant pressure (liquid water) + p_ref::Float64 # reference pressure + theta_ref::Float64 # potential temperature at surface + L_00::Float64 # latent heat of evaporation at 0 K + r_t::Float64 # total water mixing ratio + theta_0::Vector{Float64} # potential temperature + theta_rho_0::Vector{Float64} # density potential temperature + p_0::Vector{Float64} # pressure + q_v0::Vector{Float64} # vapor mixing ratio + q_l0::Vector{Float64} # liquid water mixing ratio + nz::Int64 # resolution for vertical hydrostatic distribution + height::Float64 # maximal height in vertical hydrostatic distribution + + function WarmBubbleSetup(; g = 9.81, + c_pd = 1004.0, c_vd = 717.0, R_d = c_pd - c_vd, + c_pv = 1885, c_vv = 1424.0, R_v = c_pv - c_vv, + c_pl = 4186.0, + p_ref = 100_000.0, theta_ref = 289.8486, + L_00 = 2.5e6 + (c_pl - c_pv) * 273.15, + r_t = 0.020, + nz, height) + theta_0, + theta_rho_0, + p_0, + q_v0, + q_l0 = calculate_hydrostatic_base_state(g, c_pd, R_d, c_pv, R_v, c_pl, L_00, p_ref, + theta_ref, r_t, nz, height) + new(g, c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_ref, theta_ref, L_00, r_t, theta_0, + theta_rho_0, p_0, q_v0, q_l0, nz, height) + end +end + +function get_c_p(setup::WarmBubbleSetup) + return (setup.c_pd, setup.c_pv, setup.c_pl) +end + +function get_c_v(setup::WarmBubbleSetup) + return (setup.c_vd, setup.c_vv, setup.c_pl) # c_pl = c_vl for liquid phase +end + +############################################################################### +# source terms +@inline function source_terms_geopotential(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) + @unpack g = setup + _, rho_v2 = u + rho = density(u, equations) + + return SVector(zero(eltype(u)), -g * rho, -g * rho_v2, + zero(eltype(u)), zero(eltype(u)), zero(eltype(u))) +end + +# condensation, equation by Rutledge and Hobbs (1983) +@inline function source_terms_phase_change_rh(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, R_d, c_pv, R_v, c_pl, L_00 = setup + eps = R_d / R_v + rho_v1, rho_v2, rho_e, rho_qd, rho_qv, rho_ql = u + rho = density(u, equations) + T = temperature(u, equations) + p = pressure(u, equations) + + q_vs = getqvs(eps, p, T) + q_v = rho_qv / rho + lhv = L_00 - (c_pl - c_pv) * T + rcond = (q_v - q_vs) / (1.0 + q_vs * lhv^2 / (c_pd * R_v * T^2)) + #qcond(mgs) = Min( Max( 0.0, tmp ), (qvap(mgs)-qvs(mgs)) ) + return SVector(zero(eltype(u)), zero(eltype(u)), zero(eltype(u)), + zero(eltype(u)), -rcond, rcond) +end + +@inline function source_terms_phase_change_fb(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) + + # This source term models the phase chance between could water and vapor. + @unpack R_v = setup + rho_v1, rho_v2, rho_e, rho_qd, rho_qv, rho_ql = u + rho = density(u, equations) + T = temperature(u, equations) + + T_C = T - 273.15 + # saturation vapor pressure + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + + # saturation density of vapor + rho_star_qv = p_vs / (R_v * T) + + # Fisher-Burgmeister-Function + a = rho_star_qv - rho_qv + b = rho - rho_qv - rho_qd + + # saturation control factor + # < 1: stronger saturation effect + # > 1: weaker saturation effect + C = 94.0 + + Q_ph = (a + b - sqrt(a^2 + b^2)) * C + + return SVector(zero(eltype(u)), zero(eltype(u)), zero(eltype(u)), + zero(eltype(u)), Q_ph, -Q_ph) +end + +@inline function (setup::WarmBubbleSetup)(u, x, t, + equations::CompressibleMoistEulerEquations2D) + return source_terms_geopotential(u, setup, equations) + + source_terms_phase_change_fb(u, setup, equations) +end + +############################################################################### +# Initial condition + +# Use the precomputed vertical distribution in `WarmBubbleSetup` to interpolate and perturb +# the initial solution +@inline function (setup::WarmBubbleSetup)(x, t, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_vd, R_d, c_vv, R_v, c_pl, L_00, r_t, p_ref, theta_0, + theta_rho_0, p_0, q_v0, q_l0, nz, height = setup + eps = R_d / R_v + tolerance = 0.0001 # tolerance for iterative scheme + + # get value at current z via interpolation + dz = height / (nz - 1) + k_l = convert(Int, floor(x[2] / dz)) + 1 + z_l = (k_l - 1) * dz + z_u = k_l * dz + if k_l == nz + k_u = nz + else + k_u = k_l + 1 + end + theta, theta_rho, p, q_v, q_l = map((a, b) -> (z_u - x[2]) * a / dz + + (x[2] - z_l) * b / dz, + [ + theta_0[k_l], + theta_rho_0[k_l], + p_0[k_l], + q_v0[k_l], + q_l0[k_l], + ], + [ + theta_0[k_u], + theta_rho_0[k_u], + p_0[k_u], + q_v0[k_u], + q_l0[k_u], + ]) + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + exner = (p / p_ref)^(R_d / c_pd) + T = theta * exner + + if r <= radius + theta_perturbation = 2 * cospi(0.5 * r / radius)^2 + theta_tmp = theta + converged = false + steps = 0 + while !converged + steps = steps + 1 + theta = ((theta_perturbation / 300.0) + (1.0 + r_t) / (1.0 + q_v)) * + theta_rho * (1.0 + q_v) / (1.0 + inv(eps) * q_v) + T = theta * exner + q_v = getqvs(eps, p, T) + q_l = r_t - q_v + if abs(theta - theta_tmp) < tolerance + converged = true + if q_l < 0.0 + error("calculate_hydrostatic_base_state: liquid ratio negative!") + end + else + theta_tmp = theta + if steps > 100 + error("initial condition perturbation: not converging!") + end + end + end + end + + rho_d = p / (R_d * T * (1.0 + q_v * inv(eps))) # p = rho_d * R_d * T * (1 + q_v / eps) + rho_v = rho_d * q_v + rho_l = rho_d * q_l + rho = rho_d + rho_l + rho_v + v1 = 0.0 + v2 = 0.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_e = (c_vd * rho_d + c_vv * rho_v + c_pl * rho_l) * T + + L_00 * rho_v + + 0.5 * rho * (v1^2 + v2^2) + return SVector(rho_v1, rho_v2, rho_e, rho_d, rho_v, rho_l) +end + +# Initial condition Knoth +struct AtmossphereLayers{RealT <: Real} + setup::WarmBubbleSetup + # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql + LayerData::Matrix{RealT} + total_height::RealT + preciseness::Int + layers::Int + ground_state::NTuple{2, RealT} + equivalentpotential_temperature::RealT + mixing_ratios::NTuple{2, RealT} +end + +function AtmossphereLayers(setup; total_height = 10010.0, preciseness = 10, + ground_state = (1.4, 100000.0), + equivalentpotential_temperature = 320, + mixing_ratios = (0.02, 0.02), RealT = Float64) + @unpack c_pd, R_d, c_pv, R_v, c_pl = setup + rho0, p0 = ground_state + r_t0, r_v0 = mixing_ratios + theta_e0 = equivalentpotential_temperature + + rho_qv0 = rho0 * r_v0 + T0 = theta_e0 + y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] + + n = convert(Int, total_height / preciseness) + dz = 0.01 + LayerData = zeros(RealT, n + 1, 4) + + F = generate_function_of_y(dz, y0, r_t0, theta_e0, setup) + sol = nlsolve(F, y0) + p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero + + rho_d = rho / (1 + r_t) + rho_ql = rho - rho_d - rho_qv + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + LayerData[1, :] = [rho, rho_theta, rho_qv, rho_ql] + for i in (1:n) + #println("Layer", i) + y0 = deepcopy(sol.zero) + dz = preciseness + F = generate_function_of_y(dz, y0, r_t0, theta_e0, setup) + sol = nlsolve(F, y0) + p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero + rho_d = rho / (1 + r_t) + rho_ql = rho - rho_d - rho_qv + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + LayerData[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] + end + + return AtmossphereLayers{RealT}(setup, LayerData, total_height, dz, n, ground_state, + theta_e0, mixing_ratios) +end + +function moist_state(y, dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) + @unpack g, c_pd, R_d, c_pv, R_v, c_pl, p_ref, theta_ref, L_00 = setup + (p, rho, T, r_t, r_v, rho_qv, theta_e) = y + p0 = y0[1] + + F = zeros(7, 1) + rho_d = rho / (1 + r_t) + p_d = R_d * rho_d * T + T_C = T - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + L = L_00 - (c_pl - c_pv) * T + + F[1] = (p - p0) / dz + g * rho + F[2] = p - (R_d * rho_d + R_v * rho_qv) * T + # H = 1 is assumed + F[3] = (theta_e - + T * (p_d / p_ref)^(-R_d / (c_pd + c_pl * r_t)) * + exp(L * r_v / ((c_pd + c_pl * r_t) * T))) + F[4] = r_t - r_t0 + F[5] = rho_qv - rho_d * r_v + F[6] = theta_e - theta_e0 + a = p_vs / (R_v * T) - rho_qv + b = rho - rho_qv - rho_d + # H=1 => phi=0 + F[7] = a + b - sqrt(a * a + b * b) + + return F +end + +function generate_function_of_y(dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) + function function_of_y(y) + return moist_state(y, dz, y0, r_t0, theta_e0, setup) + end +end + +function initial_condition_moist_bubble(x, t, setup::WarmBubbleSetup, + AtmosphereLayers::AtmossphereLayers) + @unpack LayerData, preciseness, total_height = AtmosphereLayers + dz = preciseness + z = x[2] + if (z > total_height && !(isapprox(z, total_height))) + error("The atmossphere does not match the simulation domain") + end + n = convert(Int, floor((z + eps()) / dz)) + 1 + z_l = (n - 1) * dz + (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] + z_r = n * dz + if (z_l == total_height) + z_r = z_l + dz + n = n - 1 + end + (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n + 1, :] + rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz + rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / + dz + rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz + rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz + + rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, + setup) + + v1 = 0.0 + v2 = 0.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) + + rho_qd = rho - rho_qv - rho_ql + + return SVector(rho_v1, rho_v2, rho_E, rho_qd, rho_qv, rho_ql) +end + +function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbleSetup) + @unpack c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_ref, L_00 = setup + xc = 10000.0 + zc = 2000.0 + rc = 2000.0 + Δθ = 2.0 + + r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + p_loc = p_ref * (R_d * rho_theta / p_ref)^(1 / (1 - kappa_M)) + T_loc = p_loc / (R_d * rho_d + R_v * rho_qv) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + + p_v = rho_qv * R_v * T_loc + p_d = p_loc - p_v + T_C = T_loc - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + H = p_v / p_vs + r_v = rho_qv / rho_d + r_l = rho_ql / rho_d + r_t = r_v + r_l + + # Aequivalentpotential temperature + a = T_loc * (p_ref / p_d)^(R_d / (c_pd + r_t * c_pl)) + b = H^(-r_v * R_v / c_pd) + L_v = L_00 + (c_pv - c_pl) * T_loc + c = exp(L_v * r_v / ((c_pd + r_t * c_pl) * T_loc)) + aeq_pot = (a * b * c) + + # Assume pressure stays constant + if (r < rc && Δθ > 0) + kappa = 1 - c_vd / c_pd + # Calculate background density potential temperature + θ_dens = rho_theta / rho * (p_loc / p_ref)^(kappa_M - kappa) + # Calculate perturbed density potential temperature + θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300) + rt = (rho_qv + rho_ql) / rho_d + rv = rho_qv / rho_d + # Calculate moist potential temperature + θ_loc = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rv) + # Adjust varuables until the temperature is met + if rt > 0 + while true + T_loc = θ_loc * (p_loc / p_ref)^kappa + T_C = T_loc - 273.15 + # SaturVapor + pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + rho_d_new = (p_loc - pvs) / (R_d * T_loc) + rvs = pvs / (R_v * rho_d_new * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + if abs(θ_new - θ_loc) <= θ_loc * 1.0e-12 + break + else + θ_loc = θ_new + end + end + else + rvs = 0 + T_loc = θ_loc * (p_loc / p_ref)^kappa + rho_d_new = p_loc / (R_d * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + end + rho_qv = rvs * rho_d_new + rho_ql = (rt - rvs) * rho_d_new + rho = rho_d_new * (1 + rt) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * θ_dens_new * (p_loc / p_ref)^(kappa - kappa_M) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + end + return SVector(rho, rho_e, rho_qv, rho_ql) +end + +############################################################################### +# semidiscretization of the compressible moist Euler equations +nelements_z = 100 +polydeg = 1 +height = 10_000.0 + +warm_bubble_setup = WarmBubbleSetup(; nz = nelements_z * (polydeg + 1), height = height) + +AtmossphereData = AtmossphereLayers(warm_bubble_setup) + +# Create the initial condition with the initial data set +function initial_condition_moist(x, t, equations) + return initial_condition_moist_bubble(x, t, warm_bubble_setup, AtmossphereData) +end + +equations = CompressibleMoistEulerEquations2D(; c_p = get_c_p(warm_bubble_setup), + c_v = get_c_v(warm_bubble_setup), + L_00 = warm_bubble_setup.L_00) + +boundary_condition = (x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +basis = LobattoLegendreBasis(polydeg) + +#surface_flux = FluxLMARS(340.0) +surface_flux = flux_lax_friedrichs +#volume_flux = flux_chandrashekar +#volume_flux = flux_ranocha +volume_flux = flux_kennedy_gruber +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, height) + +cells_per_dimension = (2 * nelements_z, nelements_z) + +# Create curved mesh with 64 x 32 elements +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + boundary_conditions = boundary_condition, + source_terms = warm_bubble_setup) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1200.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +#solution_variables = cons2aeqpot +solution_variables = cons2prim + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +#analysis_callback = AnalysisCallback(semi, interval=analysis_interval, extra_analysis_errors=(:entropy_conservation_error, ), extra_analysis_integrals=(entropy, energy_total, Trixi.saturation_pressure)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out_bf_struct200_p1_cfl02_lf_kennedy_gruber", + solution_variables = solution_variables) + +stepsize_callback = StepsizeCallback(cfl = 0.2) + +callbacks = CallbackSet(summary_callback, + #analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, + maxiters = 1e6, + callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index da7359999c5..a16e6ae59db 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -150,6 +150,7 @@ export AcousticPerturbationEquations2D, CompressibleEulerEquations3D, CompressibleEulerMulticomponentEquations1D, CompressibleEulerMulticomponentEquations2D, + CompressibleMoistEulerEquations2D, CompressibleEulerEquationsQuasi1D, IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D, @@ -219,7 +220,7 @@ export initial_condition_eoc_test_coupled_euler_gravity, export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, cons2entropy, entropy2cons export density, pressure, density_pressure, velocity, global_mean_vars, - equilibrium_distribution, waterheight_pressure + equilibrium_distribution, waterheight_pressure, temperature export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, enstrophy diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 60fce222f21..d02ceccf56b 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -102,195 +102,167 @@ end RealT end -function varnames(::typeof(cons2cons), - equations::CompressibleEulerMulticomponentEquations2D) - cons = ("rho_v1", "rho_v2", "rho_e") - rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) - return (cons..., rhos...) -end - -function varnames(::typeof(cons2prim), - equations::CompressibleEulerMulticomponentEquations2D) - prim = ("v1", "v2", "p") - rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) - return (prim..., rhos...) -end - -# Set initial conditions at physical location `x` for time `t` - """ - initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D) + temperature(u, equations::CompressibleEulerMulticomponentEquations2D) -A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +Calculate temperature. """ -function initial_condition_convergence_test(x, t, - equations::CompressibleEulerMulticomponentEquations2D) - c = 2 - A = 0.1 - L = 2 - f = 1 / L - omega = 2 * pi * f - ini = c + A * sin(omega * (x[1] + x[2] - t)) - - v1 = 1.0 - v2 = 1.0 - - rho = ini - - # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1) - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho - for i in eachcomponent(equations)) - - prim1 = rho * v1 - prim2 = rho * v2 - prim3 = rho^2 - - prim_other = SVector{3, real(equations)}(prim1, prim2, prim3) +@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) + rho_v1, rho_v2, rho_e = u + @unpack cv = equations - return vcat(prim_other, prim_rho) -end + rho = density(u, equations) + help1 = zero(rho) -""" - source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D) + # compute weighted average of cv + # normalization by rho not required, cancels below + for i in eachcomponent(equations) + help1 += u[i + 3] * cv[i] + end -Source terms used for convergence tests in combination with -[`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). -""" -@inline function source_terms_convergence_test(u, x, t, - equations::CompressibleEulerMulticomponentEquations2D) - # Same settings as in `initial_condition` - c = 2 - A = 0.1 - L = 2 - f = 1 / L - omega = 2 * pi * f - - gamma = totalgamma(u, equations) - - x1, x2 = x - si, co = sincos((x1 + x2 - t) * omega) - tmp1 = co * A * omega - tmp2 = si * A - tmp3 = gamma - 1 - tmp4 = (2 * c - 1) * tmp3 - tmp5 = (2 * tmp2 * gamma - 2 * tmp2 + tmp4 + 1) * tmp1 - tmp6 = tmp2 + c - - # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 - du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - tmp1 - for i in eachcomponent(equations)) - - du1 = tmp5 - du2 = tmp5 - du3 = 2 * ((tmp6 - 1.0) * tmp3 + tmp6 * gamma) * tmp1 - - du_other = SVector{3, real(equations)}(du1, du2, du3) - - return vcat(du_other, du_rho) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + T = (rho_e - 0.5 * rho * v_square) / help1 + return T end """ - initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D) + totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) -A for multicomponent adapted weak blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +Function that calculates the total gamma out of all partial gammas using the +partial density fractions as well as the partial specific heats at constant volume. """ -function initial_condition_weak_blast_wave(x, t, - equations::CompressibleEulerMulticomponentEquations2D) - # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) - # Set up polar coordinates - inicenter = SVector(0.0, 0.0) - x_norm = x[1] - inicenter[1] - y_norm = x[2] - inicenter[2] - r = sqrt(x_norm^2 + y_norm^2) - phi = atan(y_norm, x_norm) - sin_phi, cos_phi = sincos(phi) - - prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - 1.0 : - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - 1.1691 - for i in eachcomponent(equations)) +@inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, cp = equations - 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 + help1 = zero(u[1]) + help2 = zero(u[1]) - prim_other = SVector{3, real(equations)}(v1, v2, p) + # compute weighted averages of cp and cv + # normalization by total rho not required, would cancel below + for i in eachcomponent(equations) + help1 += u[i + 3] * cp[i] + help2 += u[i + 3] * cv[i] + end - return prim2cons(vcat(prim_other, prim_rho), equations) + return help1 / help2 end -# Calculate 1D flux for a single point -@inline function flux(u, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) +""" +cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + +Convert conservative variables to entropy. +""" +@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations rho_v1, rho_v2, rho_e = u rho = density(u, equations) + # Multicomponent stuff + help1 = zero(rho) + gas_constant = zero(rho) + for i in eachcomponent(equations) + help1 += u[i + 3] * cv[i] + gas_constant += gas_constants[i] * (u[i + 3] / rho) + end + v1 = rho_v1 / rho v2 = rho_v2 / rho - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + v_square = v1^2 + v2^2 + p = pressure(u, equations) + rho_p = rho / p + T = temperature(u, equations) - if orientation == 1 - f_rho = densities(u, v1, equations) - f1 = rho_v1 * v1 + p - f2 = rho_v2 * v1 - f3 = (rho_e + p) * v1 - else - f_rho = densities(u, v2, equations) - f1 = rho_v1 * v2 - f2 = rho_v2 * v2 + p - f3 = (rho_e + p) * v2 - end + entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 3])) - + v_square / (2 * T)) + for i in eachcomponent(equations)) - f_other = SVector{3, real(equations)}(f1, f2, f3) + w1 = gas_constant * v1 * rho_p + w2 = gas_constant * v2 * rho_p + w3 = gas_constant * (-rho_p) - return vcat(f_other, f_rho) + entrop_other = SVector{3, real(equations)}(w1, w2, w3) + + return vcat(entrop_other, entrop_rho) end -# Calculate 1D flux for a single point -@inline function flux(u, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u +""" +entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) - rho = density(u, equations) +Convert entropy variables to conservative variables +""" +@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) + @unpack gammas, gas_constants, cp, cv = equations + T = -1 / w[3] + v1 = w[1] * T + v2 = w[2] * T + v_squared = v1^2 + v2^2 + cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - + cv[i] * + (1 - log(T)) + + v_squared / + (2 * T)) / + gas_constants[i] - + 1) + for i in eachcomponent(equations)) - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + u2 = rho * v2 + gamma = help1 / help2 + u3 = p / (gamma - 1) + 0.5 * rho * v_squared + cons_other = SVector{3, real(equations)}(u1, u2, u3) + return vcat(cons_other, cons_rho) +end - f_rho = densities(u, v_normal, equations) - f1 = rho_v1 * v_normal + p * normal_direction[1] - f2 = rho_v2 * v_normal + p * normal_direction[2] - f3 = (rho_e + p) * v_normal +""" + total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) - f_other = SVector{3, real(equations)}(f1, f2, f3) +Calculate total entropy. +""" +@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations + T = temperature(u, equations) - return vcat(f_other, f_rho) + total_entropy = zero(u[1]) + for i in eachcomponent(equations) + total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) + end + + return total_entropy +end + +""" + density_gas_constant(u, equations::CompressibleEulerMulticomponentEquations2D{2}) + +Function that calculates overall density times overall gas constant. +""" +@inline function density_gas_constant(u, + equations::CompressibleEulerMulticomponentEquations2D) + @unpack gas_constants = equations + help = zero(u[1]) + for i in eachcomponent(equations) + help += u[i + 3] * gas_constants[i] + end + return help end """ - flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D) + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerMulticomponentEquations2D) Adaption of the entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) @@ -326,17 +298,12 @@ Adaption of the entropy conserving two-point flux by v_sum = v1_avg + v2_avg enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) @@ -371,27 +338,12 @@ Adaption of the entropy conserving two-point flux by return vcat(f_other, f_rho) end -""" - flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, - equations::CompressibleEulerMulticomponentEquations2D) - -Adaption of the entropy conserving and kinetic energy preserving two-point flux by -- Hendrik Ranocha (2018) - Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods - for Hyperbolic Balance Laws - [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743) -See also -- Hendrik Ranocha (2020) - Entropy Conserving and Kinetic Energy Preserving Numerical Methods for - the Euler Equations Using Summation-by-Parts Operators - [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42) -""" -@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) +@inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) # Unpack left and right state @unpack gammas, gas_constants, cv = equations - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) @@ -403,414 +355,45 @@ See also rho_ll = density(u_ll, equations) rho_rr = density(u_rr, equations) - # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) - inv_gamma_minus_one = 1 / (gamma - 1) - # extract velocities v1_ll = rho_v1_ll / rho_ll - v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) v2_ll = rho_v2_ll / rho_ll - v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) - - # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) - for i in eachcomponent(equations) - enth_ll += u_ll[i + 3] * gas_constants[i] - enth_rr += u_rr[i + 3] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] - end - - # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr - p_ll = T_ll * enth_ll - p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) - inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - - f_rho_sum = zero(T_rr) - if orientation == 1 - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg + p_avg - f2 = f_rho_sum * v2_avg - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) - else - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg - f2 = f_rho_sum * v2_avg + p_avg - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v2_rr + p_rr * v2_ll) - end - - # momentum and energy flux - f_other = SVector{3, real(equations)}(f1, f2, f3) - - return vcat(f_other, f_rho) -end - -@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - # Unpack left and right state - @unpack gammas, gas_constants, cv = equations - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], - u_rr[i + 3]) - for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + - u_rr[i + 3]) - for i in eachcomponent(equations)) - - # Iterating over all partial densities - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - - # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) - inv_gamma_minus_one = 1 / (gamma - 1) - - # extract velocities - v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_ll = rho_v2_ll / rho_ll v2_rr = rho_v2_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) - v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] - v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] - - # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) - for i in eachcomponent(equations) - enth_ll += u_ll[i + 3] * gas_constants[i] - enth_rr += u_rr[i + 3] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] - end - - # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr - p_ll = T_ll * enth_ll - p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) - inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - - f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * - (v_dot_n_ll + v_dot_n_rr) - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] - f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) - - # momentum and energy flux - f_other = SVector(f1, f2, f3) - - return vcat(f_other, f_rho) -end - -# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation -@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - - # Get the density and gas gamma - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - gamma_ll = totalgamma(u_ll, equations) - gamma_rr = totalgamma(u_rr, equations) - - # Get the velocities based on direction - if orientation == 1 - v_ll = rho_v1_ll / rho_ll - v_rr = rho_v1_rr / rho_rr - else # orientation == 2 - v_ll = rho_v2_ll / rho_ll - v_rr = rho_v2_rr / rho_rr - end - - # Compute the sound speeds on the left and right - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) - c_ll = sqrt(gamma_ll * p_ll / rho_ll) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) - c_rr = sqrt(gamma_rr * p_rr / rho_rr) - - λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) -end - -@inline function max_abs_speeds(u, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - v1 = rho_v1 / rho - v2 = rho_v2 / rho - - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2 + v2^2)) - c = sqrt(gamma * p / rho) - - return (abs(v1) + c, abs(v2) + c) -end - -@inline function rotate_to_x(u, normal_vector, - equations::CompressibleEulerMulticomponentEquations2D) - # cos and sin of the angle between the x-axis and the normalized normal_vector are - # the normalized vector's x and y coordinates respectively (see unit circle). - c = normal_vector[1] - s = normal_vector[2] - - # Apply the 2D rotation matrix with normal and tangent directions of the form - # [ n_1 n_2 0 0; - # t_1 t_2 0 0; - # 0 0 1 0 - # 0 0 0 1] - # where t_1 = -n_2 and t_2 = n_1 - - densities = @view u[4:end] - return SVector(c * u[1] + s * u[2], - -s * u[1] + c * u[2], - u[3], - densities...) -end - -# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction -# has been normalized prior to this back-rotation of the state vector -@inline function rotate_from_x(u, normal_vector, - equations::CompressibleEulerMulticomponentEquations2D) - # cos and sin of the angle between the x-axis and the normalized normal_vector are - # the normalized vector's x and y coordinates respectively (see unit circle). - c = normal_vector[1] - s = normal_vector[2] - - # Apply the 2D back-rotation matrix with normal and tangent directions of the form - # [ n_1 t_1 0 0; - # n_2 t_2 0 0; - # 0 0 1 0; - # 0 0 0 1 ] - # where t_1 = -n_2 and t_2 = n_1 - - densities = @view u[4:end] - return SVector(c * u[1] - s * u[2], - s * u[1] + c * u[2], - u[3], - densities...) -end - -# Convert conservative variables to primitive -@inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - prim_rho = SVector{ncomponents(equations), real(equations)}(u[i + 3] - for i in eachcomponent(equations)) - - rho = density(u, equations) - v1 = rho_v1 / rho - v2 = rho_v2 / rho - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) - prim_other = SVector{3, real(equations)}(v1, v2, p) - - return vcat(prim_other, prim_rho) -end - -# Convert conservative variables to entropy -@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - - # Multicomponent stuff - help1 = zero(rho) - gas_constant = zero(rho) - for i in eachcomponent(equations) - help1 += u[i + 3] * cv[i] - gas_constant += gas_constants[i] * (u[i + 3] / rho) - end - - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_square = v1^2 + v2^2 - gamma = totalgamma(u, equations) - - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) - s = log(p) - gamma * log(rho) - log(gas_constant) - rho_p = rho / p - T = (rho_e - 0.5 * rho * v_square) / (help1) - - entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * - (1 - log(T)) + - gas_constants[i] * - (1 + log(u[i + 3])) - - v_square / (2 * T)) - for i in eachcomponent(equations)) - - w1 = gas_constant * v1 * rho_p - w2 = gas_constant * v2 * rho_p - w3 = gas_constant * (-rho_p) - - entrop_other = SVector{3, real(equations)}(w1, w2, w3) - - return vcat(entrop_other, entrop_rho) -end - -# Convert entropy variables to conservative variables -@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) - @unpack gammas, gas_constants, cp, cv = equations - T = -1 / w[3] - v1 = w[1] * T - v2 = w[2] * T - v_squared = v1^2 + v2^2 - cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - - cv[i] * - (1 - log(T)) + - v_squared / - (2 * T)) / - gas_constants[i] - - 1) - for i in eachcomponent(equations)) - - rho = zero(cons_rho[1]) - help1 = zero(cons_rho[1]) - help2 = zero(cons_rho[1]) - p = zero(cons_rho[1]) - for i in eachcomponent(equations) - rho += cons_rho[i] - help1 += cons_rho[i] * cv[i] * gammas[i] - help2 += cons_rho[i] * cv[i] - p += cons_rho[i] * gas_constants[i] * T - end - u1 = rho * v1 - u2 = rho * v2 - gamma = help1 / help2 - u3 = p / (gamma - 1) + 0.5 * rho * v_squared - cons_other = SVector{3, real(equations)}(u1, u2, u3) - return vcat(cons_other, cons_rho) -end - -# Convert primitive to conservative variables -@inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas = equations - v1, v2, p = prim - - cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i + 3] - for i in eachcomponent(equations)) - rho = density(prim, equations) - gamma = totalgamma(prim, equations) - - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2) - - cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) - - return vcat(cons_other, cons_rho) -end - -@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - rho = density(u, equations) - T = temperature(u, equations) - - total_entropy = zero(u[1]) - for i in eachcomponent(equations) - total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) - end - - return total_entropy -end - -@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - help1 = zero(rho) + v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg + v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v_sum = v1_avg + v2_avg + enth = zero(v_sum) for i in eachcomponent(equations) - help1 += u[i + 3] * cv[i] + enth += rhok_avg[i] * gas_constants[i] end - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_square = v1^2 + v2^2 - T = (rho_e - 0.5 * rho * v_square) / help1 - - return T -end - -""" - totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) - -Function that calculates the total gamma out of all partial gammas using the -partial density fractions as well as the partial specific heats at constant volume. -""" -@inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas = equations + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) - help1 = zero(u[1]) - help2 = zero(u[1]) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v_dot_n_avg + for i in eachcomponent(equations)) + help1 = zero(T_ll) + help2 = zero(T_rr) for i in eachcomponent(equations) - help1 += u[i + 3] * cv[i] * gammas[i] - help2 += u[i + 3] * cv[i] + help1 += f_rho[i] * cv[i] + help2 += f_rho[i] end - return help1 / help2 -end - -@inline function density_pressure(u, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - gamma = totalgamma(u, equations) - rho_times_p = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) - - return rho_times_p -end - -@inline function density(u, equations::CompressibleEulerMulticomponentEquations2D) - rho = zero(u[1]) + f1 = (help2) * v1_avg + normal_direction[1] * enth / T + f2 = (help2) * v2_avg + normal_direction[2] * enth / T + f3 = ((help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + + v1_avg * f1 + v2_avg * f2) - for i in eachcomponent(equations) - rho += u[i + 3] - end - - return rho -end + f_other = SVector{3, real(equations)}(f1, f2, f3) -@inline function densities(u, v, equations::CompressibleEulerMulticomponentEquations2D) - return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v - for i in eachcomponent(equations)) + return vcat(f_other, f_rho) end end # @muladd diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl new file mode 100644 index 00000000000..1b80367ea10 --- /dev/null +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -0,0 +1,940 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function varnames(::typeof(cons2cons), + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + cons = ("rho_v1", "rho_v2", "rho_e") + rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) + return (cons..., rhos...) +end + +function varnames(::typeof(cons2prim), + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + prim = ("v1", "v2", "p") + rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) + return (prim..., rhos...) +end + +# Calculate 1D flux for a single point +@inline function flux(u, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = pressure(u, equations) + + if orientation == 1 + f_rho = u[4:end] .* v1 + f1 = rho_v1 * v1 + p + f2 = rho_v2 * v1 + f3 = (rho_e + p) * v1 + else + f_rho = u[4:end] .* v2 + f1 = rho_v1 * v2 + f2 = rho_v2 * v2 + p + f3 = (rho_e + p) * v2 + end + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +# Calculate 1D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + p = pressure(u, equations) + + f_rho = u[4:end] .* v_normal + f1 = rho_v1 * v_normal + p * normal_direction[1] + f2 = rho_v2 * v_normal + p * normal_direction[2] + f3 = (rho_e + p) * v_normal + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function rotate_to_x(u, normal_vector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ n_1 n_2 0 0; + # t_1 t_2 0 0; + # 0 0 1 0 + # 0 0 0 1] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] + s * u[2], + -s * u[1] + c * u[2], + u[3], + densities...) +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this back-rotation of the state vector +@inline function rotate_from_x(u, normal_vector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ n_1 t_1 0 0; + # n_2 t_2 0 0; + # 0 0 1 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] - s * u[2], + s * u[1] + c * u[2], + u[3], + densities...) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Determine the boundary numerical surface flux for a slip wall condition. +Imposes a zero normal velocity at the wall. +Density is taken from the internal solution state and pressure is computed as an +exact solution of a 1D Riemann problem. Further details about this boundary state +are available in the paper: +- J. J. W. van der Vegt and H. van der Ven (2002) + Slip flow boundary conditions in discontinuous Galerkin discretizations of + the Euler equations of gas dynamics + [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1) + +Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book +- Eleuterio F. Toro (2009) + Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + 3rd edition + [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + +Should be used together with [`UnstructuredMesh2D`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, surface_flux_function, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # rotate the internal solution state + u_local = rotate_to_x(u_inner, normal, equations) + + # compute the primitive variables + v_normal, v_tangent, p_local = cons2prim(u_local, equations) + + rho_local = density(u_local, equations) + gamma = totalgamma(u_inner, equations) + + # Get the solution of the pressure Riemann problem + if v_normal <= 0.0 + sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2 / ((gamma + 1) * rho_local) + B = p_local * (gamma - 1) / (gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + f_v = SVector{2, real(equations)}(p_star * normal_direction[1], + p_star * normal_direction[2]) + f_other = zeros(SVector{ncomponents(equations) + 1, real(equations)}) + + return vcat(f_v, f_other) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Should be used together with [`StructuredMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, + equations) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, + equations) + end + + return boundary_flux +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + prim_rho = SVector{ncomponents(equations), real(equations)}(u[i + 3] + for i in eachcomponent(equations)) + + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = pressure(u, equations) + prim_other = SVector{3, real(equations)}(v1, v2, p) + + return vcat(prim_other, prim_rho) +end + +# Convert primitive to conservative variables +@inline function prim2cons(prim, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + v1, v2, p = prim + + cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i + 3] + for i in eachcomponent(equations)) + rho = density(prim, equations) + gamma = totalgamma(prim, equations) + + rho_v1 = rho * v1 + rho_v2 = rho * v2 + + rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + + cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) + + return vcat(cons_other, cons_rho) +end + +""" + density(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates the overall density. +""" +@inline function density(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho = zero(u[1]) + for i in eachcomponent(equations) + rho += u[i + 3] + end + return rho +end + +""" + totalgamma(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates the total gamma out of all partial gammas. +This function has to be implemented for subtypes! +""" +function totalgamma end + +""" + temperature(u, equations::CompressibleEulerMulticomponentEquations2D) + +Calculate temperature. +This function has to be implemented for subtypes! +""" +function temperature end + +""" + pressure(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Calculates overall pressure. +""" +@inline function pressure(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + gamma = totalgamma(u, equations) + + p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + return p +end + +""" + density_pressure(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Calculates overall density times overall pressure. +""" +@inline function density_pressure(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + gamma = totalgamma(u, equations) + + rho_times_p = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) + return rho_times_p +end + +""" + density_gas_constant(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates overall density times overall gas constant. +This has to be implemented for subtypes! +""" +function density_gas_constant end + +""" + initial_condition_convergence_test(x, t, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function initial_condition_convergence_test(x, t, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + c = 2 + A = 0.1 + L = 2 + f = 1 / L + omega = 2 * pi * f + ini = c + A * sin(omega * (x[1] + x[2] - t)) + + v1 = 1.0 + v2 = 1.0 + + rho = ini + + # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + rho + for i in eachcomponent(equations)) + + prim1 = rho * v1 + prim2 = rho * v2 + prim3 = rho^2 + + prim_other = SVector{3, real(equations)}(prim1, prim2, prim3) + + return vcat(prim_other, prim_rho) +end + +""" + source_terms_convergence_test(u, x, t, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Same settings as in `initial_condition` + c = 2 + A = 0.1 + L = 2 + f = 1 / L + omega = 2 * pi * f + + gamma = totalgamma(u, equations) + + x1, x2 = x + si, co = sincos((x1 + x2 - t) * omega) + tmp1 = co * A * omega + tmp2 = si * A + tmp3 = gamma - 1 + tmp4 = (2 * c - 1) * tmp3 + tmp5 = (2 * tmp2 * gamma - 2 * tmp2 + tmp4 + 1) * tmp1 + tmp6 = tmp2 + c + + # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 + du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + tmp1 + for i in eachcomponent(equations)) + + du1 = tmp5 + du2 = tmp5 + du3 = 2 * ((tmp6 - 1.0) * tmp3 + tmp6 * gamma) * tmp1 + + du_other = SVector{3, real(equations)}(du1, du2, du3) + + return vcat(du_other, du_rho) +end + +""" + initial_condition_weak_blast_wave(x, t, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +A for multicomponent adapted weak blast wave taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_weak_blast_wave(x, t, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? + 2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + 1.0 : + 2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + 1.1691 + for i in eachcomponent(equations)) + + 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 + + prim_other = SVector{3, real(equations)}(v1, v2, p) + + return prim2cons(vcat(prim_other, prim_rho), equations) +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + + # Get the density and gas gamma + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + gamma_ll = totalgamma(u_ll, equations) + gamma_rr = totalgamma(u_rr, equations) + + # Get the velocities based on direction + if orientation == 1 + v_ll = rho_v1_ll / rho_ll + v_rr = rho_v1_rr / rho_rr + else # orientation == 2 + v_ll = rho_v2_ll / rho_ll + v_rr = rho_v2_rr / rho_rr + end + + # Compute the sound speeds on the left and right + p_ll = pressure(u_ll, equations) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + p_rr = pressure(u_rr, equations) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + + # Get the density and gas gamma + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + gamma_ll = totalgamma(u_ll, equations) + gamma_rr = totalgamma(u_rr, equations) + + # Calculate normal velocities and sound speed + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + v_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2]) + v_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2]) + + # Compute the sound speeds on the left and right + p_ll = pressure(u_ll, equations) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + p_rr = pressure(u_rr, equations) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speeds(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2 = u + + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + gamma = totalgamma(u, equations) + p = pressure(u, equations) + c = sqrt(gamma * p / rho) + + return (abs(v1) + c, abs(v2) + c) +end + +""" +flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Adaption of the entropy conserving and kinetic energy preserving two-point flux by +- Hendrik Ranocha (2018) +Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods +for Hyperbolic Balance Laws +[PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743) +See also +- Hendrik Ranocha (2020) +Entropy Conserving and Kinetic Energy Preserving Numerical Methods for +the Euler Equations Using Summation-by-Parts Operators +[Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42) +""" +@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # helpful variables + enth_ll = density_gas_constant(u_ll, equations) + enth_rr = density_gas_constant(u_rr, equations) + + # temperature and pressure + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + if orientation == 1 + f_rho = rhok_mean .* v1_avg + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg + f2 = f_rho_sum * v2_avg + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + else + f_rho = rhok_mean .* v2_avg + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + f2 = f_rho_sum * v2_avg + p_avg + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + end + + # momentum and energy flux + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # helpful variables + enth_ll = density_gas_constant(u_ll, equations) + enth_rr = density_gas_constant(u_rr, equations) + + # temperature and pressure + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + f_rho = rhok_mean .* (0.5 * (v_dot_n_ll + v_dot_n_rr)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] + f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +""" + flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +This flux is is a modification of the original kinetic energy preserving two-point flux by +- Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) + Kinetic energy and entropy preserving schemes for compressible flows + by split convective forms + [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058) + +The modification is in the energy flux to guarantee pressure equilibrium and was developed by +- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) + Preventing spurious pressure oscillations in split convective form discretizations for + compressible flows + [DOI: 10.1016/j.jcp.2020.110060](https://doi.org/10.1016/j.jcp.2020.110060) +""" +@inline function flux_shima_etal(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # Average each factor of products in flux + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + kin_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + f_rho = f_rho .* v1_avg + rho_avg_v = rho_avg * v1_avg + pv1_avg = 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + f1 = rho_avg_v * v1_avg + p_avg + f2 = rho_avg_v * v2_avg + f3 = p_avg * v1_avg * inv_gamma_minus_one + rho_avg_v * kin_avg + pv1_avg + else + f_rho = f_rho .* v2_avg + rho_avg_v = rho_avg * v2_avg + pv2_avg = 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + f1 = rho_avg_v * v1_avg + f2 = rho_avg_v * v2_avg + p_avg + f3 = p_avg * v2_avg * inv_gamma_minus_one + rho_avg_v * kin_avg + pv2_avg + end + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # Average each factor of products in flux + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr) + p_avg = 0.5 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on normal_direction + f_rho = f_rho .* v_dot_n_avg + rho_avg_v = rho_avg * v_dot_n_avg + f1 = rho_avg_v * v1_avg + p_avg * normal_direction[1] + f2 = rho_avg_v * v2_avg + p_avg * normal_direction[2] + f3 = (rho_avg_v * velocity_square_avg + + p_avg * v_dot_n_avg * inv_gamma_minus_one + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +""" + flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Kinetic energy preserving two-point flux by +- Kennedy and Gruber (2008) + Reduced aliasing formulations of the convective terms within the + Navier-Stokes equations for a compressible fluid + [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020) +""" +@inline function flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_e_ll = u_ll[3] + rho_e_rr = u_rr[3] + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Average each factor of products in flux + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + f_rho = f_rho .* v1_avg + f1 = rho_avg * v1_avg * v1_avg + p_avg + f2 = rho_avg * v1_avg * v2_avg + f3 = (rho_avg * e_avg + p_avg) * v1_avg + else + f_rho = f_rho .* v2_avg + f1 = rho_avg * v2_avg * v1_avg + f2 = rho_avg * v2_avg * v2_avg + p_avg + f3 = (rho_avg * e_avg + p_avg) * v2_avg + end + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_e_ll = u_ll[3] + rho_e_rr = u_rr[3] + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Average each factor of products in flux + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + p_avg = 0.5 * (p_ll + p_rr) + e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on normal_direction + f_rho = f_rho .* v_dot_n_avg + rho_avg_v = rho_avg * v_dot_n_avg + f1 = rho_avg_v * v1_avg + p_avg * normal_direction[1] + f2 = rho_avg_v * v2_avg + p_avg * normal_direction[2] + f3 = rho_avg_v * e_avg + p_avg * v_dot_n_avg + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3 = v * u_ll + f3 = f3 + p_ll * v + # density fluxes + f_rho = SVector{ncomponents(equations), + real(equations)}(u_ll[i + 3] * v + for i in eachcomponent(equations)) + else + f1, f2, f3 = v * u_rr + f3 = f3 + p_rr * v + # density fluxes + f_rho = SVector{ncomponents(equations), + real(equations)}(u_rr[i + 3] * v + for i in eachcomponent(equations)) + end + + if orientation == 1 + f1 = f1 + p + else # orientation == 2 + f2 = f2 + p + end + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Note that this is the same as computing v_ll and v_rr with a normalized normal vector + # and then multiplying v by `norm_` again, but this version is slightly faster. + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + # We treat the energy term analogous to the potential temperature term in the paper + # by Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3 = u_ll * v + f3 = f3 + p_ll * v + # density fluxes + f_rho = SVector{ncomponents(equations), + real(equations)}(u_ll[i + 3] * v + for i in eachcomponent(equations)) + else + f1, f2, f3 = u_rr * v + f3 = f3 + p_rr * v + # density fluxes + f_rho = SVector{ncomponents(equations), + real(equations)}(u_rr[i + 3] * v + for i in eachcomponent(equations)) + end + f1 = f1 + p * normal_direction[1] + f2 = f2 + p * normal_direction[2] + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end +end # @muladd diff --git a/src/equations/compressible_euler_multicomponent_moist_2d.jl b/src/equations/compressible_euler_multicomponent_moist_2d.jl new file mode 100644 index 00000000000..1e3e2b28d9f --- /dev/null +++ b/src/equations/compressible_euler_multicomponent_moist_2d.jl @@ -0,0 +1,401 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" +CompressibleMoistEulerEquations2D(;gammas, gas_constants, c_p, c_v, L_00) + +Similar to `CompressibleEulerMulticomponentEquations2D` but including latent heat of +vaporization +```math +... +``` + +Either the specific heat ratios `gammas` and the gas constants `gas_constants` in should be +passed as tuples, or the specific heats at constant volume `cv` and the specific heats at +constant pressure `cp`. The other quantities are then calculated considering a calorically +perfect gas. +""" +struct CompressibleMoistEulerEquations2D{NVARS, NCOMP, RealT <: Real} <: + AbstractCompressibleEulerMulticomponentEquations{2, NVARS, NCOMP} + gammas::SVector{NCOMP, RealT} + gas_constants::SVector{NCOMP, RealT} + c_p::SVector{NCOMP, RealT} + c_v::SVector{NCOMP, RealT} + L_00::RealT # latent heat of evaporation at 0 K + + function CompressibleMoistEulerEquations2D{NVARS, NCOMP, RealT}(gammas::SVector{NCOMP, + RealT}, + gas_constants::SVector{NCOMP, + RealT}, + c_p::SVector{NCOMP, + RealT}, + c_v::SVector{NCOMP, + RealT}, + L_00::RealT) where { + NVARS, + NCOMP, + RealT <: + Real + } + new(gammas, gas_constants, c_p, c_v, L_00) + end +end + +function CompressibleMoistEulerEquations2D(; + gammas = nothing, gas_constants = nothing, + c_p = nothing, c_v = nothing, L_00) + if gammas !== nothing && gas_constants !== nothing + _gammas = promote(gammas...) + _gas_constants = promote(gas_constants...) + RealT = promote_type(eltype(_gammas), eltype(_gas_constants), + typeof(gas_constants[1] / (gammas[1] - 1))) + _c_v = _gas_constants ./ (_gammas .- 1) + _c_p = _gas_constants + _gas_constants ./ (_gammas .- 1) + elseif c_p !== nothing && c_v !== nothing + _c_p = promote(c_p...) + _c_v = promote(c_v...) + RealT = promote_type(eltype(_c_p), eltype(_c_v), typeof(c_p[1] / c_v[1])) + _gas_constants = _c_p .- _c_v + _gammas = _c_p ./ _c_v + else + throw(DimensionMismatch("Either `gammas` and `gas_constants` or `c_p` and `c_v` \ + have to be filled with at least one value")) + end + + NVARS = length(_gammas) + 3 + NCOMP = length(_gammas) + + __gammas = SVector(map(RealT, _gammas)) + __gas_constants = SVector(map(RealT, _gas_constants)) + __c_p = SVector(map(RealT, _c_p)) + __c_v = SVector(map(RealT, _c_v)) + + return CompressibleMoistEulerEquations2D{NVARS, NCOMP, RealT}(__gammas, + __gas_constants, + __c_p, __c_v, L_00) +end + +@inline function Base.real(::CompressibleMoistEulerEquations2D{NVARS, NCOMP, + RealT}) where {NVARS, + NCOMP, + RealT} + RealT +end + +""" + totalgamma(u, equations::CompressibleMoistEulerEquations2D) + +Function that calculates the total gamma out of all partial gammas using the +partial density fractions as well as the partial specific heats at constant volume. +""" +@inline function totalgamma(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_v, c_p = equations + + help1 = zero(u[1]) + help2 = zero(u[1]) + + # compute weighted averages of cp and cv + # normalization by total rho not required, would cancel below + for i in eachcomponent(equations) + help1 += u[i + 3] * c_p[i] + help2 += u[i + 3] * c_v[i] + end + + return help1 / help2 +end + +""" +pressure(u, equations::CompressibleMoistEulerEquations2D) + +Calculate pressure. This differs from the calculation in `AbstractCompressibleEulerMulticomponentEquations` in that the latent heat is accounted for. +""" +@inline function pressure(u, equations::CompressibleMoistEulerEquations2D) + @unpack L_00 = equations + rho_v1, rho_v2, rho_e, rho_d, rho_v = u + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + gamma = totalgamma(u, equations) + p = (gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - L_00 * rho_v) + return p +end + +""" + density_pressure(u, equations::CompressibleMoistEulerEquations2D) + +Calculates overall density times overall pressure. +""" +@inline function density_pressure(u, + equations::CompressibleMoistEulerEquations2D) + rho = density(u, equations) + p = pressure(u, equations) + return rho * p +end + +""" +temperature(u, equations::CompressibleMoistEulerEquations2D) + +Calculate temperature. Account for latent heat. +""" +@inline function temperature(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_v, L_00 = equations + rho_v1, rho_v2, rho_e, rho_d, rho_v = u + + rho = density(u, equations) + help1 = zero(rho) + + # compute weighted average of cv + # normalization by rho not required, cancels below + for i in eachcomponent(equations) + help1 += u[i + 3] * c_v[i] + end + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + T = (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - L_00 * rho_v) / help1 + return T +end + +""" + density_gas_constant(u, equations::CompressibleMoistEulerEquations2D) + +Function that calculates overall density times overall gas constant. +""" +@inline function density_gas_constant(u, + equations::CompressibleMoistEulerEquations2D) + @unpack gas_constants = equations + help = zero(u[1]) + for i in eachcomponent(equations) + help += u[i + 3] * gas_constants[i] + end + return help +end + +""" +cons2entropy(u, equations::CompressibleMoistEulerEquations2D) + +Convert conservative variables to entropy. +""" +@inline function cons2entropy(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_v, gas_constants = equations + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + # Multicomponent stuff + help1 = zero(rho) + gas_constant = zero(rho) + for i in eachcomponent(equations) + help1 += u[i + 3] * c_v[i] + gas_constant += gas_constants[i] * (u[i + 3] / rho) + end + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + p = pressure(u, equations) + rho_p = rho / p + T = temperature(u, equations) + + entrop_rho = SVector{ncomponents(equations), real(equations)}((c_v[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 3])) - + v_square / (2 * T)) + for i in eachcomponent(equations)) + + w1 = gas_constant * v1 * rho_p + w2 = gas_constant * v2 * rho_p + w3 = gas_constant * (-rho_p) + + entrop_other = SVector{3, real(equations)}(w1, w2, w3) + + return vcat(entrop_other, entrop_rho) +end + +""" +entropy2cons(w, equations::CompressibleMoistEulerEquations2D) + +Convert entropy variables to conservative variables +""" +@inline function entropy2cons(w, equations::CompressibleMoistEulerEquations2D) + @unpack gammas, gas_constants, cv = equations + T = -1 / w[3] + v1 = w[1] * T + v2 = w[2] * T + v_squared = v1^2 + v2^2 + cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - + cv[i] * + (1 - log(T)) + + v_squared / + (2 * T)) / + gas_constants[i] - + 1) + for i in eachcomponent(equations)) + + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + u2 = rho * v2 + gamma = help1 / help2 + u3 = p / (gamma - 1) + 0.5 * rho * v_squared + L_00 * cons_rho[2] + cons_other = SVector{3, real(equations)}(u1, u2, u3) + return vcat(cons_other, cons_rho) +end + +""" + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleMoistEulerEquations2D) + +Adaption of the entropy conserving two-point flux by +- Ayoub Gouasmi, Karthik Duraisamy (2020) + "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" + [arXiv:1904.00972v3](https://arxiv.org/abs/1904.00972) [math.NA] 4 Feb 2020 +""" +@inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, + equations::CompressibleMoistEulerEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, cv, L_00 = equations + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v_sum = v1_avg + v2_avg + + enth = zero(v_sum) + for i in eachcomponent(equations) + enth += rhok_avg[i] * gas_constants[i] + end + + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + + # Calculate fluxes depending on orientation + help1 = zero(T_ll) + help2 = zero(T_rr) + if orientation == 1 + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + help1 += f_rho[i] * cv[i] + help2 += f_rho[i] + end + f1 = (help2) * v1_avg + enth / T + f2 = (help2) * v2_avg + # Account for latent heat. This is the only difference compared to + # AbstractCompressibleEulerMulticomponentEquations + f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + v2_avg * f2 + L_00 * f_rho[2] + else + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + help1 += f_rho[i] * cv[i] + help2 += f_rho[i] + end + f1 = (help2) * v1_avg + f2 = (help2) * v2_avg + enth / T + # Account for latent heat. This is the only difference compared to + # AbstractCompressibleEulerMulticomponentEquations + f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + v2_avg * f2 + L_00 * f_rho[2] + end + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, c_v, L_00 = equations + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg + v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v_sum = v1_avg + v2_avg + + enth = zero(v_sum) + for i in eachcomponent(equations) + enth += rhok_avg[i] * gas_constants[i] + end + + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v_dot_n_avg + for i in eachcomponent(equations)) + + help1 = zero(T_ll) + help2 = zero(T_rr) + for i in eachcomponent(equations) + help1 += f_rho[i] * c_v[i] + help2 += f_rho[i] + end + + f1 = (help2) * v1_avg + normal_direction[1] * enth / T + f2 = (help2) * v2_avg + normal_direction[2] * enth / T + + # Account for latent heat. This is the only difference compared to + # AbstractCompressibleEulerMulticomponentEquations + f3 = ((help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + + v1_avg * f1 + v2_avg * f2 + L_00 * f_rho[2]) + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8f476cf6f16..e6e95c1a182 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -418,24 +418,15 @@ include("compressible_euler_quasi_1d.jl") # CompressibleEulerMulticomponentEquations abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: AbstractEquations{NDIMS, NVARS} end -include("compressible_euler_multicomponent_1d.jl") -include("compressible_euler_multicomponent_2d.jl") - -# PolytropicEulerEquations -abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <: - AbstractEquations{NDIMS, NVARS} end -include("polytropic_euler_2d.jl") # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, - NCOMP}) where { - NDIMS, - NVARS, - NCOMP - } + NCOMP}) where + {NDIMS, NVARS, NCOMP} NCOMP end + """ eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) @@ -447,6 +438,16 @@ In particular, not the components themselves are returned. Base.OneTo(ncomponents(equations)) end +include("compressible_euler_multicomponent_1d.jl") +include("compressible_euler_multicomponent_abstract_2d.jl") +include("compressible_euler_multicomponent_2d.jl") +include("compressible_euler_multicomponent_moist_2d.jl") + +# PolytropicEulerEquations +abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end +include("polytropic_euler_2d.jl") + # Ideal MHD abstract type AbstractIdealGlmMhdEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end @@ -462,11 +463,8 @@ include("ideal_glm_mhd_multicomponent_2d.jl") # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractIdealGlmMhdMulticomponentEquations{NDIMS, NVARS, - NCOMP}) where { - NDIMS, - NVARS, - NCOMP - } + NCOMP}) where + {NDIMS, NVARS, NCOMP} NCOMP end """ diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 64a1faf05b8..fda831fe15b 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -291,6 +291,36 @@ end end end +#= +@trixi_testset "elixir_eulermultimoist_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermultimoist_warm_bubble.jl"), + l2=[ + 1.5123651627525257e-5, + 1.51236516273878e-5, + 2.4544918394022538e-5, + 5.904791661362391e-6, + 1.1809583322724782e-5, + ], + linf=[ + 8.393471747591974e-5, + 8.393471748258108e-5, + 0.00015028562494778797, + 3.504466610437795e-5, + 7.00893322087559e-5, + ], + cells_per_dimension=(32, 16), + tspan=(0.0, 10.0)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end +=# + @trixi_testset "elixir_euler_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), # Expected errors are exactly the same as with TreeMesh! diff --git a/test/test_unit.jl b/test/test_unit.jl index 03a78f6918a..593b9be6b06 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -619,14 +619,17 @@ end end @timed_testset "Consistency check for single point flux: CEMCE" begin - equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), - gas_constants = (0.4, 0.4)) + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.2), + gas_constants = (0.4, 0.6)) u = SVector(0.1, -0.5, 1.0, 1.0, 2.0) orientations = [1, 2] - for orientation in orientations - @test flux(u, orientation, equations) ≈ - flux_ranocha(u, u, orientation, equations) + fluxes = [flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340)] + for f_test in fluxes + for orientation in orientations + @test flux(u, orientation, equations) ≈ + f_test(u, u, orientation, equations) + end end end @@ -1400,25 +1403,26 @@ end @testset "FluxRotated vs. direct implementation" begin @timed_testset "CompressibleEulerMulticomponentEquations2D" begin - equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.2), gas_constants = (0.4, - 0.4)) + 0.6)) normal_directions = [SVector(1.0, 0.0), SVector(0.0, 1.0), SVector(0.5, -0.5), SVector(-1.2, 0.3)] u_values = [SVector(0.1, -0.5, 1.0, 1.0, 2.0), SVector(-0.1, -0.3, 1.2, 1.3, 1.4)] + fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, + FluxLMARS(340)] - f_std = flux - f_rot = FluxRotated(f_std) - println(typeof(f_std)) - println(typeof(f_rot)) - for u in u_values, - normal_direction in normal_directions + for f_std in fluxes + f_rot = FluxRotated(f_std) + for u_ll in u_values, u_rr in u_values, + normal_direction in normal_directions - @test f_rot(u, normal_direction, equations) ≈ - f_std(u, normal_direction, equations) + @test f_rot(u_ll, u_rr, normal_direction, equations) ≈ + f_std(u_ll, u_rr, normal_direction, equations) + end end end